20

In one dimension, you can map your infinite interval to a finite interval using integration by substitution, e.g. $$\int_a^b f(x)\,\mathrm dx \quad=\quad \int_{u^{-1}(a)}^{u^{-1}(b)}f(u(t))u'(t)\,\mathrm dt$$ Where $u(x)$ is some function that goes off to infinity in some finite range, e.g. $\tan(x)$: $$\int_{-\infty}^{\infty}f(x)\,\mathrm dx \quad=\quad ...


15

There is a straightforward solution with only two passes through the data: First compute $$K := \max_i\; a_i,$$ which tells you that, if there are $n$ terms, then $$\sum_i e^{a_i} \le n e^K.$$ Since you presumably don't have $n$ anywhere near as large as even $10^{20}$, you should have no worry about overflowing in the computation of $$\tau := \sum_i e^{...


11

The standard way of doing it is to extract from the expression for $f(x)$ an exponential prefactor, transform that to $e^{-x^2}$, and then use Gaussian quadrature rules (or Gauss Kronrod) with this as a weight. If $f$ is smooth, this usually gives excellent results. In $R^3$, the same works with weight $e^{-|x|^2}$, and appropriate cubature formulas can be ...


10

To keep precision while you add doubles together you need to use Kahan Summation, this is the software equivalent to having a carry register. This is fine for most values, but if you are getting overflow then you are hitting the limits of IEEE 754 double-precision which would be about $e^{709.783}$. At this point you need a new representation. You can ...


10

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


9

The long thermalization time that you're running into is a generic problem that typically goes under the name "critical slowing down" and is common to the local-update scheme that you're using (you update by locally changing a single spin at a time). Once you realize that, the way out is to do better sampling - local updates are out so you have to invent ...


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


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

For one-dimensional quadrature, you can check the book on Quadpack (a golden oldie but still very relevant in one-dimensional quadrature) and the techniques used in the algorithm QAGI, an automatic integrator for an infinite range. Another technique is the double-exponential quadrature formula, nicely implemented for an infinite interval by Ooura. For ...


7

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

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

In principle one could preselect the number of samples to be drawn from each sub-distribution, then visit each sub distribution only once and draw than number of points. That is Find the random set $<n_1, n_2, \dots, n_k>$ such that $n = \sum_{i=1}^k n_i$ and respecting the weights. I believe that you do this by drawing a Poisson distribution a ...


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


6

If you remember how quadrature worked, then you'll also know one way how to approximate infinite integrals. Namely: for quadrature, you approximate the function $f(x)$ you want to integrate by something similar, say a polynomial $\tilde f(x)$ (or a piecewise polynomial) for which you can write down the integral analytically. You get $\tilde f$ from $f$ by ...


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


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

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

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

If you want to use Monte Carlo integration, you could start by using importance sampling with a sampler that roughly approximates your integrand. The better your sampler matches your integrand, the less variance in your integral estimates. It doesn't matter than your domain is infinite as long as your sampler has the same domain.


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

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


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