
#0.0以上以上1.0以下の一様分布から乱数を生成するプログラム
import random

print(random.uniform(0.0,1.0))

0.5439912915686401


#一次関数関数f(x)=xを 0<=x<=1 で積分するプログラム
import random

# N：繰り返しの回数
N = 100
sum = 0
ave = 0

# 積分する関数
def function(x):
  return x

for i in range(N):
  s = random.uniform(0.0,1.0)
  sum += function(s)

#期待値の計算と出力
ave = sum / N
print(ave)

0.5247123710670295


%%writefile MC_sph.cpp

//=====================//
// モンテカルロ法による半径１の球の体積を求めるプログラム
// 乱数の範囲は 0 <= x,y <= 1 
//
//  S = {x>0, y>0, x^2+y^2<1}
//  P(x,y) = 4 / \pi (x^2+y^2<1)
//           0       (x^2+y^2>=1)
//  \int_S dxdy f(x,y) P(x,y) = 1 / (\pi/4) x \int_S dxdy f(x,y)
//=====================//

#include <stdlib.h>
#include <iostream>
#include <random>

int main() {
  //=====================//
  //  In c++, the <random> header is necessary.
  //  the random-distribution can be obtained by two steps.
  //  ENGINE      : genetates random values
  //  DISTRIBUTION: molds them as a target distribution
  //
  //  random_device: integer random value used as a seed
  //  mt19937      : pseudo-random value
  //=====================//

  int loop_N = 1000; // 何回体積を計算するか
  int niter = 10000; // １回の計算につき何回乱数を生成するか（期待値の和Nに対応）
  int rand_max = 10000000;

  double ave = 0.;
  double sav = 0.;
  double err = 0.;

  for (int k = 0; k < loop_N; k++) {
    std::random_device rd;
    std::mt19937 mt(rd());
    std::uniform_int_distribution<int> rand_dist(0, rand_max);

    int n_in = 0;
    double sum_z = 0.0;

    for (int iter = 0; iter < niter; iter++) {
      double x = (double)rand_dist(mt) / rand_max;
      double y = (double)rand_dist(mt) / rand_max;

      if (x * x + y * y < 1.) {
        n_in += 1;
        sum_z += sqrt(1. - x * x - y * y);
      }
    }
    ave += sum_z / n_in / loop_N * 2.0 * M_PI;
    sav += sum_z * sum_z / n_in / n_in / loop_N * 2.0 * M_PI * 2.0 * M_PI;
  }
  err = sqrt(sav - ave * ave);
  
  std::cout << "公式の結果：" << 4. * M_PI / 3. << std::endl;
  std::cout << "モンテカルロ法による評価：（中心値、誤差）＝（" << ave << "、 " << err << "）" << std::endl;

  return 0;
}

Overwriting MC_sph.cpp


!g++ MC_sph.cpp -o RUN.x; ./RUN.x

公式の結果：4.18879
モンテカルロ法による評価：（中心値、誤差）＝（4.18867、 0.00524855）


%%writefile Metropolis_gauss.cpp

//=====================//
//  Metropolis sampling
//
//  target: P(x) = exp(-0.5*x^2) / Z
//  action = 0.5*x^2
//=====================//

#include <stdlib.h>
#include <iostream>
#include <map>
#include <random>
#include <vector>

int main() {
  //=====================//
  //  In c++, the <random> header is necessary.
  //  the random-distribution can be obtained by two steps.
  //  ENGINE      : genetates random values
  //  DISTRIBUTION: molds them as a target distribution
  //
  //  random_device: integer random value used as a seed
  //  mt19937      : pseudo-random value
  //=====================//
  int niter = 100000;
  int rand_max = 10000000;
  double step_size = 0.5;

  //=====================//
  // Initial input
  // x=100 is far from distribution
  // -> Need termalization
  //=====================//
  double x = 0.;
  int naccept = 0;

  std::vector<double> list(niter);

  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_int_distribution<int> rand_dist(0, rand_max);

  for (int iter = 1; iter < niter + 1; iter++) {
    double backup_x = x;
    double action_init = 0.5 * x * x;

    //=====================//
    // Uniform distribution
    // [1] [-0.5:0.5] : satisfy the detailed balance (c=0.5)
    // [2] [-0.5:1.0] : break
    //=====================//
    //[1]
    double dx = (double)rand_dist(mt) / rand_max;
    dx = (dx - 0.5) * step_size * 2.;
    x += dx;

    //[2]
    // double dx = (double)rand_dist(mt) / rand_max * 0.75;
    // dx = (dx - 0.25) * 2.;
    // x += dx;

    double action_fin = 0.5 * x * x;

    //=====================//
    // Metropolis test
    //=====================//
    double metropolis = (double)rand_dist(mt) / rand_max;
    if (std::exp(action_init - action_fin) > metropolis) {
      naccept += 1;
    } else {
      x = backup_x;
    }
    list[iter - 1] = x;
  }

  //=====================//
  // Histogram
  //=====================//
  double range = 5.25;
  double step = 0.1;

  std::map<double, int> histogram;

  for (double i = -range; i <= range; i += step) {
    int number = 0;
    for (auto& pl : list) {
      if (i - step * 0.5 < pl && pl <= i + step * 0.5) {
        ++number;
      }
    }
    histogram.insert(std::make_pair(i, number));
  }

  // for (auto itr = histogram.begin(); itr != histogram.end(); ++itr) {
  //  std::cout << itr->first << " " << itr->second << std::endl;
  //}

  //=====================//
  // Thermalization
  //=====================//
  // for (int iter = 1; iter < niter + 1; iter++) {
  //   std::cout << iter << " " << list[iter - 1] << std::endl;
  // } 

  return 0;
}

Writing Metropolis_gauss.cpp


!g++ Metropolis_gauss.cpp -o RUN.x; ./RUN.x


!g++ --version

g++ (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


%%writefile main.cpp
//=====================//
//  Hybrid Monte-Carlo -Multi variables
//
//  target: P(\phi) = exp(-J[\phi]) / Z where \phi is NxN hermitian matrix
//  action = N Tr { 0.5 x \phi^2 + 0.25 x \phi^4 }
//
//  Leapfrog method
//  [1] intro n = 0
//  phi_ij( 0.5 * dt ) = phi_ij( 0 ) + pi_ij( 0 ) * 0.5 * dt
//
//  [2] main  n = 1 ~ Nt-1
//  pi_ij( n * dt )
//  = pi_ij( (n-1) * dt ) - \pdif H / \pdif phi_ij( (n-0.5) * dt ) * dt
//
//  phi_ij( (n+0.5) * dt )
//  = phi_ij( (n-0.5) * dt ) + pi_ij( n * dt ) * dt
//
//  [3] outro n= Nt
//  pi_ij( Nt * dt )
//  = pi_ij( (Nt-1) * dt ) - \pdif H / \pdif phi_ij( (Nt-0.5) * dt ) * dt
//
//  phi_ij( Nt * dt )
//  = phi_ij( (Nt-0.5) * dt ) + pi_ij( Nt * dt ) * dt * 0.5
//
//=====================//
#include <stdlib.h>
#include <algorithm>
#include <cmath>
#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <vector>

//=====================//
//  rand_max: uniform_int_distribution # 変更しないでください
//  niter: sample number
//  ntau : step number of leapfrog     # 変更しないでください
//  dtau : step size of leapfrog       # 変更しないでください       
//  ndim : number of variables
//  ninit: 1 -> new config., 0 -> old config.
//=====================//
const int rand_max = 1000000;
const int niter = 50000;
const int ntau = 20;
const double dtau = 0.05;
const int ndim = 2;
const int ninit = 1;

int BoxMuller(double& p, double& q, int r_max) {
  //=====================//
  //  In c++, the <random> header is necessary.
  //  the random-distribution can be obtained by two steps.
  //  ENGINE      : genetates random values
  //  DISTRIBUTION: molds them as a target distribution
  //
  //  random_device: integer random value used as a seed
  //  mt19937      : pseudo-random value
  //=====================//
  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_int_distribution<int> rand_dist(0, r_max);

  double pi = 2.0 * asin(1.0);

  double x = (double)rand_dist(mt) / r_max;
  double y = (double)rand_dist(mt) / r_max;

  p = sqrt(-2. * log(x)) * sin(2. * pi * y);
  q = sqrt(-2. * log(x)) * cos(2. * pi * y);

  return 0;
}

// 作用関数の計算
double calc_action(const std::complex<double> phi[ndim][ndim]) {
  double action = 0.;

  // \phi^2 パート
  std::complex<double> phi2[ndim][ndim];
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      phi2[idim][jdim] = std::complex<double>(0., 0.);
    }
  }
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      for (int kdim = 0; kdim != ndim; kdim++) {
        phi2[idim][jdim] += phi[idim][kdim] * phi[kdim][jdim];
      }
    }
  }
  for (int idim = 0; idim != ndim; idim++) {
    action += 0.5 * phi2[idim][idim].real();
  }

  // \phi^4 パート
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      action += 0.25 * (phi2[idim][jdim] * phi2[jdim][idim]).real();
    }
  }

  // 全体に掛ける因子
  action = action * ndim;

  return action;
}

// ハイブリッドモンテカルロ法に必要な計算（ハミルトニアンの微分）
int calc_delH(const std::complex<double> phi[ndim][ndim],
              std::complex<double> (&delH)[ndim][ndim]) {
  // \phi^2 と \phi^3 の初期化
  std::complex<double> phi2[ndim][ndim];
  std::complex<double> phi3[ndim][ndim];
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      phi2[idim][jdim] = std::complex<double>(0., 0.);
      phi3[idim][jdim] = std::complex<double>(0., 0.);
    }
  }
  // \phi^2 パート
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      for (int kdim = 0; kdim != ndim; kdim++) {
        phi2[idim][jdim] += phi[idim][kdim] * phi[kdim][jdim];
      }
    }
  }
  // \phi^3 パート
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      for (int kdim = 0; kdim != ndim; kdim++) {
        phi3[idim][jdim] += phi2[idim][kdim] * phi[kdim][jdim];
      }
    }
  }

  // ハミルトニアンの微分（delH）の初期化
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      delH[idim][jdim] = std::complex<double>(0., 0.);
    }
  }
  // delH パート
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      delH[idim][jdim] = (double)ndim * (phi[idim][jdim] + phi3[idim][jdim]);
    }
  }

  return 0;
}

// ハイブリッドモンテカルロ法のための計算（ハミルトニアンの計算）
double calc_H(const std::complex<double> phi[ndim][ndim],
              const std::complex<double> pi[ndim][ndim]) {
  // ポテンシャルエネルギー
  double H = calc_action(phi);

  // 運動エネルギー
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      H += 0.5 * (pi[idim][jdim] * pi[jdim][idim]).real();
    }
  }

  return H;
}

// ハイブリッドモンテカルロ法のための計算（分子動力学）
int Molecular_Dynamics(std::complex<double> (&phi)[ndim][ndim], double& H_ini,
                       double& H_fin) {
  // initialize pi
  std::complex<double> pi[ndim][ndim];
  for (int idim = 0; idim != (ndim - 1); idim++) {
    for (int jdim = idim; jdim != ndim; jdim++) {
      double r1, r2;
      BoxMuller(r1, r2, rand_max);
      pi[idim][jdim] =
          r1 / sqrt(2.) + r2 / sqrt(2.) * std::complex<double>(0., 1.);
      pi[jdim][idim] =
          r1 / sqrt(2.) - r2 / sqrt(2.) * std::complex<double>(0., 1.);
    }
  }
  for (int idim = 0; idim != ndim; idim++) {
    double r1, r2;
    BoxMuller(r1, r2, rand_max);
    pi[idim][idim] = r1;
  }

  // initialize delH
  std::complex<double> delH[ndim][ndim];
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      delH[idim][jdim] = std::complex<double>(0., 0.);
    }
  }

  H_ini = calc_H(phi, pi);

  // First step
  // CAUTION: 0.5
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      phi[idim][jdim] = phi[idim][jdim] + pi[idim][jdim] * 0.5 * dtau;
    }
  }
  // 2 ~ Ntau step
  for (int step = 1; step != ntau; step++) {
    calc_delH(phi, delH);

    for (int idim = 0; idim != ndim; idim++) {
      for (int jdim = 0; jdim != ndim; jdim++) {
        pi[idim][jdim] = pi[idim][jdim] - delH[idim][jdim] * dtau;
        phi[idim][jdim] = phi[idim][jdim] + pi[idim][jdim] * dtau;
      }
    }
  }
  // Final step
  // CAUTION: 0.5
  calc_delH(phi, delH);
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      pi[idim][jdim] = pi[idim][jdim] - delH[idim][jdim] * dtau;
      phi[idim][jdim] = phi[idim][jdim] + pi[idim][jdim] * 0.5 * dtau;
    }
  }

  H_fin = calc_H(phi, pi);

  return 0;
}

// ジャックナイフ用のテンプレートクラス
template <typename T>
class Jackknife {
 public:
  Jackknife(int n, T m) : N((double)n), M(m){};
  T operator()(const T& t) const {
    if (N == 1) return t;
    return (M * N - t) * (1. / (N - 1.));
  };

 private:
  double N;
  T M;
};

int main() {
  //=====================//
  //  In c++, the <random> header is necessary.
  //  the random-distribution can be obtained by two steps.
  //  ENGINE      : genetates random values
  //  DISTRIBUTION: molds them as a target distribution
  //
  //  random_device: integer random value used as a seed
  //  mt19937      : pseudo-random value
  //=====================//
  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_int_distribution<int> rand_dist(0, rand_max);

  //=====================//
  // 配位（configuration）の初期化
  // MCMC法では、生成した乱数の組みを「配位」と呼びます。
  // この配位生成は通常非常に時間がかかるため、時短のために既に生成した乱数を初期値とする場合が
  // 多いです。
  // このプログラムでは、ninit が１の時は全て０で初期化、０の時は読み込んで初期化します。
  // デフォルトは１です。  
  //=====================//

  std::complex<double> phi[ndim][ndim];
  if (ninit == 1) {
    for (int idim = 0; idim != ndim; idim++) {
      for (int jdim = 0; jdim != ndim; jdim++) {
        phi[idim][jdim] = std::complex<double>(0., 0.);
      }
    }
  } else {
    std::ifstream configufile("configuration.dat");

    for (int idim = 0; idim != ndim; idim++) {
      for (int jdim = 0; jdim != ndim; jdim++) {
        configufile >> phi[idim][jdim];
      }
    }
    configufile.close();
  }

  int naccept = 0;
  int nconf = 0;
  double sum_action = 0.;

  // confirmation of Thermalization
  std::vector<double> list(niter);

  // allocate the measured action
  std::vector<double> action(niter);

  for (int iter = 0; iter != niter; iter++) {
    std::complex<double> backup_phi[ndim][ndim];
    for (int idim = 0; idim != ndim; idim++) {
      for (int jdim = 0; jdim != ndim; jdim++) {
        backup_phi[idim][jdim] = phi[idim][jdim];
      }
    }

    double H_ini, H_fin;
    Molecular_Dynamics(phi, H_ini, H_fin);
    double metropolis = (double)rand_dist(mt) / rand_max;

    if (std::exp(H_ini - H_fin) > metropolis) {
      naccept += 1;
    } else {
      for (int idim = 0; idim != ndim; idim++) {
        for (int jdim = 0; jdim != ndim; jdim++) {
          phi[idim][jdim] = backup_phi[idim][jdim];
        }
      }
    }
    list[iter] = calc_action(phi);

    if (iter >= 20000 && iter % 100 == 0) {
      nconf += 1;
      sum_action += calc_action(phi);
      action.push_back(calc_action(phi));
    }
  }

  //=====================//
  // Thermalization（熱化）
  // 生成した乱数の組みが初期値によらないこと（熱化）を確認する必要があります。
  // このプログラムでは、乱数の値をテキストファイルに書き出しています。
  // 初期値の近くで値が大きくぶれていれば、その値は初期値に強く影響されており、大きな誤差のもと
  // になります。実際の計算では、その様な箇所は使いません。
  //=====================//
  std::ofstream action_monitor("action_monitor.dat");
  for (int iter = 0; iter != niter; iter++) {
    action_monitor << iter << " " << std::setprecision(16) << list[iter]
                   << std::endl;
  }
  action_monitor.close();

  //=====================//
  // Save final configuration
  //=====================//
  std::ofstream configfile("configuration.dat");
  for (int idim = 0; idim != ndim; idim++) {
    for (int jdim = 0; jdim != ndim; jdim++) {
      configfile << phi[idim][jdim] << std::endl;
    }
  }
  configfile.close();

  //=====================//
  // 平均値と誤差の計算
  // MCMC法では、自分自身の情報を使って乱数の生成を行っているので、データ間に相関を持ってしまっ
  // ています。そのため、相関を加味した誤差評価（ここではジャックナイフ法）を行うことになります。 
  //=====================//
  double ave_action = sum_action / action.size();
  transform(action.begin(), action.end(), action.begin(),
            Jackknife<double>(action.size(), ave_action));

  double ave = 0.;
  double sav = 0.;
  for (auto paction : action) ave += paction / action.size();
  for (auto paction : action)
    sav += (paction - ave) * (paction - ave) * (1. - 1. / action.size());
  std::cout << ave << " " << sqrt(sav) << std::endl;

  return 0;
}

Writing main.cpp


!g++ main.cpp -o MatrixModel.x; ./MatrixModel.x

0.00865388 0.000632883

スマホで始める素粒子原子核の物理¶

もっと知りたい人のために¶

数値シミュレーションについて¶

素粒子原子核の物理を記述する「量子場の理論」とは？¶

「スケール」と「物質の階層構造」¶

ドラマチックなミクロの世界¶

登場人物ー素粒子ー¶

脚本ー物理法則ー¶

実演ー計算ー¶

数値計算技術「マルコフ連鎖モンテカルロ法」とは？¶

コンピュータで積分の計算をする？？¶

積分の計算¶

確率と期待値¶

期待値と積分¶

モンテカルロ法と数値積分¶

乱数と確率分布¶

乱数と数値積分¶

マルコフ連鎖モンテカルロ法¶

メトロポリス法¶

MCMC法のまとめと応用例¶

実践！スマホで素粒子原子核物理¶

あとがき¶

参考文献¶