冠层作图与随机化误差

Question

import math
import pylab as plt
import numpy
from numpy import sqrt
from scipy.integrate import quad
import random

numpy.seterr(divide='ignore', invalid='ignore')

#this integrates sqrt(1-x^2) from 0 to 1 and stores the value in q1area
def integrand (x):
return sqrt(1-x**2)
q1area, err = quad(integrand,0,1)

print "This program estimates the convergence of Pi to a ratio of one."

#code enters a while loop that runs until the user decides to exit the code (choice 4)
while True:
    print "Please choose from one of the five following options:"
    print " 1. 10^1\n 2. 10^2\n 3. 10^3\n 4. Exit"
    choice = int(raw_input())
    options = {1,2,3,4}

if choice == 1:
    #sets the x and y bounds on the graph
    plt.xlim([0,12])
    plt.ylim([-5,5])

    #sets the graph for the program to store random points on - sqrt(1-x^2)
    x = numpy.linspace(0,12,500)
    y = numpy.sqrt(1-x**2)
    #the program says that x and y need to be of the same dimension
    #if z doesn't have the x*0 term. That term seems to fix the error
    z = 1+x*0

    xcord = []
    ycord = []
    under = []
    above = []
    pratiolist = []
    yvalues = []
    xvalues = range(1,11)

    #appends a random point on the graph
    for i in range(10):
        xcord.append(random.random())
        ycord.append(random.random())

    #checks to see if the point is under or above the curve then assigns it to 
    #the corresponding list
    for j in ycord:
        if (j <= q1area):
            under.append(1)

        else:
            above.append(1)

        #defines two variables value as the length of a list
        punder = len(under)
        pabove = len(above)

        #checks for division by zero, adds one if true
        if pabove == 0:
            pabove = pabove + 1

        #calculates pratio as points above and below, creates list, appends values
        pratio = punder / float(pabove)
        pratiolist.append(pratio)

    #runs through pratiolist and calculates ratio to pi for each point
    for k in pratiolist:
        rtpi = k / float(math.pi)
        yvalues.append(rtpi)

    #displays graph
    plt.scatter(xvalues,yvalues,c='b')
    plt.plot(x,z,'g')
    plt.show()

if choice == 2:
    plt.xlim([0,120])
    plt.ylim([-5,5])

    x = numpy.linspace(0,120,500)
    y = numpy.sqrt(1-x**2)
    z = 1+x*0

    xcord = []
    ycord = []
    under = []
    above = []
    pratiolist = []
    yvalues = []
    xvalues = range(1,101)

    for i in range(100):
        xcord.append(random.random())
        ycord.append(random.random())

    for j in ycord:
        if (j <= q1area):
            under.append(1)

        else:
            above.append(1)

        punder = len(under)
        pabove = len(above)

        if pabove == 0:
            pabove = pabove + 1

        pratio = punder / float(pabove)
        pratiolist.append(pratio)

    for k in pratiolist:
        rtpi = k / float(math.pi)
        yvalues.append(rtpi)

    plt.scatter(xvalues,yvalues,c='b')
    plt.plot(x,z,'g')
    plt.show()

if choice == 3:
    plt.xlim([0,1100])
    plt.ylim([-5,5])

    x = numpy.linspace(0,1100,500)
    y = numpy.sqrt(1-x**2)
    z = 1+x*0

    xcord = []
    ycord = []
    under = []
    above = []
    pratiolist = []
    yvalues = []
    xvalues = range(1,1001)

    for i in range(1000):
        xcord.append(random.random())
        ycord.append(random.random())

    for j in ycord:
        if (j <= q1area):
            under.append(1)

        else:
            above.append(1)

        punder = len(under)
        pabove = len(above)

        if pabove == 0:
            pabove = pabove + 1

        pratio = punder / float(pabove)
        pratiolist.append(pratio)

    for k in pratiolist:
        rtpi = k / float(math.pi)
        yvalues.append(rtpi)

    plt.scatter(xvalues,yvalues,c='b')
    plt.plot(x,z,'g')
    plt.show()

if choice == 4:
    break

while choice not in options:
    print "Not a valid choice!\n"

一些空白是关闭的，我真的不知道如何制表符多行。无论如何，这个代码取10、100或1000个点，并在y=0和1，x=0和1之间的sqrt（1-x^2）图上随机断言它们。你知道吗

然后取这些点的比率，并将其与Pi的实际值进行比较，得到一个收敛图，该图将等于1（要估计的点越多越好）。你知道吗

我有两个问题。第一个可能很明显…第二张图上的那些点并不像我想的那样随机。它们似乎以线性方式上升并重新开始。你知道吗

第二个问题是，我希望这个代码循环，让用户反复选择三个选项中的一个，直到他们选择退出。每一次，这应该清除图表，使一个新的随机一切再次。这不是我的行为。有人能给我指出正确的方向吗？你知道吗

Answer 1

对于问题2：这样做的原因是计算曲线上的点与曲线下的点的比率的方法。当pabove随着时间的推移而增加时，比率会逐步变化。随着步数的增加，这应该会变得越来越不明显。
对于问题1：请参阅如何处理无休止的输入循环的代码

一些备注（后面代码中的许多其他提示）：
-您似乎忽略了绘制随机点的命令（图1）。我在代码中留下了数组定义，否则这些定义将是未使用/多余的。
-将y值与积分q1area进行比较。如果只看尺寸（长度到面积），这是不对的。我已经将比较替换为它真正应该是什么-将点位置（使用x和y）与曲线上的点进行比较。曲线q1area下面积的度量实际上根本没有使用。你知道吗

以下是我无法运行或测试的代码建议，因为我没有安装所有模块：

import math
import pylab as plt
import numpy
from numpy import sqrt
from scipy.integrate import quad
from random import random

### numpy.seterr(divide='ignore', invalid='ignore')

# function to integrate
def integrand (x):
    return sqrt(1-x**2)

# calculate integral with error estimate
q1area, error_estimate = quad(integrand, 0, 1)

print "This program estimates the convergence of Pi to a ratio of one."
pi = math.pi

while True:
    choice = 0
    options = [1, 2, 3, 4]
    while choice not in options:
        print "Please choose from one of the four following options:"
        print " 1. 10^1 steps\n 2. 10^2 steps\n 3. 10^3 steps\n 4. Exit"
        choice = int(raw_input())
    if choice == 4:
        break

    n = 10**choice # number of random samples

    # set the x and y bounds on the graph
    plt.xlim([0, 1.2*n])
    plt.ylim([-5, 5])

    # set the graph for the program to store random points on - sqrt(1-x^2)
    x = numpy.linspace(0, 1.2*n, 500)
    y = numpy.sqrt(1 - x**2)
    # x and z need to be of the same dimension
    # if z doesn't have the x*0 term. That term seems to fix the error
    z = 1 + x*0 # constant 1 over the range of x values

    xcoord = []
    ycoord = []
    pratiolist = []
    xvalues = range(1, n+1)
    yvalues = []

    # append n random points on the graph
    for i in range(n):
        xcoord.append(random())
        ycoord.append(random())

        # count if the point at (xcoord[i],ycoord[i]) is under or above the curve
        if ycoord[i] <= integrand(xcoord[i])::
            punder += 1
        else:
            pabove += 1

        # calculates pratio as points above and below, creates list, appends values
        if pabove > 0:
           pratiolist.append(punder / float(pabove))
        else:
           pratiolist.append(punder)

    # run through pratiolist and calculates ratio to pi for each point
    for k in pratiolist:
        rtpi = k / float(pi)
        yvalues.append(rtpi)

    # display graph
    plt.scatter(xvalues, yvalues, c='b')
    plt.plot(x, z, 'g')
    plt.show()

您将看到，您不必像开始时那样使用那么多的数组，例如用于计数。
主要的比较是使用integrand()函数，这样当您更改要集成的函数时就不必更改代码。你知道吗