Tuesday, February 10, 2015

Estimating Pi number with Monte Carlo Simulation in Java and Kdb+ q

package com.denizstij.montecarlo.sim;
 
import java.util.Random;
 
/**
*  This class estimates Pi with monte carlo simulation.
* 
 *  We assume we have a square with dimension 1x1 which contains a disk in middle with radious of 0.5
*  where it touches the square on 4 edges. We generate N number of pairs (x,y) between 0 and 1.
 *  Then we check if each pair is in the circle or not. 4 times ratio of total pairs in circle
*  to total simulation is an estimation of Pi. 
 * 
 *  pi=4*(total pairs in circle)/total pairs
* 
 * @author Deniz Turan, http://denizstij.blogspot.co.uk, denizstij AT gmail.com
*
 */
public class PIEstimationWithMonteCarloSimulation {
       // number of total simulation
       private static final long TOTAL_SIMULATION=10000000L;
      
       private double estimatePI(){            
              double pi=0;
              Random randomGen= new Random();
              long inCircleCounter=0;
              for (int i=0;i<TOTAL_SIMULATION;i++){
                     // generate pairs
                     double x = randomGen.nextDouble();
                     double y = randomGen.nextDouble();
                     // check if the pair  in the circle
                     inCircleCounter+=isInCircle(x,y);
              }
              double ratio=(double)inCircleCounter/TOTAL_SIMULATION;
              pi=ratio*4;
              return pi;
       }
      
       /*
       *  Check given two points in a circle: x^2+y^2<=1
       */
       private int isInCircle(double x, double y) {
              double z=Math.sqrt(x*x+y*y);            
              return z<=1?1:0;
       }
 
       public static void main(String[] args) {
              PIEstimationWithMonteCarloSimulation est= new PIEstimationWithMonteCarloSimulation();
              double estimatedPI = est.estimatePI();
              System.out.println("Estimated PI:"+estimatedPI);
       }
}
And now in Kdb+ q :
piEst:{[N]
        sqrArea:4;                                  / area of square with length 2 and containing a circle with radius 1
        xy:flip {2?1.0} each til N;                 / xy points
        p:sqrt sum xy*xy;                           / sqrt(x^2+y^2) <=1 to be in a circle
        totalIncircle: sum p<=1;
        piEst:sqrArea*totalIncircle%N;                    / that’s the estimation of PI
        :piEst
}
q)piEst[1000000]
3.141476
Above code can be reduced to one line such as:
q)N:10000000
q)(1%N)*4*sum 1>=sqrt sum ({x*x}') each flip {2?1.0} each til N
3.141842

1 comment:

Sarah Jordan said...

Excellent post. You have shared some wonderful tips. I completely agree with you that it is important for any blogger to help their visitors. Once your visitors find value in your content, they will come back for more Monte Carlo Simulation Stock Trading Systems .