モンテカルロ法
ついに課題を好きな言語で書いてよい科目がでてきました。わーい。Windowsが死んでるので、MonoDevelopで書き始めたんですが、見えない文字を打ってめんどくさい事になってしまったので(どうにかならんかな...)、全部emacsで書いたよ。
Random.Next でLINQを使ってランダム値を取ってくるみたいなコードがC#デザインパターン - Observerにありましたが、奇特なコードを書くと読んでもらえない(読めない)かもしれないので、普通に読めそうな感じに書きました。先生はJavaは読める、という情報はあります。
(πは一応指定が。RANDOM_MAXはstdlib.hのRAND_MAXだっけ?と同じ値をセット)
その1
任意の点の数を受け取って、それについてモンテカルロ法を10回行い、その平均値を出力するコード。
using System; namespace MonteCarloMethod { class MainClass { private static double Calc( double points ) { const int RANDOM_MAX = 21478643; double r, s, n = 0; Random rand = new Random(); for( int i = 0; i < points; i++ ) { r = ( double )( 1.0 / ( RANDOM_MAX + 1.0 ) ) * rand.Next( RANDOM_MAX ); s = ( double )( 1.0 / ( RANDOM_MAX + 1.0 ) ) * rand.Next( RANDOM_MAX ); if( Math.Pow( r, 2 ) + Math.Pow( s, 2 ) < 1.0 ) n++; } return ( double ) 4 * n / points; } public static void Main(string[] args) { const double PI = 3.14159; Console.WriteLine( "input points number: " ); double points = ( double )int.Parse( Console.ReadLine() ); double sum = 0.0; for( int i = 0; i < 10; i++ ) { sum += Calc( points ); } double avg = ( double ) sum /10; Console.WriteLine( string.Format( "Avg. = {0}", avg ) ); Console.WriteLine( string.Format( "( PI - Avg. ) / PI = {0}", ( PI - avg ) / PI ).ToString() ); } } }
実行例
πに近づいてく。
133051082075:console inohiro$ mono Main1.exe input points number: 1000000 Avg. = 3.1419532 ( PI - Avg. ) / PI = -0.000115610248313684
その2
-1<=x<=1, -1<=y<=1, x*x+y*y<=1 な条件を満たすx, yを適当に作って、それをCSVに吐くコード。(Excelとかで散布図にして、半径1の単位円ぽくなってるよね、みたいな確認をする)。無駄にクラスとか定義してる。
using System; namespace MonteCarloMethod { public class Point { public Point( double x, double y) { this.x = x; this.y = y; } public double x { get; set; } public double y { get; set; } } class MainClass { private static Point[] Calc() { const int RANDOM_MAX = 21478643; double x, y; Point[] points = new Point[200]; Random rand = new Random(); int i = 0; while( true ) { x = ( double )( 1.0 / ( RANDOM_MAX + 1.0 ) * rand.Next( -100000000, 100000000 ) ); // 適当 y = ( double )( 1.0 / ( RANDOM_MAX + 1.0 ) * rand.Next( -100000000, 100000000 ) ); // 適当 if( ( -1.0 <= x ) && ( x <= 1.0 ) ) { if( ( -1.0 <= y ) && ( y <= 1.0 ) ) { if( Math.Pow( x, 2 ) + Math.Pow( y, 2 ) <= 1.0 ) { points[i] = new Point( x, y ); i++; if( i == 200 ) break; } } } } return points; } public static void Main( string[] args ) { Point[] points = Calc(); using( System.IO.StreamWriter stWriter = new System.IO.StreamWriter( @"data.csv" ) ) { for( int i = 0; i < 200; i++ ) stWriter.WriteLine( "{0}, {1}", ( points[i].x ).ToString(), ( points[i].y ).ToString() ); stWriter.Close(); } } } }