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

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

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

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

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

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

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

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

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

If your PDF is bounded, you could try approximating its inverse with a high-degree polynomial interpolant. This is usually considered a bad thing, but that's just a myth. Some things to keep in mind: Instead of using an equispaced grid, interpolate at the Chebyshev nodes of the first kind, i.e. $x_i = \cos\left(\pi\frac{2i-1}{2N}\right)$, for $F(x)$ ...


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


3

Well I guess that this question could take place in both computational science and physics, depending on the kind of physicist or computational scientist stubbling upon it. Yet for your answer, as I saw some solid MC algorithm, I advise you to put your rejection goal (or acceptation) as a constant inside your code and vary your trial size proportionnaly to ...


3

I found that the grid summing method beat out both the Cubature package and the adaptive Monte Carlo schemes from the GSL package (standard, VEGAS, and MISER) when it came to both speed an accuracy (in fact, Cubature's hcubature adaptive methods wouldn't converge for abstol<1e-1, showing that it really cannot handle these kinds of problems). For future ...


3

If you want to use MT, you can use SFMT as your PRNG and SFMT jump to generate multiple streams. You can simply initialise MT with one seed, and then jump ahead by e.g. $1 \cdot 10^{60}$, $2 \cdot 10^{60}$, $3 \cdot 10^{60}$ … steps to generate multiple streams. Jumping is somewhat expensive, but you only need to do it once when you initialise your PRNGs.


3

If I understand what information you have available, what you want is not possible: the information available to you is not enough to determine the entropy. It's not even enough to approximate the entropy. It sounds like you have a way to sample from the distribution $p(\cdot)$, and you have a way to compute the ratio $p(a_1)/p(a_2)$ for any pair of ...


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