11

This is the stochastic root-finding problem, as in The stochastic root-finding problem: Overview, solutions, and open questions.


9

Source-to-source transformation is considered the gold standard in terms of performance. OO approaches seem to be almost as good, in that there are more OO packages out there, and performance is not mentioned as a significant drawback. If you find an OO library you like for the language you're working in, I'd use it first, and then figure out later if you ...


9

Like you say, using the Mersenne Twister for parallel computations is almost always done incorrectly, as the correct method is tricky to implement. By far the easiest and best answer would be to move away from the Mersenne Twister entirely, and use something like the PCG family, which provides multiple streams out of the box. The Mersenne Twister is known ...


7

There are several Bayesian optimization techniques you could try. Easiest are based on Gaussian process: Harold J. Kushner. A new method of locating the maximum of an arbitrary multipeak curve in the presence of noise. Journal of Basic Engineering, pages 86:97–106, March 1964. J. Mockus. The Bayesian approach to global optimization. Lecture Notes in Control ...


7

Our Matlab package SnobFit was created precisely for this purpose. No assumption about the distribution of the noise is needed. Moreover, function values can be supplied through text files, thus you can apply it to functions implemented in any system able to write a text file. See http://www.mat.univie.ac.at/~neum/software/snobfit/ SnobFit had been ...


7

Disclaimer: I wrote my PhD thesis on adaptive quadrature, so this answer will be severely biased towards my own work. GSL's QAGS is the old QUADPACK integrator, and it is not entirely robust, especially in the presence of singularities. This usually leads to users requesting far more digits of accuracy than they actually need, thus making the integration ...


7

Your lattice consists of 5 x 5 x 5 = 125 spins, so your number of Montecarlo steps to reach equilibrium should be >> 125, because you randomly picking a site and flipping it, so random numbers should uniformly generated so that it will cover whole lattice. For much finer measurement of thermodynamic quantities, you should take more number of points between ...


6

If I understood your question correctly: Subdivide the interval $[0,1]$ in $N$ segments, each having a width proportional to your probabilities $p_i$ (where $\sum_{i=1}^N p_i = 1$). Use a pseudo-random-number-generator using the uniform probability density distribution on $[0,1]$ to generate the number $x$. Determine in which segment $x$ is located; this ...


6

As part of the lab the developed and maintains StochKit. I am happy to hear that it is highly recommended in the previous answers. However, I wanted to update everyone. There already is a python wrapper for StochKit: GillesPy. I would also recommend you checkout StochSS, This is an fully functional modeling and simulation IDE that uses StochKit and as ...


5

Monte Carlo methods can in general not compete with adaptive quadrature unless you have a high dimensional integral where you cannot afford the combinatorial explosion of quadrature points with the dimension. The reason is relatively easy to understand. Take, for example, just $\int_{[0,1]^n} f(x)\; d^nx$ where $n$ is the dimension of the problem. Let's say,...


5

The key is to take the differences $\Delta x$, $\Delta y$, and $\Delta z$ separately before beginning. Given the edge vectors ${\bf a}$, ${\bf b}$, and ${\bf c}$ that define a unit cell that has one corner at the origin, then your tilt factors are $b_1$, $c_1$, and $c_2$, where $n_m$ defines the $m^{th}$ component of vector ${\bf n}$. Now, for each ...


5

Some of my research involves solving large scale stochastic partial differential equations. In which case, traditional monte carlo approximation of the integral of interest converges too slowly for it to be worthwhile in a practical sense... i.e. I don't want to have to run 100 times more simulations just to get a decimal point more accuracy to the integral....


5

A fantastic introductory book (final year undergraduate to graduate level) for simulation of electron and holes in semiconductor is, Fundamentals of Carrier Transport by M. Lundstrom. The book covers difference techniques of solving the Boltzmann transport equation and it has a chapter on Monte-Carlo. It should provide you with an overview of the relevant ...


5

You haven't specified the distribution of $x(t)$. I'll assume that you want to use a complex normal distribution, since that choice makes it reasonably easy to solve the problem and because this assumption is quite common in signal processing. I'll also discretize the problem so that you're generating a vector $X$ of $N$ entries with a specified complex ...


5

TLDR Use Python to manage/modify your input and coral your output, and use HDF5 to organize/store your data. As complex as it might seem at first it'll still be simpler than SQL-anything. Longer answer + Example I personally use a combination of Python scripting and the HDF5 file format to deal with these kinds of situations. Python scripting can handle the ...


5

You definitely want to derandomize your program during development. Otherwise you will not be able to debug it since problems are not reproducible. At the same time, once you know the algorithm is working, you need to run it for multiple seeds or with different random number generators to ensure that your results (such as ensemble averages, standard ...


5

This is possible (see [1]) but uncommon, as it requires Monte Carlo moves that alter the current conformations by a very small perturbation. In that setting of "small" Metropolis MC moves, it is usually easier (both in theory and in practice) to just use Molecular Dynamics instead. [1] Kikuchi, K., Yoshida, M., Maekawa, T., & Watanabe, H. (1991). ...


5

It is in some cases possible to map the dynamics obtained in MC simulations to other (more realistic) dynamics, especially for the case of dense colloidal suspensions. The following two papers talk about the problems and caveats of performing such a mapping: http://aip.scitation.org/doi/10.1063/1.3414827 (spherical particles) http://aip.scitation.org/doi/10....


5

I sincerely thank @Daniel Shapero for directing me towards this answer. Discontinuity in the specific heat or susceptibility curves to be visible significantly, you should take much more finer measurements for large number of sweeps, say, I ran the simulation for 1024 steps for system to reach equilibrium 1024 steps for sparse averaging/to measure specific ...


4

I would use the first option and would use a synchronous AC run before (using the GPU), to detect collisions, execute a step of a hexagonal AC whose rule is the value of the center cell = Sum (neighbors), This CA must have seven states should be initiated with randomly selected cell, and their status verified before running the update rule for each gpu. ...


4

Advantages of traditional Monte-Carlo integration over quasi-Monte Carlo integration are discussed in Kocis and Whiten's paper here. They list the following reasons: The error bound of qmc methods $\mathcal{O}(\log(N)^{d}/N)$ is "theoretically" better than the $\mathcal{O}(N^{-1/2})$ bound given by naive Monte-Carlo. However, for values of $N$ that are ...


4

Monte Carlo simulations are the method of choice for computation of electron scattering. Tricks like importance sampling are used sometimes, so you might say it's not plain old Monte Carlo. But the main point is probably that an inherently stochastic process is simulated here, while you are asking only about using Monte Carlo for integration. Because nobody ...


4

James Spall's SPSA algorithm (short for Stochastic Perturbation Simulated Annealing, if I recall correctly) has been designed for exactly this kind of problem. He has a couple of papers where he uses it for problems like the one you describe.


4

If you really want the variables to be called pi1, pi2, ..., you can build the variable names using eval: for k=1:1:4 ... eval(['pi' num2str(k) ' = 4*N2/(N(k)*R.^2)']) ... end But a more typical thing to do would be to use an array: for k=1:1:4 ... pi(k) = 4*N2/(N(k)*R.^2) ... end If you're not planning on using the names of the pi variables ...


4

I'm in the middle of doing this for myself. I think the most appropriate analogue to the Gaussian would be the heat kernel in hyperbolic space. Fortunately, this has been figured out before: https://www.math.uni-bielefeld.de/~grigor/nog.pdf (also available in a Bulletin of the London Mathematical Society). If you use the standard decay ($e^{-dist^2/constant}...


4

As noted in the comments by Kirill, the y-axes of the two plots are very different. And if they are rescaled accordingly, the boxes will certainly look very similar, if not identical. Therefore, it is very reasonable to conclude that the raw simulation result in data.dat coming from your C++ code is correct, no matter what seed has been used for random-...


4

What you're looking for goes under the name of quasi-Monte Carlo (QMC) sequences. Quasi Monte Carlo sequences are "more random than random", i.e. they fill high dimensional spaces better than random sequences tend to do, which somewhat matches a more intuitive form of randomness. But more importantly, since they "fill space well" (a ...


4

As I understand, your ultimate goal is to solve an inverse problem (i.e., infer some parameters from given data / observations). To this end, you want to apply Bayesian Inference, which relates the posterior (i.e., the probability distribution of the unknown parameters) to the likelihood (i.e., the probability model of observing some values given the ...


4

Here is a corrected and slightly improved code, set up here for calculation of the full torus volume to verify the result. import numpy as np from scipy import random def func(x,y,z): return x**2 + (np.sqrt(y**2 + z**2)-5)**2 def MCvolume(N=1000): #Limits of Integration x_min = -2 x_max = 2 y_min = -7 y_max = 7 #1 z_min = -...


3

Here is a very good introductory book chapter on Monte Carlo simulations of classical spin systems, which also introduces and discusses the critical exponents: W. Janke, Monte Carlo Methods in Classical Statistical Physics, invited lectures, in: Computational Many-Particle Physics, Wilhelm & Else Heraeus Summerschool, Greifswald, Germany, 18-29 ...


Only top voted, non community-wiki answers of a minimum length are eligible