モンテカルロ積分法を使用して、n球の体積を見つけ、それを分析結果と比較したいと思います。これを10 ^ 6ポイントと10 ^ 9ポイントで実行したいのですが、10 ^ 6ポイントでは多少機能しますが（n = 2（円）、n = 3（球）、n = 12の場合は約1分かかります）、 10 ^ 9ポイントはひどく遅いです。
MCメソッドの簡単な説明：半径r = 1のn球の体積を見つけるために、n球を完全に含む簡単に知られている体積（たとえば、辺の長さが2 * rのn立方体）を想像します。次に、n立方体の一様分布点からサンプリングし、その点が球内にあるかどうかを確認します。そのように生成されたn球内のすべての点を数えます。V_sphere / V_cubeの比率は、N_inside / N_totalとして概算できるため、V_sphere = V_cube * N_inside / N_total
関数は次のとおりです。
def hyp_sphere_mc(d,samples):
inside = 0 #number of points inside sphere
sum = 0 #sum of squared components
for j in range(0,samples):
x2 = np.random.uniform(0,1,d)**2 #generate random point in d-dimensions
sum = np.sum(x2) #sum its components
if(sum < 1.0):
inside += 1 #count points inside sphere
V = ((2)**d)*(float(inside)/samples) #V = V_cube * N_inside/N_total
V_true = float(math.pi**(float(d)/2))/math.gamma(float(d)/2 + 1) #analytical result
ERR = (float(abs(V_true-V))/V_true)*100 #relative Error
print "Numerical:", V, "\t" , "Exact: ", V_true, "\t", "Error: ", ERR
問題は、反復ごとに新しいランダム配列を生成することであり、特に10 ^ 9回の反復がある場合、これには多くの時間がかかります。これをスピードアップする方法はありますか？
ループを次のように置き換えることができます。
inside = np.sum(np.sum(np.random.rand(samples,d)**2,1)<1)
numpyを使用するときは、ループを回避するようにしてください。アイデアは、行列内のすべてのサンプルを一度に生成し、その後のすべての操作をベクトル化できるということです。
この記事はインターネットから収集されたものであり、転載の際にはソースを示してください。
侵害の場合は、連絡してください[email protected]
コメントを追加