0

所以这就是我已经拥有的:

import numpy as np
import matplotlib.pyplot as plt

def monteCarloPi(n):
    np.random.seed()                 #seed the random number generator
    y = np.random.rand(n)*2 - 1      #n random samples on (-1,1)
    x = np.linspace(-1,1,n)          #x axis to plot against

    square = np.array([x,y])         #collecting axes as a single object

    mask1 = ((x**2 + y**2) < 1)      #filters

    hits = np.sum(mask1)             #calculating approximation
    ratio = hits/n
    pi_approx = ratio * 4

    return pi_approx

这是我想做的事情:

x = np.arange(100,1000)
y = monteCarloPi(x)
plt.scatter(x,y)

但是,当我运行上述代码块时,出现以下错误:

---------------------------------------------------------------------
TypeError                           Traceback (most recent call last)
<ipython-input-52-bf4dcedaa309> in <module>()
      1 x = np.arange(100,1000)
----> 2 y = monteCarloPi(x)
      3 plt.scatter(x,y)

    <ipython-input-51-8d5b36e22d4b> in monteCarloPi(n)
      1 def monteCarloPi(n):
      2     np.random.seed()                 #seed the random number generator
----> 3     y = np.random.rand(n)*2 - 1      #n random samples on (-1,1)
      4     x = np.linspace(-1,1,n)          #x axis to plot against
      5 

mtrand.pyx in mtrand.RandomState.rand()

mtrand.pyx in mtrand.RandomState.random_sample()

mtrand.pyx in mtrand.cont0_array()

TypeError: only integer scalar arrays can be converted to a scalar index

根据我对广播如何工作的理解numpy,这应该有效。我可以只使用一个for循环,但随着样本数量的增加,它会变得非常慢。

停止

4

1 回答 1

0

这是一种选择,其中基于最大样本量,然后在start>0(不包括错误处理)发生子采样。

import numpy as np
import matplotlib.pyplot as plt

def monteCarloPi(n,start=0,stride=1):
    np.random.seed()                 # seed the random number generator
    y = np.random.rand(n)*2 - 1      # n random samples on (-1,1)
    x = np.linspace(-1,1,n)          # x axis to plot against
    mask = ( x**2 + y**2 ) < 1       # masking
    samples = {}
    inds = arange(n)
    for k in range(n-start,n+1,stride):
        sub_inds = np.random.choice(inds,k,replace=False)
        sub_mask = mask[sub_inds]
        sub_hits = np.sum(sub_mask)
        ratio    = sub_hits/n
        pi_approx = ratio * 4
        samples[k]=pi_approx
    return pi_approx

这仍然需要一个 for 循环,但它在方法内部被快速处理,因为您是从一个大的随机样本中进行二次抽样。要重新创建您的原始呼叫(从n=100n=1000[注意,我要去n=1000这里]):

estimates = monteCarloPi(1000,start=900)
plt.plot(estimates.keys(),estimates.values())

您当然可以传递原始的x=arange(100,1001),但随后需要在方法中进行错误检查(以确保传递了数组或列表),然后n等于x( n=x[-1]) 的最后一个元素,最后是循环将在x( for k in x:) 的元素上完成。

于 2019-02-19T18:51:23.483 回答