STATA入门10随机模拟.pdf

《STATA入门10随机模拟.pdf》由会员分享，可在线阅读，更多相关《STATA入门10随机模拟.pdf（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、10 随机模拟只要你自己试试模拟随机现象几次，就会加强对概率的了解，比读很多页的数理统计和概率论的文章还有用。学习模拟，不仅是为了解模拟本身，也是为更了解概率而了解模拟。10.1伪随机数生成（ 0，1）之间均匀分布的伪随机数的函数为uniform() di uniform() di uniform() di uniform() 每次都得到一个大于零小于1 的随机数。如果要生成一位数的随机数（即0，1，2，3，4，5，6，7，8，9） ，可以取小数点后第一位数，通常用下面的命令di int(10*uniform() 两位随机数 (0-99)则取小数点后两位小数，即di int(100*unifo

2、rm() 任意均匀分布随机数（a,b）由下述函数得到a+(b-a)*uniform() 任意均匀分布整数随机数（a,b）由下述函数得到a+int(b-a)*uniform()也可以同时生成多个随机数，然后将该随机数赋给某个变量。要注意的是， 电脑中给出的随机数不是真正的随机数，而是伪随机数， 因为它是按照一定的规律生成的。如果给定基于生成伪随机数的初始数值（即set seed # ） ，则对相同的初始数值，生成的伪随机数序列完全一样。*=begin=clear set obs 10 gen x1=uniform() gen x2=uniform() /注意到 x1 与 x2 不一样set se

3、ed 1234 gen y1=uniform() set seed 1234 gen y2=uniform() gen y3=uniform() /注意到 y1 与 y2 一样，但均与y3 不同set seed 5634 gen z1=uniform() set seed 1234 gen z2=uniform()/注意到 z2 与 y1,y2 一样，但z1 与 z2 不同list *=end=10.2简单模拟利用随机数字表或者电脑软件中的随机数字，来模仿机遇现象，叫模拟（simulation ） 、一旦有了可靠的概率模型，模拟是找出复杂事件发生概率的有效工具。一个事件在重复结果中发生的比例，

4、迟早会接近它的概率，所以模拟可以对概率做适当的估计。例 1：如何执行模拟掷一枚硬币10 次，结果中会出现至少3 个连续正面或者至少3 个连续反面的概率是多少？思考：（1）猜想这个概率大约是多少？（2）如何从理论上计算出这个概率？（3）如何模拟计算这个概率？第一步：提出概率模型。每一次掷，正面和反面的概率各为0.5 投掷之间， 彼此是独立的。也就是说，知道某一次掷出的结果，不会改变任何其他次所掷结果的概率。第二步：分配随机数字以代表不同的结果。随机数字表中的0-9 每个数字出现的概率都是0.1 每个数字模拟掷一次硬币的结果。奇数代表正面，偶数代表反面。第三步：模拟多次重复。生成 10 个随机数字

5、记录开心的事件（至少连续三个正面或反面）是否发生，如果发生，记为1，否则为0 重复 10 次（或者 100，1000，1000000 次） ，计算概率 =开心事件发生/总重复次数。真正的概率是0.826。大部分的人认为连续正面或反面不太容易发生。但模拟结果足以修正我们直觉错误。*=begin=capt program drop seq3 program seq3,rclass /rclass选项表示计算结果将由return返回到 r() version 9 drop _all /清空所有数据，不能用clear set obs 10/ 将生成 10 个观察值tempvar x y z /设定 x

6、,y,z为临时变量gen x =int(10*uniform() /产生 10 个随机变量，可能为0，1，9 gen y =(mod( x ,2)=0)/ 如果生成的随机变量为奇数，则y=0; 为偶数， y=1 gen z =0/ 生成 Z=0 forvalues i=3/10 replace z =1 if y = y _n -1 & y = y _n -2 in i / 连续三个变量相等时z=1 sum zreturn scalar max=r(max) /z取 1 表示高兴的事发生（三连续），否则失败end simulate max=r(max),reps(10000) nodots:s

7、eq3 / 重复 1 万次，取平均结果sum *=end= 由于上述命令要不停生成变量，进行变量代换，所以计算速度较慢，如果使用矩阵，则计算速度会大大加快，命令也更简捷。*MATA*=begin=mata A=uniform(10000,10):0.5 /每行为一次试验, 每列为某次试验中硬币的正反面结果B=J(10000,1,0) / 记录每次试验的结果, 先记为 0. 若高兴结果出现再改为1for (j=1;j=rows(A);j+) /第 j 次试验for (i=3;i=cols(A);i+) / 第 j 次试验中第i 次硬币出现结果if ( Aj,(i-2,i-1,i)=(0,0,0)

8、 | Aj,(i-2,i-1,i)=(1,1,1) ) Bj,1=1 / 若连续为正面或者连续为反面, 则记为 1break mean(B) / 求开心事件的次数end *=end= 10.3复杂模拟例 2：我们要女儿任务： 一对夫妇计划生孩子生到有女儿才停，或者生了三个就停，他们拥有女儿的概率是多大？思考：理论上，概率是多少？第一步：概率模型每一个孩子是女孩的概率是0.49 ，且各个孩子的性别是互相独立的。第二步：分配数字00 ，01 ，02 ，。， 49= 女孩49 ，50 ，51 ，。， 99= 男孩第三步：模拟从随机数表中生成一对一对的数字，直到有了女儿，或已有3个孩子。重复 100

9、次。6905 16 48 17 8717 648987 男女女女女男女男男男用数学可以计算出来，有女孩的真正概率是0.867 *=begin=capt program drop girl program girl, rclass drop _all set obs 3gen x=int(100*uniform() gen y=(x0.49 & A1,20.49 & A1,30.49 scalar girl=0 /若三个全为男孩，生女孩数为0 end simulate girl ,reps(10000) nodots :girl tab _sim *=end= *=begin=mata A=un

10、iform(10000,3):0.49 B=J(10000,1,1) for (i=1;i=rows(A);i+) if (Ai,.=(0,0,0) Bi,1=0 /若三个全为男孩，生女孩数为0 mean(B) end *=end= 10.4多阶段模拟例 3：肾脏移植张三肾脏不行了，正在等待换肾，他的医师提供了和他情况类似的病人资料。“换肾后的五年存活率：撑过手术的有0.9 ，术后存活的人中有0.6 移植成功， 0.4 还得回去洗肾；五年后，有新肾的人 70% 仍然活着，而洗肾的只有一半仍活着。”张三希望知道，他能活过五年的概率，然后再决定是否换肾。思考：计算理论概率第一步：采用概率树将信息组

11、织起来。第二步：分配数字阶段 1：0= 死亡1-9= 存活阶段 2：0-5= 移植成功6-9= 仍需洗肾阶段 3，成功0-6= 存活五年7-9= 死亡阶段 3：洗肾0-4= 存活五年5-9= 死亡第三步：模拟数学计算结果为0.558。*=begin=/* 换肾后的五年存活率：撑过手术牛0.9 ，术后存活的人中有0.6 移植成功， 0.4 还得回去洗肾；五年后，有新肾的人70%仍然活着，而洗肾的只有一半仍活着。计算换肾的人活过五年的概率。 */ capt program drop survprogram surv,rclass drop _all set obs 1 gen z=1 gen x1

12、=int(10*uniform() if x1=0 replace z=0 else gen x2=int(10*uniform() if x26 replace z=0 else gen x4=int(10*uniform() if x44 replace z=0 sum zreturn scalar max=r(max) end simulate max=r(max),reps(10000) nodots:surv sum *=end= 更简捷的方式*=begin= capt program drop surv program surv scalar z=1 if uniform()0.1

13、 scalar z=0 else if uniform()=0.7 scalar z=0 else if uniform()=0.5 scalar z=0 end simulate z,reps(10000) nodots:surv sum *=end= 10.5商店案例任务：设某种商品每周的需求量X 是 10,30上均匀分布的随机变量，而某店进货数量为10,30中的某一整数。商店每销售一单位商品可获利500 元；若供大于求则削价处理，每处理一单位商店亏损100 元；若供不应求则可从外部调剂供应，此时每一单位商店仅获利300 元。为使商店所获利润期望值不少于9280 元，试确定最少进货量。解：

14、由题设，随机变量X 的概率密度为设商店进货量为a 则商店利润Z 为11030200 xfx其它500100500300XaXZaXa10306001001030020030100XaXaXaaXE ZZ x fx dxdxZ x dxdx故知最小进货量为21 个单位可使期望利润大于9280 元。*=end= captu program drop goods program goodsscalar x=10+int(20*uniform() if 1=x scalar z=500*x-100*(1-x) else scalar z=500*1+300*(x-1) end set more off

15、 quietly forvalues i=10/30 simulate z,rep(1000) nodots: goods i quietly sum scalar zi=r(mean) scalar list *=end= 1002102160010020130020020130010020115020020aaaxa dxax dxxaxxax227.5350525092807.5350403003622.5650262652026332.5aaaaaaa10.6练习亚洲随机甲虫的命运（Asian stochastic beetle） ，这种昆虫的繁殖模式是：（1）20%的虫还没有生雌幼虫

16、之前就死掉了，30%生 1 只雌虫， 50%生 2 只雌虫。（ 2）个别雌虫的繁殖情况互相独立。问：亚洲随机甲虫的前途会：繁殖很快？勉强保持数目？逐渐灭绝？生日问题： 只要一间屋里有23 个人， 则至少有两人同一天生日的概率会超过1/2。试模拟两个班 60 名学生同一天过生日的概率，并用你们班，或者前几届后几届班组验证之。古罗马时代的赌博： 掷 4 块绵羊距骨是古罗马最受欢迎的赌博，掷多次后骨头的四个面挖概率如下：窄而平的一面0.1 宽而凹的一面0.4 宽而凸的一面0.4 窄而凹的一面0.1 掷 4 块距骨最好的结果叫“维纳斯”，这是朝上的四个面都不一样的情况，估计掷出 “维纳斯“的概率。中国

17、会有多少人成为可怜的单身汉？不少中国人（尤其是农村）有重男轻女思想，他们总是想尽办法生男孩，性别失调可能会导致越来越多的可怜的男性单身汉。假设所有夫妇都直到生出男孩， 或者生够3 个孩子才停止生育，则一个家庭平均会有几个孩子？多少个男孩？多少个女孩？商店的平均获利：一商店经销某种商品每周进货的数量X 与顾客对该种商品的需求量Y 是相互独立的随机变量，且都服从区间 10,20 上的均匀分布。 商店每售出一单位商品可获得利润 1000 元，若需求量超过了进货量，商店可从其它商店调剂供应，这时每单位商品可获利润 500 元，试计算商店经销该种商品所获利润的期望值。设随机变量Z 表示商店每周所获利润，

18、则,1020110201000XYfx yfx fyxy其它100010005001000500,YZXYXXYXXYYXE ZZ x fx y dxdy1220202010101020102021011000 *50010010510203510502DDyyydxdyxy dxdydyydxdyxy dxyy dyyYdy202310203210110 10315550220000yyyyy10.7附录亚洲随机甲虫的命运（Asian stochastic beetle） ，这种昆虫的繁殖模式是：（1）20%的虫还没有生雌幼虫之前就死掉了，30%生 1 只雌虫， 50%生 2 只雌虫。（ 2

19、）个别雌虫的繁殖情况互相独立。问：亚洲随机甲虫的前途会：繁殖很快？勉强保持数目？逐渐灭绝？*=begin=capture prog drop bbprog bb,rclass local beetle=1while beetle0 drop _all set obs beetle tempvar x y gen x=uniform() gen y=1 replace y=2 if x0.8 quietly sum y local beetle=r(sum) noisily display as error beetle end set more off quietly bb 10 *=end=

20、 生日问题： 只要一间屋里有23 个人， 则至少有两人同一天生日的概率会超过1/2。试模拟两个班 60 名学生同一天过生日的概率，并用你们班，或者前几届后几届班组验证之。*=begin= captu prog drop birth prog birth drop _all set obs 60 tempvar y gen y= int(365* uniform() sort yscalar z=0 forvalues i=1/59 if yi=yi+1 scalar z=1 continue, break end simulate birth z,reps(100) sum *=end= 古罗

21、马时代的赌博： 掷 4 块绵羊距骨是古罗马最受欢迎的赌博，掷多次后骨头的四个面挖概率如下：窄而平的一面0.1 宽而凹的一面0.4 宽而凸的一面0.4 窄而凹的一面0.1 掷 4 块距骨最好的结果叫“维纳斯”，这是朝上的四个面都不一样的情况，估计掷出 “维纳斯“的概率。*=begin= captu prog drop wnsprog wnsdrop _all set obs 4tempvar x y gen y= uniform()recode y (0/0.1=1) (0.1/0.2=2) (0.2/0.6=3) (0.6/1=4), gen(x)sort xscalar z=1forvalu

22、es i=1/4 if xi=xi+1 scalar z=0 continue, break end simulate wns z,reps(1000)sum *=end=中国会有多少人成为可怜的单身汉？不少中国人（尤其是农村）有重男轻女思想，他们总是想尽办法生男孩，性别失调可能会导致越来越多的可怜的男性单身汉。假设所有夫妇都直到生出男孩， 或者生够3 个孩子才停止生育，则一个家庭平均会有几个孩子？多少个男孩？多少个女孩？*=begin= captu prog drop childprog child drop _all set obs 3 tempvar y gen y=int(100*un

23、iform() scalar boy=0scalar girl=0if y149 scalar girl=1 if y249 scalar girl=2 /第二个为女孩if y3=49 sum bscalar boy=r(sum) replace g=1 if y=y scalar z=1000*y else scalar z=500*(x+y) end simulate pf z,reps(10000) sum *=end= 结果约为： 14153.12 当进货量和销售量均取整数时，即X=10 ，11 ，。 20 ；而 y=10 ，11 ，。， 20 *=begin= capture pro

24、gram drop pf prog pf scalar x=10+int(11*uniform() scalar y=10+int(11*uniform() if x=y scalar z=1000*y else scalar z=500*(x+y) end simulate pf z,reps(10000) sum *=end= 结果约为： 14102.65 clear mat A=J(11,11,.) forvalue i=1/11 forvalue j=1/11 if i=j mat Ai,j=1000*(9+i) else if i=j mat Ai,j=1000*(10+(i-1)/10) else if ij mat Ai,j=500*(20+(i-1)/10+(j-1)/10) mat a=J(1,101,1) mat C=a*A*a/10201 mat list C 15841.584 *=end= 请思考，上述程序有什么问题？为什么答案有如此大的差距？