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


13

A more structured way of providing a basis or quadrature (which may replace MC in many cases) in multiple dimensions is that of sparse grids, which combines some family of one dimensional rules of varying order in such a way as to have merely exponential growth in dimension, $2^d$, rather than having it be that dimension is an exponent of the resolution $N^d$...


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

I would warmly recommend to use StochKit, which is really extremely reliable and efficient. It implements the Gillespie's direct method as well as its most popular optimized variants. Furthermore, it has a built-in mechanism for selecting the appropriate simulation scheme depending on your network. However, it is written in C++. If you really need to ...


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

(Disclaimer: I haven't used any of the packages or code below.) Finding any implementation of Gillespie's method in Python was a bit of a challenge; the most fruitful search terms seemed to be "tau leap" or "kinetic Monte Carlo." This blog post implements Gillespie's algorithm, though it's not clear that it's efficient. One of the commenters mentions some ...


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

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

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

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

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


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

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

For the Ising model, the relevant quantities are properties such as the correlation length, the magnetic susceptibility, and the heat capacity. What makes finding these exponents challenging is that the calculations have to be performed in the neighborhood of the critical point: you need to be close enough to the critical temperature so that you can see the ...


6

As a general rule, it's easy to understand why regular grids can't go much beyond 3 or 4-dimensional problems: in d dimensions, if you want to have a minimum of N points per coordinate direction, you'll get N^d points overall. Even for relatively nice functions in 1d, you need at least N=10 grid points to resolve them at all, so the overall number of points ...


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

As you may already know if you're asking this question, Gillespie's method is often used in computational systems biology. There are a number of software system that provide implementations; you may want to browse the SBML Software Guide to find some of them (specifically, the ones that are known to support SBML, something that I realize may not be relevant ...


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

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

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


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