统计计算 课后习题答案汇 田霞 第1-6章

习题1.11、计算球的体积，相对误差限为1%，则测量半径R时允许的相对误差限为多少？解：体积公式为v4x3，若自变量x有绝对误差，函数的相对误差限为3e(f(x))f(x)f(x)f(x)3xe(f(x)) 3e(x)1%，e(f(x))f(x)f(x)f(x)3xr r r对误差限为0.00333。2、某圆柱体高度h的近似值h*=20cm，半径r的近似值r*=5cm，已知hh*

0.2cm,rr*

0.1cm，求圆柱体体积的绝对误差限与相对误差限。vv*

r2hr*2h*

hh*

rr*

hh*

rr*vhvrvhvrr*2hh*2r*h*rvhvrvhvrvvv*

r*20.22r*h*0.125。相对误差限v*

4%。52203xx*=0.236，误差限为0.5105，则该近似值有几位有效数字？xx*

0.5105，x*0.236100，因此近似值有6位有效数字。4、圆周率=3.14159，它的具有4位有效数字的近似值*为多少？解：=3.14159=0.31459101，*0.51014，所以0.5103*0.5103，即3.14109*3.14209。5、已知x为正数，x的相对误差限为，求lnx的绝对误差限和相对误差限。xxx*x*解：绝对误差限：lnxlnx ，相对误差限：

。lnxlnx*lnx*xxlnxlnx*lnx*xx*x*lnx*lnx*xn(x*)nn(x*)n1(xx*)n(xx*)解：(x*)n

。(x*)n x*7、自变量x的四舍五入的近似值x*=0.1753，则该近似值的误差限为多少？解：xx*0.0000a0.5104，a0,1,2,3,4，或xx*0.0000(10a)0.5104，a5,6,7,8,9因此绝对误差限为0.5104。8、3.141，3.14，3.15分别作为的近似值，确定它们各有几位有效数字。解：3.141=0.3141×101，3.14=0.314×101，3.15=0.315×101，xx*3.14159263.1410.00059260.0050.5102，因此有3位有效数字。xx*3.14159263.140.00159260.0050.5102，因此有3位有效数字。xx*3.14159263.150.008407350.050.5101，因此有2位有效数字。n9、已知积分s11xnn0x5

(n0,1,)，分别使用s0和s10作为初值，计算积分sn的值。1解：对该积分使用分部积分法可得递推公式：11s 11

xndx 1

xn1(x55)dx

1(xn1 5

xn1)dx15s ，n0x5

0x5

0 x5

n n111s 1 dxln60.1823，1s 111

x10dx1

，取这两个数的中值00x5 5

66 10

0x5 551(11≈0.01667s10作的近似值。266 55（1）先从s0作为初值开始计算。算法为：0 s0.18230s15s ,k2,,10代码：importnumpyasnps0=0.1823n=10forkinrange(n+1):s1=1/(k+1)-5*s0

kk k1print(k+1,s0)s0=s1输出结果：10.182320.0885000000000000230.05749999999999988540.0458333333333338950.0208333333333305460.095833333333347317-0.3125000000000699481.70535714285749259-8.4017857142874631042.12003968254843（2）先从s10作为初值开始计算。算法为：s

s100.016671k,k1,2代码：t0=0.01667n=10forkinrange(n,0,-1):t1=1/(5*k)-t0/5print(k,t1)t0=t1输出结果：100.01666690.0188890222222222280.02122219555555555670.02432698946031745860.0284679354412698450.03430641291174603440.04313871741765079630.05803892318313650520.088392215363372710.18232155692732546

k1k 5习题1.21、求方程xxe0的根，精度要求为10-5。解：两边取对数，得：xlnx1，求出方程xlnx10的根即可。令f(x)xlnx1，

f(x)lnx1，

f(x)1x

，牛顿法迭代公式为：x x

xklnxk1。下面选择初值。因为f(1)1，f(2)0.38629，而f(x)1kk

lnxk1 x迭代的轮数i输出值误差11.771850.2281521.763240.0086131.763220.00002在区间[1,2]为正的，因此在区间[1,2]f迭代的轮数i输出值误差11.771850.2281521.763240.0086131.763220.00002步骤：x0=2n=20；第二步，定义函数和函数的导数；n，设置终止条件为nn0.00001。输出结果：

，其中0迭代值1.7718483274489238绝对误差为0.228151672551076161迭代值1.7632362113366402绝对误差为0.008612116112283632迭代值1.7632228343842757绝对误差为1.3376952364474448e-052、用牛顿迭代法求方程x33x10的根，精度要求为10-4。解：令f(x)x33x1，f(x)3x23，f(x)6x，牛顿法迭代公式为：x33x1xk

xk

k k f(2)1k3x23k

f(x)6x在区迭代的轮数i输出值误差11.888880.11111间[1,2]为正的，因此在区间[1,2]f(2)f(x)迭代的轮数i输出值误差11.888880.1111121.879450.0094331.879390.00006因此根为1.87939。步骤：x0=2n=50；第二步，定义函数和函数的导数；第三步，按照迭代公式计算n，设置终止条件为nn0.0001。输出结果：1迭代值1.8888888888888888绝对误差为0.11111111111111116

，2迭代值1.879451566951567绝对误差为0.009437321937321833迭代值1.879385244836671绝对误差为6.632211489598916e-053(x)xc(x2)（1c{x}收敛于3（2）c取何值时收敛速度最快。33（1）(x)1x(3)1c3（2）(3)0时收敛速度最快，因此c1 。23

1，因此

1c034yx1+x2-x3，x1是区间[-2,5]是区间[2,6]是区间[-5,2]y[3,5,2]。轮数解集函数值最大值最优解1{(3,5,2)2,5,2(4,5,2)，(3,4,2)，(3,6,2)，(3,5,1)}6，5，7，5，7，75(2,5,2)2{(2,5,2)（15,2(3,5,2)(2,4,2)，(2,6,2)，(2,5,1)}5，4，6，4，6，64（1,5,2）3{(1,5,2)0,5,2(2,5,2)，(1,4,2)，(1,6,2)，(1,5,1)}4，3，5，3，5，53(0,5,2)4{(0,5,2)（1,5,2，(1,5,2)，(0,4,2)，(0,6,2)，(0,5,1)}3，2，4，2，4，42(-1,5,2)5{(-1,5,2)-2,5,2，(05,2)，(-1,4,2)，(-1,6,2)，(-1,5,1)}2，1，3，1，3，31(-2,5,2)6{(-2,5,2)-1,5,2，(-,4,2)，1，2，0，2，20(-2,4,2)(-2,6,2)，(-2,5,1)}7{(-2,4,2)-1,4,2，(-,3,2)，(-2,5,2)，(-2,4,1)}0，1，-1，1，1-1(-2,3,2)8{(-2,3,2)-1,3,2，(-,2,2)，(-2,4,2)，(-2,3,1)}-1，0，-2，0，0-2(-2,2,2)}9{(-2,2,2)-1,2,2，(-,3,2)，(-2,2,1)}-2，-1，-1，-1-2(-2,2,2)输出结果：7[2,5,2]6[1,5,2]5[0,5,2]4[-1,5,2]3[-2,5,2]2[-2,4,2]1[-2,3,2]0[-2,2,2]-1[-2,2,2]最优解为：-2取得最优解时自变量为[-2,2,2]5yx1+x2-x3，x1是区间[-2,5]是区间[2,6]是区间[-5,2]y[3,5,2]。解：轮数解集函数值最大值最优解1{(3,5,2)2,5,2(4,5,2)，(3,4,2)，(36,2)（3,5,1）}6，5，7，5，7，77(2,5,2)2{(4,5,2)3,5,2(5,5,2)，(4,4,2)，(46,2)（4,5,1）}7，6，8，6，8，88(5,5,2)3{(5,5,2)4,5,2(5,4,2)，(5,6,2)（55,1）}8，7，7，9，99(5,6,2)4{(5,6,2)4,6,2，(5,5,2)，（5,6,1）}9，8，8，1010(5,6,1)5{(5,6,1)4,6,1(5,5,1)，（5,6,0，(5,6,2)}10，9，9，11，911(5,6,0)6{(5,6,0)4,6,0，(5,5,0)，（5,6,-1，(5,6,1)}1012(5,6,-1)7{(5,6,-1)（4,6,-1，(5,5,-1)，13(5,6,-2)（5,6,-2，(5,6,0)}118{(5,6,-2)（4,6,-2，(5,5,-2)，(5,6,-3)，(5,6,-1)}1214(5,6,-3)}9{(5,6,-3)，(4,6,-3)，(5,5,-3)，(5,6,-4)，(5,6,-2)}1315(5,6,-4)10{(5,6,-4)，(4,6,-4)，(5,5,-4)，(5,6,-5)，(5,6,-3)}1416(5,6,-5)11{(5,6,-5)，(4,6,-5)，(5,5,-5)，(5,6,-4)}16，15，15，1516(5,6,-5)输出结果：7[4,5,2]8[5,5,2]9[5,6,2]10[5,6,1]11[5,6,0]12[5,6,-1]13[5,6,-2]14[5,6,-3]15[5,6,-4]16[5,6,-5]16[5,6,-5]最优解为：16取得最优解时自变量为[5,6,-5]6yx1+x2-x3，x1是区间[-2,5]是区间[2,6]是区间[-5,2]y取值最小的解。初值在范围内随机选取。输出结果：0[-1,4,-5]-1[-2,4,-5]-2[-2,3,-5]-3[-2,2,-5]-4[-2,2,-5]最优解为：-5取得最优解时自变量为[-2,2,-5]分析：此时的初值为[-2,2,-3]，得到的最小值是-5。7、使用梯度下降法求函数yx22x5的最小值，初值为3，学习率为0.4。x=-14。输出结果：迭代次数0-0.200000000000000184.64迭代次数1-0.84000000000000014.0256迭代次数2-0.9684.001024迭代次数3-0.99364.00004096迭代次数4-0.998724.0000016384迭代次数5-0.9997444.000000065536迭代次数6-0.99994884.00000000262144迭代次数7-0.999989764.0000000001048575迭代次数8-0.9999979524.000000000004194迭代次数9-0.99999959044.0000000000001688zx2y2的最小值，初值为（1，3，学习率为0.4。解：很容易看出二元函数的最小值为0，此时(x,y）=(0,0)。输出结果为：迭代次数00.6000000000000000.2000000000000007.42638585493563e-9迭代次数10.1200000000000000.04000000000000007.42638585493563e-9迭代次数20.02400000000000000.007999999999999997.42638585493563e-9迭代次数30.004799999999999990.001600000000000007.42638585493563e-9迭代次数40.0009599999999999980.0003200000000000007.42638585493563e-9迭代次数50.0001920000000000006.39999999999999e-57.42638585493563e-9迭代次数63.83999999999999e-51.28000000000000e-57.42638585493563e-9结果分析：可以看出极小值为0，此时(x,y）=(0,0)9、使用梯度下降法求二元函数z1x21y2的最小值选取初值为(2,2)，学习率3 2α=6/5。解：很容易看出二元函数的最小值为0，此时(x,y）=(0,0)。输出结果：迭代次数00.400000000000000-0.4000000000000000.133333333333333迭代次数10.08000000000000010.08000000000000000.00533333333333333迭代次数20.0160000000000000-0.01600000000000000.000213333333333333迭代次数30.003200000000000000.003200000000000008.53333333333334e-6迭代次数40.000640000000000001-0.0006400000000000003.41333333333334e-7迭代次数50.0001280000000000000.0001280000000000001.36533333333333e-8迭代次数62.56000000000001e-5-2.56000000000000e-55.46133333333334e-10结果分析：可以看出极小值为0，此时(x,y）=(0,0)10z2(x1)2y2α=0.4。解：当(x,y）=(1,0)。二元函数取得最小值，为0。输出结果：迭代次数00.4000000000000000.4000000000000000.880000000000000迭代次数11.360000000000000.08000000000000000.265600000000000迭代次数20.7840000000000000.01600000000000000.0935680000000001迭代次数31.129600000000000.003200000000000000.0336025600000000迭代次数40.9222400000000000.0006399999999999990.0120936448000000迭代次数51.046656000000000.0001280000000000000.00435358105600001迭代次数60.9720064000000002.55999999999999e-50.00156728393728001迭代次数71.016796160000005.11999999999999e-60.000564222007705597迭代次数80.9899223040000001.02400000000000e-60.000203119914385407迭代次数91.006046617600002.04799999999999e-77.31231688432032e-5迭代次数100.9963720294400004.09599999999999e-82.63243407701307e-5迭代次数111.002176782336008.19199999999997e-99.47676267670981e-6迭代次数120.9986939305984001.63839999999999e-93.41163456359394e-6迭代次数131.000783641640963.27679999999999e-101.22818844289275e-6迭代次数140.9995298150154246.55359999999997e-114.42147839441398e-7迭代次数151.000282110990751.31071999999999e-111.59173222198927e-7迭代次数160.9998307334055532.62143999999999e-125.73023599916286e-8迭代次数171.000101559956675.24287999999998e-132.06288495970134e-8迭代次数180.9999390640259991.04857599999999e-137.42638585493563e-9结果分析：可以看出极小值为0，此时(x,y）=(1,0)第二章2.11x03455，n=10000~1之间的随机数。x0n1000，下面进行迭代。5534870mod1000，870/10000.870；x255870850mod1000，850/10000.850；55850750mod1000，750/10000.750；x455750250mod1000，250/10000.250；55250750mod1000，750/10000.750；x55750250mod1000r250/10000.2502。6 62、使用平方取中法给出由2563生成随机数序列的前4个数。解：5263的位数为m=4，平方后为6568969，不够8位，左边补零，2m位数是06568969，中间4位数是5689，5689就是由1234得到的随机数。5689的平方为32364721，正好8位，左边不需要补零，中间4位数是3647，3647就是由5689得到的随机数。3647的平方为13300609，正好8位，左边不需要补零，中间4位数是3006，3006就是由3647得到的随机数。3006平方为9036036，不够8位，左边补零，为09036036，中间4位数是0360，360就是由3006得到的随机数。4个随机数为5689,3647,3006,360。3、初值x0为11，乘因子a为91，c为4，n为131。使用混合同余法求5个0-1之间的随机数。解：x011，a91,c4,n131，下面进行迭代。91114100588mod131，88/1310.672；x291884801221mod131，21/1310.160；91214191581mod131，81/1310.618；x491814737539mod131，39/1310.298；x591394355316mod131，r516/1310.122。4、n=5时，求有限域的所有的本原元素及非零元素的阶。解：243442232个本原元素。习题2.21f(x1n

1xn

,0x1

的随机数的抽样公式。解：XBeta(1,1)，其分布函数为n

F(x)

x111xn0n

1dxxn,0x1，1x<010x1F(x)XnXYnY服从（0，1）上的均匀分布，则X服从Beta(1,1)。因此抽样公式可取为XYn。n2、使用逆变换法给出生成密度函数为f(x)nxn1,0x1的随机数的抽样公式。解：随机变量X服从Beta(n,1)，其分布函数为F(x)xnxn1dxxn,0x1，01x<010x1F(x)XnXYn1YXBeta(n,1)XYn。3、使用逆变换法给出生成密度函数为f(x)2x,0x1的随机数的抽样公式。解：随机变量X的密度函数为f(x)2x,0x1，即Beta(2，1)分布，其分布函F(xx2xdxx2,0x1。因此YF(xX2X0

Y。当Y服从（0，1）上的均匀分布，则Xf(x2x,0x1

。因此抽样公式可取为X Y。4、f(x2(1x0x1的随机数的抽样公式。解：随机变量X的密度函数为f(x)2(1x),0x1，即Beta(1，2)分布，其分布函数为

F(x)x2(1x)dx1(1x)2,0x10

。因此0x1时，YF(x)1(1X)2,X1

1YYX的密f(x)2(1x0x1

X1

1Y。5、使用逆变换法给出生成密度函数为f(x)

2 ,0x1的随机数的抽样公1x2解：随机变量X的密度函数为f(x)

2 ,0x1，其分布函数为1x2F(x)x 2 dx2arcsinx0x1x<0时，F(x)=0x11x21x201x2因此0x1时，YF(x)2arcsinxXsinyY（0，1）上的均匀分 2布则X的密度函数为f(x)2 ,0x1因此抽样公式可取为Xsiny。11x26、f(x式。

1(1x2)

2,xR的随机数的抽样公解：随机变量X的密度函数为f(x)

1(1x2)

,xR

，其分布函数为xF(x) 1 dx1arctanx1xR。因此x1x2) 2YF(x)1arctanx1,Xtan(Y)。当Y服从（0，1）上的均匀分布，则X 2 2的密度函数为f(x)1 ,xR(1x2)

。因此抽样公式可取为Xtan(Y)。27f(x)cosx0x的随机数的抽样公式。2解：随机变量X的密度函数为f(x)cosx,0x2

，其分布函数为F(x)xcosxdxsinx0x0xF(x)sinxXarcsiny。0 2 2当Y服从（0，1）上的均匀分布，则X的密度函数为f(x)cosx,0x 。因2此抽样公式可取为Xarcsiny。8、使用逆变换法给出生成密度函数为f(x)

2x/m,0xm

的随机数的抽样公式。

2(1x)/(1m),mx1Xf(x)

2x/m,0xm

，即三角分布，其

2(1x)/(1m),mx1x02 x2xdxx,0xmF(x

0m m 。2m2xdxx21xdx11x)

,mx10m

m1m1,

1mx1YF(x)

2x,0xmxm2

,X

Y0Ym

Y（0，1(1x)

,mx1

1

(1m)(1Y),mY1 1m1）上的均匀分布，则X服从三角分布。9、使用舍选法生成密度函数为f(x)20x(1x)3,0x1的随机数。解：使用求导数的方法，可求出密度函数的最大值为M=f(1)135。4 64伪代码为：While(Y小于等于f(x)=20x(1-x)3)1 2{UUU1 2XU,Y135U1Z=X}

64 2即可求出Z。平均迭代轮数为135/64。10、使用舍选法生成密度函数为 1 1

的随机数。3p(x)3

()2

x2ex,x0f(x)

1e2

x/2,x

0，

h(x)

1x2e

x

0，0

h(x)1，

1 1 2 4L( f(x)h(x)dx)1(

exx2dx)1

>1。3 023

()2计算过程中用到伽玛函数x1exdx)。0当U1UX2lnU1f(x)1/2的指数分布。xWhile(Y大于 xe2)1 2{UU(0,1),UU(0,1),1 2X2lnU1,YU2，if Y小于等于跳出}

xxe2，Z=X即可求出Z，效率为1 。L 411、使用舍选法生成密度函数为p(x)1xexe5,x5的随机数。6f(x1ex/2x5，2

h(x)

xx e2,x2

5，0h(x)1，当x5 1e dxe ，因此f(x)1e e ,x5， x/2 5/2 5/2x/2时，52 2L(f(xh(x)dx)1(xe5/2exdx)12e5/2>1。 54 3因为YF(x)x1e5/2ex/2dx1e(5x)/2,x5，X52ln(1Y)，而当1-Y服52X52lnY当U1UX52lnU1f(x。While(Y大于

xxe2)21 2{UU(0,1),UU(0,1),1 2X52lnU1,YU2，if Y小于等于跳出}

xxe2，2Z=X即可求出Z。效率为131。L 2e5/212、已知Z的密度函数为p(z)1ez,zR，使用复合法生成服从该密度函数的2随机数。f(x1ex，x0,f(x1exx0p(x)1exxRf(x和1 2 2 2 2 1f2(xy<0.5时，Xln1(xy0.5(01)2Xln(1f2(xX。伪代码如下：R1U(0,1)yR1{ify<0.5:R2U(0,1)Xln2R2else:R3U(0,1)3Xln(12R)}3Z=X即可求出Z。习题2.31、取值为有限个的离散型分布的随机数离散型随机变量X的取值为1，2，3，4，5，6，对应的概率分别为0.03,0.07,0.02,0.08,0.5,0.3，由该分布生成100个随机数。XP{XP{X0.07,P{X0.02,P{X4}0.08,P{X5}0.5,P{X6}0.3，分布函数为0, x1F(x)P{Xk},k2,,m

0.03,1x20.1,2x3，即F(x)0.12,3x40.2,4x50.7,5x6x6由定理3知，F(k1)YF(k)，此时X取值X取值如下：Y0.032,0.03Y0.13,0.1Y0.12X4,0.12Y0.25,0.2Y0.76,0.7Y1改进：为了减少判断次数，将概率按照从大到小的顺序排列：0.5，0.3，0.08，0.07，0.03，0.02，相应的X的取值顺序也改为5，6，4，2，1，3。此时分布函数F的取值也相应的改为0，0.5，0.8，0.88，0.95，0.98，1。取Y为服从U（0，1）的随机数，由定理3知，F(k1)YF(k)，此时X取值如下：Y0.56,0.5Y0.84,0.8Y0.88X2,0.88Y0.951,0.95Y0.983,0.98Y1（2）实验步骤：未改进的算法：1）按照概率，计算分布函数F的值；2）生成服从U(0,1)的随机数Y；3）判断F(i1)YF(i)是否成立，如果成立的话，X=X[i-1]。改进后的算法：1）pX的取值也按照概率排序；2）按照排序后的概率，计算分布函数F的值；3）生成服从U(0,1)的随机数Y；4）判断F(i1)YF(i)是否成立，如果成立的话，X=X[i-1]。(3)输出结果：1）未改进的输出结果：产生的取值为有限个值的离散型分布随机数为:[5,2,5,4,5,5,2,5,4,6,5,6,5,5,5,5,5,6,5,5,6,2,6,6,5,6,1,6,3,6,5,6,2,5,4,6,5,6,5,2,6,5,5,5,6,5,3,5,5,6,5,5,6,5,5,5,5,5,5,5,5,5,5,5,5,4,5,5,1,6,6,5,6,2,5,6,6,3,6,5,5,5,5,5,5,5,6,1,5,6,5,1,4,5,6,5,5,5,6,4]2）改进后的输出结果：产生的取值为有限个值的离散型分布随机数为:[5,6,5,6,4,5,5,5,5,5,5,5,5,5,6,5,5,5,5,1,5,6,4,5,5,5,6,4,6,5,1,6,4,6,6,5,6,4,5,5,6,5,5,6,5,5,2,6,5,6,6,6,5,2,4,2,5,6,2,2,6,6,6,1,5,2,5,2,5,5,5,5,5,5,5,5,5,5,5,5,2,6,6,5,6,5,4,4,6,5,3,5,5,6,5,5,5,6,5,4]2、随机排列（不放回）对自然排列1，2，……，m进行随机排列。（1）原理在1，2，……，m随机选择一个数，记为k1，将其和m进行对换；再从1，2，……，m-1随机选择一个数，记为k2，将其和倒数第2个数进行对换，依次进行下去可得到1，2，……，m的一个随机排列。（2）实验步骤1）将1，2，……，m放入列表z中。2）k=，生成服从区间(0,mysy1，此时s1msk1hx中。xsmmk1。3k=-1y，sy1s1到m-1sk2hx中。x存储zsm-1m-1个数值为本次取到的k2。4）依次进行下去，直到取完最后一个数，此时列表x中存放的就是取出的1，2，……，m的一个随机排列。5）如果要取多个随机排列，可设参数n表示取的随机排列个数，多次取随机排列即可。（3）输出结果为：产生的随机排列为[1,6,4,3,10,8,7,2,9,5]如果只想取r个，可以把k的范围限定为k>m-r即可。3、有放回的随机排列对自然排列1，2，……，m进行r次有放回的随机排列。（1）原理随机选择k21，2，……，m的一个随机排列。（2）实验步骤1）将1，2，……，m放入列表z中。2）生成服从区间(0,)ysy1，此时实现s在1到m中随机取值。将本次取到的第s个数值k1记为h，放到空列表x中。x存储取出的数。3rx12M中有放回的取r个的随机排列。（3）输出结果：产生的随机排列为[8,9,9,10,2,5]4、生成服从负二项分布的随机数有一系列独立试验，事件A成功的概率都是p，事件A第r次成功的试验次数kXP(YmCrpr(1p)kr(krr1（帕斯卡分布XNb(r,pk（1）原理mm若1,,mp的几何分布且相互独立，那么Yi服从负二项i1m分布。所以要生成服从负二项分布的随机数，可以先生成相互独立的随机数m1,,nGpXiX服从负二项分布。i1（2）实验步骤1）生成服从U(0,1)的随机数Y；2）X

lnY

，则可以得到服从几何分布的随机数X。 ln(1p) 3）再将X累加，即可得到负二项分布。（3）输出结果：产生服从负二项分布的随机数：[26,30,21,24,20,26,27,34,36,30,28,38,27,28,23,31,30,37,31,22,41,32,29,28,27,30,27,24,41,19,20,29,37,25,35,22,23,22,31,32,28,31,32,30,33,46,41,37,39,29,40,21,35,23,35,35,17,35,33,24,31,22,27,21,34,23,32,30,22,28,29,34,22,35,25,38,22,28,27,24,26,26,33,24,31,25,40,24,28,41,33,36,34,29,32,22,21,34,14,27]（4）Python原有库代码：importnumpyprint(numpy.random.negative_binomial(12,0.4,100))习题2.51、服从F分布的随机数的生成XY2(m2(n

ZX/m为F(m,n)Y/n（1）实验步骤：1）X2(mY2(n。2）令ZX/m，由此得到服从参数为m和n的F分布的随机数。Y/n（2）输出结果：Python原有库代码：importnumpyprint(numpy.random.f(13,21,10))2、服从t分布的随机数的生成XY2(mt（1）实验步骤：1)生成两个均匀分布U(0，1)的U1和U2。

XY/m

为t(m)分布。2)利用U1和U2生成相互独立的服从标准正态分布的随机数X和卡方分布的随机数Y。3)t

XY/m

，可得到t分布的随机数。（2）输出结果（3）Python原有库代码：importnumpyasnpprint(np.random.standard_t(13,10))3、服从对数正态分布的随机数的生成X~N(,2，则YeX的密度函数为Yf(y)Y

1 (lny)22ye 22y0

y0y0Y服从对数正态分布，记为Y~LN(,2。实验步骤：1)生成两个均匀分布U(0，1)的U1和U2。2)利用U1和U2生成服从标准正态分布的随机数Y。3）由XeY，生成服从对数正态分布的随机数。4、服从柯西分布的随机数的生成Xf(x

1(1x2)

,x(,)，其分布函数为F(x)1arctanx1。使用逆变换法可以生成出服从柯西分布的随机数。 2（1）逆变换法随机变量X的密度函数为f(x

1(1x2)

,xR

，其分布函数为xF(x) 1 dx1arctanx1xR。因此x1x2) 2YF(x)1arctanx1,Xtan(Y)。当Y服从（0，1）上的均匀分布，则X 2 2f(x实验步骤：

1(1x2)

,xR

。因此抽样公式可取为Xtan(Y)。2（1）生成YU(0,1)；（2）由抽样公式Xtan(Y)，得到服从分布为柯西分布的随机数。2输出结果：Python原有库代码：importnumpyasnpprint(np.random.standard_cauchy(10))（2）其他变换实验步骤：（1）生成独立的服从标准正态分布的随机数X1,X2。（2）X

X1，服从柯西分布的随机数。X25、使用逆变换法生成密度函数为

f(x)

111xnn

,0x1

的随机数。随机变量X服从Beta(

1,1)F(xx1

1xn

1dxxn,0x1，n 0n1当x<0时，F(x)=0，x1时，F(x)=1。因此0x1时，YF(x)Xn,XYn。当Y服从（0，1）上的均匀分布，则X服从Beta(1,1)。因此抽样公式可取为XYn。n6、使用逆变换法生成密度函数为f(x)2x,0x1的随机数。Y解随机变量X的密度函数为f(x)2x,0x1即Beta(2，1)分布其分布函数为F(x)x2xdxx2,0x1。因此0x1时，YF(x)X2,X 。Y0Y服从（0，1）Xf(x)2x,0x1

。因此抽样公式可取为X Y。7、使用逆变换法生成密度函数为f(x)2(1x),0x1的随机数。解：随机变量X的密度函数为f(x)2(1x),0x1，即Beta(1，2)分布，X1

1Y

F(x)x2(1x)dx1(1x)2,0x1。因此0x101Y时，YF(x)1(1X)2,X1Y

。当Y服从（0，1）上的均匀分布，则Xf(x)2(1x0x1

。因此抽样公式可取为。1x28、1x2

2 ,0x1的随机数。1x2解：随机变量X1x2

2 ,0x1，其分布函数为F(x)x 2 dx2arcsinx0x1x<0时，F(x)=0x1时，F(x)=1。1x201x21 y因此0x1YF(x)XnXsin

。当Y服从（0，1）上的均匀分布，2则X的密度函数为f(x)2 ,0x1。因此抽样公式可取为Xsiny。1x1x29、使用逆变换法生成密度函数为f(x)cosx,0x的随机数。2解：随机变量X的密度函数为f(x)cosx,0x2

，其分布函数为F(x)xcosxdxsinx0x0xF(x)sinxXarcsiny。0 2 2当Y服从（0，1）上的均匀分布，则X的密度函数为f(x)cosx,0x 。因2此抽样公式可取为Xarcsiny。10f(x)

2x/m,0xm

的随机数，此m=0.6。

2(1x)/(1m),mx1解：随机变量X的密度函数为f(x)

2x/m,0xm

，其分布函数为2(1x)/(1m),mx1

x02 x2xdxx,0xmF(x)

0m m 。2m2xdxx21xdx11x)

,mx10m

m1m1,

1mx1YF(x)

2x,0xmxm2

,X

Y0Ym

。当Y服从（0，1(1x)

,mx1

1

(1m)(1Y),mY1 1m1）上的均匀分布，则X服从三角分布。取m=0.6，则X1

0.Y0Y0.6(0.4(1Y),0.6Y111、使用舍选法生成密度函数为f(x)20x(1x)3,0x1的随机数。解：While(Y小于等于p(x)=20x(1-x)3)1 2{UU(0,1),UU(0,1),1 2XU,Y135U1Z=X}

64 2即可求出Z。12、使用舍选法生成密度函数为

1 1

的随机数。解：While(Y大于

xxe2)

f(x)

3()23

x2ex,x01 2{UU(0,1),UU(0,1),1 2X2lnR1,YR2，if Y小于等于跳出}

xxe2，Z=X即可求出Z。效率为1L 413、使用舍选法生成密度函数为f(x)1xexe5,x5的随机数。6解：While(Y大于

xxe2)21 2{UU(0,1),UU(0,1),1 2X52lnR1,YR2，if Y小于等于跳出}

xxe2，2Z=XZ11。L e5/214、已知Z的密度函数为p(z)0.5z,0z1，使用复合法生成服从该密度函数的随机数。解：R1U(0,1){ify<0.5:XU(0,1)Else:R2U(0,1)R2X R2Z=X即可求出Z。15、已知Z的密度函数为p(z)1ez,zR，使用复合法生成服从该密度函数2的随机数。解：R1U(0,1){ify<0.5:R2U(0,1)Xln2R2Else:U3Xln(2(1R))}3Z=X即可求出Z。16、已知Z服从标准正态分布，使用舍选法生成服从该密度函数的随机数。解：YZ的密度函数为(y1)2

2y2e2，y0。先求Y的随机数。f(y)ey,y0，h(y)e

2，0h(y)1，L(

f(yh(y)dy)1(

eye

(y2

dx)1(

1e2e

y22dy)1

12e222p(y)2

y2e2

(y2e 22

0 01eze2Lh(y)f(y)当R1U(0,1)时，XlnR1的密度函数为f(x)，即参数为1的指数分布。当(x1)2U时，Y的密度函数为U(0,1)，Yh(x)e 2，两边取对数，lnY(x，(xlnY

lnY

。利用舍选法生成Z的随机数。若要生成Z，则2 2U(0,1)0.50.5(x1)2While(Y小于 2 )，{RU

U(0,1),1 2XlnR1,YlnR2，(x1)2if Y大于等于 2 ，跳出S=XR3U(0,1)ifR30.5Z=Selse:Z=-S}2e即可求出Z。效率为2e1、计算I1sinxdx。

习题3.10x解：p(x)sinx，设X服从[0，1]上的均匀分布，密度函数为f(x)1,0x1，x则此时I1sinxdxE(sinx)，所以I1

sinxi。0x x输出结果为：

ni1xi0.94588560030144252、求x24xsinxdx0先计算精确值。x24xsindx1x34xosx4sinx)13≈22.9017956。0 3 0 3输出结果：simulatevalue=22.89940087558649Exactvalue=22.901795645549043此时模拟值和精确值近似相等。3、使用随机投点法、平均值法和重要抽样法分别计算I1exdx的值，并比较方0差的大小。解：I1exdxe1≈1.71828。0（1）随机投点法。X服从f(x)0x1，y

f(x)在[0，1]上有界，0yf(x)e。正方形区域为{0x1,0ye}，区域D为D{0X1，Yf(X)}。算法：1）取服从[0，1]上的均匀分布的随机数X和区间[0，e]上的均匀分布的随机数Y；2）如果YeX，计数s加1。3）将前两步重复n次，把s/n作为定积分的近似值。pIe1

ˆeˆ的期望和方差：e eE(ˆ)E(ˆe)pe， ˆ

1 1e2p(1p) e1 1.71828。1输出结果为：

D(1)D(ˆ) n n n模拟值=1.715801844099469精确值=1.718281828459045（2）平均值法U服从[01]fu),0u1XU，1n1un所以I2ei。ni1D(ex)1e2xdx(1exdx)21e2e3，0 0 2 22)D(12

ei)

1 1 2e

3）=0.2420ne2（ni1 n 2 2 nne2（输出结果为：模拟值=1.7186132647133805精确值=1.718281828459045模拟值的标准差=0.4930330927479177（3）重要抽样法f(x)

x x2

，取g(x)1x，但是该函数在区间[0，1]上不能作e 1 x 2!g(x)2(1xyx2(1x)dx2(x1x23 03 3 23y3y1

1为抽样公式，此时h(x)

f(x)g(x)

ex2(1x)3

。步骤如下：（1g(x)21x)g(x)G(x)。33y1（2）生成n个服从[0,1]上的均匀分布的随机数y1,y2,3y1（3）利用抽样公式xx1,x2,,xn。

1求出n个服从密度函数为g(x)的随机数h(x)

f(xi),i2,,n

ˆI1

nnh(x)n计算 i

g(xi)

的平均值3

i1

i，作为积分的估计值。1D(h(x))1

e2x

dx(1exdx)20.0269，02(1x) 0n3nn)D(1n

h(x))0.02692 i 。i1n输出结果为：n模拟值=1.7169256370081443真实值=1.718281828459045模拟值的标准差=0.163768590171378三种方法中重要抽样法的方差最小，平均值法居中，随机投点法的方差最大。x214、使用平均值法求01

1 2e2dx。2101

1 x2e2dx=(1)(0)。使用泰勒展开方法求22t-t

x2 12e2dx 20

1 x22e2dx，将密度函数泰勒展开：21 x2 1

x2 x4 x6

x2n e2= (1 (1)n )2 2

2 2!22

3!23

n!2n1xx2 2 4 6 2n x21xt 2 1 t

x x x n 202

e dx=

)2)

2 2!223!23

en!2n

)dx (t (t

n t2n+12 23 225 237 2n(2n1)取前面的50项，进行计算即可。算法：（1）当0x3时，(x)的幂级数的通项表达式为x2k1

x2k1(2k1)!!1

，所以递推公x2式为：u00，uk uk2,n再计算(x)= e

2，最后计(2k1)!! 算(x)1+(x)u。2 n（2）当x3时，(x)1 1

x2 e e

11335(1)k(2k1)!!)x x3x5 x7

x2k11 x2

1 1 3 35 (2k1)!!1 e2(1 (1)k )x x2x4 x6

x2k(x)的幂级数的通项表达式为(1)k(2k1)!!，所以递推公式为：x2k1，u

(1)k(2k1)!!u

，k2,3,,n。再计算1

x2(x)= (x)= e

和(x)k=1(x)u。x

x2k

k（3）x0(x)1(x。输出结果：分布函数近似值=0.34134474897989975模拟值=0.34101616108361793精确值=0.3413447460685429第一行是使用泰勒展开求得的积分值，第二行结果是使用平均值法求得的积5、使用随机投点法、平均值法和重要抽样法分别计算I1exdx的值，并比较0方差的大小。解：I1exdx1e1=0.63210（1X服从[01f(x),0x1，yf(x在上有界，0yf(x1x0yDDXfX)}。算法：1）取服从[0，1]上的均匀分布的随机数X，Y；2）如果YeX，计数s加1。3）将前两步重复n次，把s/n作为定积分的近似值。pII，， ˆ

p(1p) e1(1e1) 0.23251输出结果为：

D)D(ˆ) n n n模拟值=0.6323333981629577精确值=0.6321205588285577（2）平均值法解：设U服从[0，1]上的均匀分布，密度函数为f(u)1,0u1，则此时1n unXeU，所以I2ei。ni1D(ex)E(e2x)(E(ex))21(1e2)(1e1)212e3e22 2 2D(ˆI)1D(ex)112e3e2=0.0328n n 2 2 n 输出结果为：模拟值=0.6323333981629577精确值=0.6321205588285577模拟值的标准差=0.49199672431939023（3）重要抽样法f(x)

xx2

，取g(x)1x，但是在区间[0，1]上不能作为密度e 1x 11y函数，经过计算，g(x)2(1x)。yx2(1x)dx(2xx2),x0

1，反函x1

（x01，此时h(x)f(x)

ex。算法：

1yg(x) 2(11y1y1）g(x)2(1xg(x1y1y2）成n个服从[0,1]上的均匀分布的随机数y1,y2,,1y

1。x1

求出n个服从密度函数为g(x)的随机数x1,x2,,xn。h(x)

f(xi),i2,,n

ˆI1

nh(x)4）计算 i

g(xi)

的平均值3

i1

i，作为积分的估计值。n1e2xn

1x 2

ˆ 1 2.191D(h(x)021x)dx(0edx)输出结果为：

2.191D(I2)D(h(i)) 。i1nnn模拟值=0.6322688592391777真实值=0.6321205588285577nnn模拟值的标准差=0.5210613930977748216f(x)5(x10x0.5，使用分层抽样法求I141x,0.5x1区间上取的抽样点数相同，且使用平均值法。

0f(xdx把区间[0,1]分成2个区间，[0,1/2]，[1/2,1]，若一共取10个点的话，分别在1525个服从各自区间上的均匀分布的随X1，X2。1I1/22(x1)dx

4(1x)dxII05

/2 1 21/22 1 1

2 1 1

52 1 510

(x1)dx( 0)5 2 n

5(i) 5(i)25(i1 11

1i1n2

5 21 4

i1

i125I2

4(1x)dx(1)

4(1) (1)

(1xi)1/2 2E(I)E(I1)E(I2)I

n2i1

5 2 i1

5i1D(I)D(15

(x1))5D(X

1)1D(X)1 25

ii1

625 1

125 11(/22x2dx(/22dx)2)11125 0 0 60001D(I

)D(251x))4D1X2 i2i

)4D(X

)4( 2x2dx(

12xdx)2)15i1 5

2 5 2

5/2

/2 60D(I)D(I1)D(I2)=0.01683算法：一共取10个点，（1）生成服从均匀分布U(0,1)的随机数U1共5个，Y1=0.5×U1服从均匀分布U(0,0.5)。（2）U(0,1)U25个，Y2=0.5+0.5×U2服从均匀分U(0.5,1)。2(Y1),i1（3）)5 i

，利用平均值法求出各个子区间上积分的平均值，作i 4(1Y),i2iIi2个估计值相加得到积分的模拟值。输出结果：模拟值=0.7471617805512791精确值=0.75本例采用平均值法求分层抽样法，也可使用重要抽样法求分层抽样法。7、使用投点法求无理数e的近似值。解：1exdxe1，长方形：0<x<1，0<y<e，曲线yex,x0,x1,x轴围成0的区域D。向长方形内随机投点，计算点落在区域D内的个数，由几何概型，e1n

p，e

1 e。e N输出结果：

1p2.72168091013009672.718281828459045分析：近似值为2.7216809101300967，精确值为2.718281828459045习题3.21、使用平均值法和控制变量法计算I1exdx的值，并比较两种方法随机误差0的方差的大小。解：精确解：I1exdxe101n1Un设U服从区间[0,1]上的均匀分布，令XeU，使用平均值法，令I1ei，ni1此时nnD(I)D(1nn

Ui)1D(U)1

12u u

21(1e22e3)0.2421 i1

n n(0edu

(0edu)) n 2 2 n ，1其中D(X)D(eU)0.242。1令YU0.5，E(Y)0,D(Y)1/12，Cov(X,Y)E(Y)E(X)EY)E(Y)1euu0.du31e0.14086，0 2 2aCov(X,Y)1.6903时，ZX1.6903YeU1.6903(U0.5)，D(Y)D(X)D(Y)D(X)D(Y)

0.9919，D(Z12)DX10.991920.2420.003872DX0.242。nˆI1n

nZ1

n(X1.6903Y)

)D(1

nZ)0.003872积分的估计值2

ii1

i1

i i， 2

i 。i1nnn使用控制变量法的误差方差小于使用平均值法的误差方差。控制变量法的算法为：nnn（1）生成服从U(0,1)的随机数U，令XeU，YU0.5，则ZX1.6903Y；（2）使用平均值法求Z的平均值即可求出积分的估计值。输出结果为：模拟值=1.71636533513561精确值=1.718281828459045使用平均值法的标准差=0.4941664070633802模拟值=1.7182858741359999精确值=1.718281828459045使用控制变量法的标准差=0.062707728651969452、使用平均值法和对偶变量法分别计算I1exdx的值，并比较两种方法随机0误差的方差的大小。解：精确解：I1exdxe10设U服从区间[0,1]上的均匀分布，令XeU，Ye1U，因为X与Y同分布，则E(Y)E(X),D(Y)D(X)。积分的估计值为11nUiUi n U1

1UI3 (eie i)ni1

2 2ni1I3)

1(D(eU)D(e1U)2Cov(eU,e1U))，4nU 12u

1u 2 e2 3，D(e

)0edu(0edu)

2e2 2 Cov(X,Y)E(XY)E(X)E(Y)1eue1udu(1eudu)2e(e1) 0 0)

1(3e25e5)0.0039133 4n

4 2 4 nY（E(Y)＝0，CovX,Y)0（1）生成服从U(0,1)的随机数U，令XeU，Ye1U，则Z(XY)/2（2）使用平均值法求Z的平均值即可求出积分的估计值。输出结果为：模拟值=1.7180476615566254精确值=1.718281828459045使用对偶变量法的标准差=0.06275545194319794习题3.31424次至少得到一个双六点哪一个概率更大？解：令事件A：一颗骰子投4次至少得到一个六点，事件B：两颗骰子投24次至少得到一个双六点。P1

(5)40.5177，因为投24次双骰子不出6现双六点概率为P(B)(35)24，可以得到P(B)1(35)240.4915。36 36输出结果：出现六点的频率p1=0.5142，出现双六点的频率p2=0.4972。投骰子问题模拟中的两次模拟结果中的理论值与模拟值十分接近，p1与0.5177p20.4914p2，符合事实，由此可验证理论计算的正确性。22000.6，运行的时候将要耗费1KW，求工厂需要多少电力能以0.999的概率保证工厂正常生产。200次伯努利实验。X为正常运行的机器数m200P(Xm)Ck(0.6)k(0.4)200k0.999，考虑到n200k0定理近似计算，变成P(Xm)(m119.5)0.999，解得m140.97141。48进行1000次模拟，计算每次运行正常的机器数量，从小到大排列，取第999个值。连续三次运行的结果为：需要的电力为140需要的电力为139需要的电力为142320解：约会问题符合几何概率，由于他们在任意时刻到达的可能性都相等，那XY,则和Y均符合[060]P0.5*40*400.5*40*400.4444。60*60连续三次运行结果为：错过的理论概率:0.4444模拟值:0.4393错过的理论概率:0.4444模拟值:0.4482错过的理论概率:0.4444模拟值:0.44324、星巴克冰粽销售量符合泊松分布，参数7，问店员每月要销售多少盒冰粽才能0.999的概率达到业绩要求。XXP(7可得到m16。思想：

mmP(Xm)k0

7ke7k!

0.999，10009991000次泊松分布999个数。先产生U1eX0；如果整数k满足rr...r

err...rr

，则Xk，X~P()。12 k连续五次输出结果为：销售的盒数15销售的盒数17销售的盒数16销售的盒数16销售的盒数16

12 kk七次模拟结果均与理论值接近，可以验证模拟结果与实际计算值吻合。1习题3.411、求

ex01x2egrateimport quadimportnumpyasnp解：使用egrateimport quadimportnumpyasnp#第1题定积分z1=lambdax:np.exp(-x)/(1+x**2)I1=quad(z1,0,1)print(I1)输出结果为(0.5247971432602705,5.826418717050591e-15) 2、计算21xy0021y2ddy200 3使用python的Scipy库中的dblquad代码：egrateimport dblquad#2yxz2=lambday,x:x*y**2I2=dblquad(z2,0,2,lambdax:0,lambdax:1)print(I2)输出结果：(0.6666666666666667,7.401486830834377e-15)该例题的积分限均为常数，积分区域为矩形。x23、计算

eD:x2y21

2sin(x2y)dxdy解：先将其化为累次积分1x21e2sin(x2y)dxdy dx

x2

x21e2sin(x2y)dy dy1

y2

x2e2sin(x2y)dx 1D:x2y21

x2

1

y2使用python的Scipy库中的dblquad，得：egrateimport dblquadimportnumpyasnp#第3题二重积分2先y后xz31=lambday,x:np.exp(-(x**2)/2)*np.sin(x**2+y)db31=lambdax:np.sqrt(1-x**2)I31=dblquad(z31,-1,1,lambdax:-db31(x),db31)print(I31)#第3题二重积分2先x后yz32=lambdax,y:np.exp(-(x**2)/2)*np.sin(x**2+y)db32=lambday:np.sqrt(1-y**2)I32=dblquad(z32,-1,1,lambday:-db32(y),db32)print(I32)输出结果：(0.5368603826989582,3.6961556038050958e-09)(0.536860382698808,1.3944401189291966e-12)4、计算x2dxdy，D：xy2,yx2D2 y2 1解： x2dxdy dy x2dx dx

xx4x2dy dx x2dy15.11x4 1D

y2

0x 1x2pythonScipydblquad计算积分：代码：egrateimport dblquadimportnumpyasnp#43yxz41=lambday,x:x**2