28

numpy を使用して、球の表面にランダムな点を生成しようとしています。ここで一様分布を説明する投稿を確認しました。ただし、球の表面にのみポイントを生成する方法についてのアイデアが必要です。これらの各球の座標 (x、y、z) と半径があります。

私はこのレベルの数学に精通しておらず、モンテカルロ シミュレーションを理解しようとしています。

どんな助けでも大歓迎です。

ありがとう、パリン

4

5 に答える 5

46

このページの最後のアプローチに基づいて、 3 つの標準正規分布から独立したサンプルで構成されるベクトルを単純に生成し、大きさが 1 になるようにベクトルを正規化できます。

import numpy as np

def sample_spherical(npoints, ndim=3):
    vec = np.random.randn(ndim, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

例えば：

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d

phi = np.linspace(0, np.pi, 20)
theta = np.linspace(0, 2 * np.pi, 40)
x = np.outer(np.sin(theta), np.cos(phi))
y = np.outer(np.sin(theta), np.sin(phi))
z = np.outer(np.cos(theta), np.ones_like(phi))

xi, yi, zi = sample_spherical(100)

fig, ax = plt.subplots(1, 1, subplot_kw={'projection':'3d', 'aspect':'equal'})
ax.plot_wireframe(x, y, z, color='k', rstride=1, cstride=1)
ax.scatter(xi, yi, zi, s=100, c='r', zorder=10)

ここに画像の説明を入力

ndim=2同じ方法を一般化して、単位円 ( ) または高次元単位超球面上の均一に分布した点を選択します。

于 2015-11-28T23:11:15.127 に答える
8

theta球面上の点は、とphi、 、の 2 つの球座標を使用して表すことができ0 < theta < 2piます0 < phi < pi。

x, y, zデカルト座標への変換式:

x = r * cos(theta) * sin(phi)
y = r * sin(theta) * sin(phi)
z = r * cos(phi)

ここでrは球の半径です。

したがって、プログラムはランダムにサンプリングthetaしphi、その範囲内で一様分布で、そこからデカルト座標を生成できます。

しかし、その後、点は球の極でより密に分散されます。ポイントが球面上に均一に分布するようにするには、 を均一な分布phiで選択する必要がありphi = acos(a)ます-1 < a < 1。

Numpy コードの場合、変数が固定値を持つことを除いて、球状ボリューム内の均一に分散されたランダム ポイントのサンプリングと同じです。radius

于 2015-11-28T22:18:22.173 に答える
3

ハードウェアに依存する別の方法は、はるかに高速になる可能性があります。

a, b, cそれぞれ-1から1の間の3つの乱数を選択してください

計算するr2 = a^2 + b^2 + c^2

r2 > 1.0 (= ポイントが球内にない) または r2 < 0.00001 (= ポイントが中心に近すぎるため、球の表面への投影中にゼロ除算が行われる) の場合、値を破棄します。 、別のランダムなセットを選択しますa, b, c

それ以外の場合は、(球の中心に相対的な) ランダムな点が得られます。

ir = R / sqrt(r2)
x = a * ir
y = b * ir
z = c * ir
于 2015-11-28T22:32:11.100 に答える
3

@Soonts との議論の後、回答で使用された 3 つのアプローチのパフォーマンスに興味を持ちました。

これが私の試みた比較です：

import numpy as np

def sample_trig(npoints):
    theta = 2*np.pi*np.random.rand(npoints)
    phi = np.arccos(2*np.random.rand(npoints)-1)
    x = np.cos(theta) * np.sin(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(phi)
    return np.array([x,y,z])

def sample_normals(npoints):
    vec = np.random.randn(3, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

def sample_reject(npoints):
    vec = np.zeros((3,npoints))
    abc = 2*np.random.rand(3,npoints)-1
    norms = np.linalg.norm(abc,axis=0) 
    mymask = norms<=1
    abc = abc[:,mymask]/norms[mymask]
    k = abc.shape[1]
    vec[:,0:k] = abc
    while k<npoints:
       abc = 2*np.random.rand(3)-1
       norm = np.linalg.norm(abc)
       if 1e-5 <= norm <= 1:  
           vec[:,k] = abc/norm
           k = k+1
    return vec

そしたら1000ポイント

In [449]: timeit sample_trig(1000)
1000 loops, best of 3: 236 µs per loop

In [450]: timeit sample_normals(1000)
10000 loops, best of 3: 172 µs per loop

In [451]: timeit sample_reject(1000)
100 loops, best of 3: 13.7 ms per loop

拒否ベースの実装では、最初にサンプルを生成npointsし、悪いサンプルを破棄し、残りのポイントを生成するためにループのみを使用したことに注意してください。直接の段階的な拒否は、より長い時間がかかる場合があるようです。ケースとの比較をより明確にするために、ゼロ除算のチェックも削除しましたsample_normals。

2 つのダイレクト メソッドからベクトル化を削除すると、それらは同じ球場に置かれます。

def sample_trig_loop(npoints):
    x = np.zeros(npoints)
    y = np.zeros(npoints)
    z = np.zeros(npoints)
    for k in range(npoints):
        theta = 2*np.pi*np.random.rand()
        phi = np.arccos(2*np.random.rand()-1)
        x[k] = np.cos(theta) * np.sin(phi)
        y[k] = np.sin(theta) * np.sin(phi)
        z[k] = np.cos(phi)
    return np.array([x,y,z])

def sample_normals_loop(npoints):
    vec = np.zeros((3,npoints))
    for k in range(npoints):
      tvec = np.random.randn(3)
      vec[:,k] = tvec/np.linalg.norm(tvec)
    return vec

In [464]: timeit sample_trig(1000)
1000 loops, best of 3: 236 µs per loop

In [465]: timeit sample_normals(1000)
10000 loops, best of 3: 173 µs per loop

In [466]: timeit sample_reject(1000)
100 loops, best of 3: 14 ms per loop

In [467]: timeit sample_trig_loop(1000)
100 loops, best of 3: 7.92 ms per loop

In [468]: timeit sample_normals_loop(1000)
100 loops, best of 3: 10.9 ms per loop
于 2015-11-29T13:47:50.610 に答える