#include #include #include using namespace std; int main() { double sigma; int niter; double x; double delta; int seed; int bins[10000]; cout << "Input sigma: " << endl; cin >> sigma; cout << "Input delta: " << endl; cin >> delta; cout << "Number of iterations: " << endl; cin >> niter; cout << "Random seed: " << endl; cin >> seed; srand(seed); // Initialize histogram double sigma2=sigma*sigma*2; int nbins = 100; for(int i = 0; i < nbins; i++) bins[i] = 0; // Start loops int nwalkers = 100000; for(int n=1; n <= nwalkers; n++){ x = double(rand())/double(RAND_MAX)*2*sigma-sigma; for(int i = 1; i <= niter; i++){ double xt = double(rand())/double(RAND_MAX)*2*delta-delta; xt += x; double w = exp(-xt*xt/sigma2)/exp(-x*x/sigma2); // double w = exp(-(xt*xt-x*x)/sigma2); if(w >= 1){ // accept the move x = xt; } else { double r = double(rand())/double(RAND_MAX); if(r < w){ x = xt; // we accept the move } } } // Add to the histogram double xhist = (x+5*sigma)/(10*sigma); int bin = xhist*nbins; bins[bin]++; } int ntotal=0; for(int i=0; i < nbins; i++) ntotal += bins[i]; cout << ntotal << " " << nwalkers << endl; double norm = 0.; for(int i=0; i < nbins; i++) { norm += bins[i]/double(nwalkers); cout << 10*sigma/nbins*i-5*sigma << " " << bins[i]/double(nwalkers) << endl; } cout << norm << endl; }