· 

モンテカルロ法によるガウス積分の推定

(著)山たー

ソースコード

//estimation of Gaussian integral using Monte Carlo Method

float count=0;
float Num=0;
int N=1000;
float X[];
float Y[];

void setup(){
  size(700,650);
  background(255);
  colorMode(RGB);
  //ガウス関数の描画
  X= new float[N];
  Y= new float[N];
  for(int i=0;i<N;i++){
    X[i]=-3+(6.0*i/1000.0);
    Y[i]=gaussian(X[i]);
    //println(X[i], Y[i]);
    noStroke();
    fill(255,0,0);
    ellipse(350+X[i]*100,550-Y[i]*200,4,4);
  }
  //グラフの描画
  strokeWeight(3);
  noFill();
  stroke(0,0,255,200);
  rect(50,349,600,201); //6*1の長方形で考える
  fill(0);
  stroke(0);
  triangle(350,300,345,309,355,309); 
  triangle(675,550,667,555,667,545); 
  //textの描画
  textSize(20);
  fill(0);
  textAlign(CENTER);
  text("O",340,570);
  text("-3",50,570);
  text("3",650,570);
  text("1",340,345);
  text("x",670,570);
  text("y",340,300);
  textSize(30);
  text("Estimation of Gaussian integral",350,50);
  text("Using Monte Carlo Method",350,90);
  text("Number of Points=",160,150);
  text("Estimated Value=",170,200);
  text("Calculated Value=√π=1.772453…",287,250);
  text("y=exp(-x^2)",470,330);
}

void draw(){
  float x=random(-3,3);
  float y=random(0,1);
  Num+=1;
  if(y<=gaussian(x)){
    count +=1;
    fill(255,0,0);
  }else if(y>gaussian(x)){
    fill(0,0,255);
  }
  noStroke();
  ellipse(350+100*x,550-200*y,5,5);
  textSize(40);
  noStroke();
  fill(255);
  rect(300,160,300,60);
  rect(300,120,300,50);
  fill(0);
  stroke(0);
  textAlign(RIGHT);
  text(int(Num),470,155);
  text(nf(6*count/Num,1,5),470,205);
  stroke(0);
  line(30,550,670,550);
  line(350,600,350,300);
}

float gaussian(float x){
  return exp(-x*x);
}