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

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

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

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 ...


3

As far as I can tell, indexing a simple cubic lattice in the way you require should be easy. If we let $i, j, k$ be the indices along the $x, y, z$ directions, we can easily map these to a linear index $n$ like so: $$ n = ((k N_y) + j) N_x + i$$ Next, we need to accomodate 8 atoms per unit cell. The easiest way to do this should be to add a stride of 8 to ...


3

Suppose, that your circle have unitary radius: $r=1$, then length of square's side equal to $l = r_a \cdot d = 2\cdot r_a$. So, area of circle equals to $S_c = \pi$ and area of square equals to $S_s = 4 r_a^2$. If you "throw" $N_t$ random points inside of square, a part of them (say, $N_c$) will fall into the circle. Relation $\frac{N_c}{N_t}$ equals to ...


3

This is discussed Section 11.1 of Kloeden and Platen, "Numerical Solution of Stochastic Differential Equations". There it states: Using the deterministic Taylor expansion it is easy to show that the ratio $$\frac{1}{\sqrt\Delta}\left\{b(\tau_n,Y_n + a\Delta + b\sqrt\Delta) - b(\tau_n,Y_n)\right\}$$ is a forward difference approximation for $b\frac{\...


3

Quasi random numbers are not statistically independent so if your algorithm requires this, you should not use them. One area where you can usually use quasi random numbers instead of pseudo random numbers is in the evaluation of integrals via monte carlo integration. In this application, the fact that quasi random numbers are not truly random is not an ...


3

The main purpose of sampling tons of pseudo-random data as opposed to non-random data is related to Runge's phenomenon for polynomial interpolation: Uniform spacing of interpolation points is often a bad idea. But choosing better interpolation points require knowledge of the function you want to interpolate (or integrate etc.). If you don't have that ...


3

Arrays allocated with brackets such as int s[i][j] are fixed to be of size i x j, whereas arrays allocated with new and int *s / int **s are not necessarily uniformly deep and can be resized with subsequent calls to new. If you're solving on a fixed grid, there isn't any need for being able to resize the array, so the two options above will have similar ...


3

You may not be able to determine the exact number of points required to obtain 1% error, but you can estimate the order of magnitude of points needed to obtain this accuracy. Monte Carlo converges at a rate $O(\frac{1}{\sqrt{N}})$, where $N$ is the sample size. This means the absolute error is bounded as $|\mu-\mu_{approx}|<\frac{C}{\sqrt{N}}$, where ...


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