inohilog

/var/log/inohiro.log

モンテカルロ法

ついに課題を好きな言語で書いてよい科目がでてきました。わーい。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();
			}
		}
	}
}