I implemented as an exercise a program to sample the statistics of the escape time of a Brownian particle in a potential well. I used the Euler-Maruyama method to numerically integrate the trajectories of the solution of the Ornstein-Uhlenbeck equation. For each trajectory I stored the escape time from the well.

#include <iostream>
#include <fstream>
#include <cmath>
#include <random>
#include <array>
#include <iterator>
#include <memory>
#include <omp.h>

typedef std::mt19937_64 rng_t;
typedef double real_t;

namespace params {
  constexpr real_t t0 = 0.0;
  constexpr real_t x0 = 0.0;
  constexpr real_t v0 = 0.0;
  constexpr real_t tau = 0.005;

  constexpr real_t beta = 1.0;
  constexpr real_t alpha = 1.0;
  constexpr int iters = 10000;

  constexpr real_t well_width = 1.0;
  constexpr real_t well_depth = 1.0;
  constexpr real_t ramp_width = 0.1;

using namespace std;
using namespace params;

constexpr auto half_width = well_width / 2;
constexpr auto g = well_depth / ramp_width;

template<typename T> T gamma(T x) {
  if (x > -half_width && x < -half_width + ramp_width)    return g;
  else if (x > half_width - ramp_width && x < half_width) return -g;
  else                                                    return 0.0;

template<typename T> class OrnsteinUhlenbeck {
    T EscapeTime() { while (IsInBounds()) EMTimestep(); return t; };
    OrnsteinUhlenbeck(shared_ptr<rng_t> r) { rng = r; };
    void EMTimestep(){
      normal_distribution<T> dW(0.0, sqrt(dt));
      t += dt;
      x += v * dt;
      v += (gamma<T>(x) - alpha * v) * dt + beta * dW(*rng);
    bool IsInBounds() { return fabs(x) < half_width; };
    T t = t0;
    T dt = tau; 
    T x = x0;
    T v = v0;
    shared_ptr<rng_t> rng;

int main() {
  rng_t master_rng;
  uniform_int_distribution<unsigned> random_seed(314); 

  auto threads = omp_get_max_threads();
  vector<shared_ptr<rng_t>> rngs(threads);
  for (auto &r : rngs) r = make_shared<rng_t>(random_seed(master_rng));

  cout << "Using " << threads << " threads." << endl;

  cout << "Simulation running..." << endl;
  array<real_t, iters> data;
  #pragma omp parallel for
  for (int n = 0; n < iters; n++) {
    auto thread = omp_get_thread_num();
    OrnsteinUhlenbeck<real_t> particle(rngs[thread]);
    data[n] = particle.EscapeTime();
  cout << "Done." << endl;

  ofstream dataFile("data.dat");
  copy(data.begin(), data.end(), ostream_iterator<real_t>(dataFile, " "));

Afterward I plotted the result with the following Python script

import numpy as np
from scipy.stats import expon
import matplotlib.pyplot as plt

data = np.loadtxt("data.dat")

plt.title("Escape Time Statistics")
plt.xlabel("escape time")

loc, scale = expon.fit(data)
_, bins, _ = plt.hist(data, bins=30, normed=True)
x = np.linspace((bins[0] + bins[1])/2, (bins[-1]+bins[-2])/2, 100)
y = expon.pdf(x, loc, scale)

However, the plot of the statistics sometimes look like this

enter image description here

sometimes like this

enter image description here

depending on the value of the seed of the random_seed generator. How come?

  • 1
    $\begingroup$ Monte-Carlo should get different answers as a function of the random seed, though it should also converge to some similar distribution as you increase the number of Monte-Carlo samples. Are you sure the number of samples you have is sufficient to reduce the variance in the results enough that these histograms should be very close? $\endgroup$ – spektr Oct 20 '17 at 17:16
  • 5
    $\begingroup$ The y-axes are different in the two plots, and the histograms look very similar to my eye, it's the yellow lines that are different. Is it possible expon.fit is causing this instead of C++'s rng? $\endgroup$ – Kirill Oct 20 '17 at 18:03
  • $\begingroup$ @C.Howard From a qualitative point of view, the match between the exponential distribution and the data in the first plot seems pretty good. Moreover, changing the seed, the plot looks either as the first or the second, nothing in between. So I rule out that the problem is caused by a small data sample. @Kirill I thought about that the moment I submitted the question. However, whatever is the problem with expon.fit, it is triggered by changing the seed. $\endgroup$ – Stefano Campanella Oct 21 '17 at 7:06
  • $\begingroup$ However, the mean of the sample is almost the same in the two cases. The variance (scale paramenter in expon.fit parlance) change by an order of magnitude, smaller in the dataset within the one not matching the exponential distribution (as can be seen by the plot). $\endgroup$ – Stefano Campanella Oct 21 '17 at 7:12
  • $\begingroup$ I don't speak C++ very well, so I'm not sure I understand how you're seeding the Mersenne Twister PRNG. I would guess that either spektr or Anton M. has the correct explanation, but just in case it might be relevant: Mersenne Twisters can have a flaw in that if the initial many-word state contains mostly zero bits, the behavior can be nonrandom for tens of thousands of iterations. $\endgroup$ – Mars Jun 28 '19 at 6:44

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-number generation. The suspect is on the exponential fit from scipy that was unable to fit the function correctly and possibly had a bug in the implementation. That incorrect fit causes the orange line to skew the appearance of the plots.

I was very curious about this behavior, but was totally unable to reproduce it on my machine. So I would conclude with a high level of confidence that this bug in non-existent in my scipy+numpy+python toolchain.

My configuration (I did about 150 runs with different seeds visually checking the resultant PDFs):

  • GCC 6.4 (MacPorts gcc6 6.4.0_0)
  • compiled with -O3: /opt/local/bin/g++-mp-6 -openmp -O3 monte.cpp -lgomp
  • Python 2.7.10, native from Mac OS 10.13.3
  • numpy: 1.8.0rc1
  • scipy: 0.13.0b1
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.