1、10隨機模擬只要你自己試試模擬隨機現(xiàn)象幾次，就會加強對概率的了解， 比讀很多頁的數(shù)理統(tǒng)計和概率論的文章還有用。學習模擬，不僅是為了解模擬本身，也是為更了解概率而了解模擬。101偽隨機數(shù)生成(0， 1)之間均勻分布的偽隨機數(shù)的函數(shù)為uniform()di uniform。di uniform。di uniform。每次都得到一個大于零小于1的隨機數(shù)。如果要生成一位數(shù)的隨機數(shù)(即 0，1，2，3, 4，5，6，乙8，9)，可以取小數(shù)點后第一位 數(shù)，通常用下面的命令di int(10*uniform()兩位隨機數(shù)(0-99)則取小數(shù)點后兩位小數(shù)，即di int(100*uniform()任意均勻分布

2、隨機數(shù)(a,b)由下述函數(shù)得到a+(b-a)*u ni form()任意均勻分布整數(shù)隨機數(shù)(a,b)由下述函數(shù)得到a+i nt(b-a)*u ni form()也可以同時生成多個隨機數(shù),然后將該隨機數(shù)賦給某個變量。要注意的是，電腦中給出的隨機數(shù)不是真正的隨機數(shù)，而是偽隨機數(shù)，因為它是按照一定的規(guī)律生成的。如果給定基于生成偽隨機數(shù)的初始數(shù)值(即 set seed #)，則對相同的初始數(shù)值，生成的偽隨機數(shù)序列完 全一樣。*=begin=clearset obs10gen x1= un iform() gen x2=un iform() set seed1234gen y1= un iform()

3、set seed1234gen y2=uniform。 gen y3=uniform。 set seed5634gen z1= uniform。set seed1234gen z2=un iform() list注意到x1與x2不一樣注意到y(tǒng)1與y2 樣，但均與 y3不同/注意到z2與y1,y2 樣，但z1與z2不同=end=10.2簡單模擬利用隨機數(shù)字表或者電腦軟件中的隨機數(shù)字，來模仿機遇現(xiàn)象，叫模擬（simulation ）、一旦有了可靠的概率模型，模擬是找出復雜事件發(fā)生概率的有效工具。一個事件在重復結果中發(fā)生的比例，遲早會接近它的概率，所以模擬可以對概率做適當?shù)墓烙嫛＠?:如何執(zhí)行模擬擲

4、一枚硬幣10次，結果中會出現(xiàn)至少 3個連續(xù)正面或者至少 3個連續(xù)反面的概率是多少？ 思考：（1）猜想這個概率大約是多少？（2）如何從理論上計算出這個概率？（3）如何模擬計算這個概率？第一步：提出概率模型。每一次擲，正面和反面的概率各為0.5投擲之間，彼此是獨立的。也就是說，知道某一次擲出的結果，不會改變任何其他次所 擲結果的概率。第二步：分配隨機數(shù)字以代表不同的結果。隨機數(shù)字表中的0-9每個數(shù)字出現(xiàn)的概率都是 0.1每個數(shù)字模擬擲一次硬幣的結果。奇數(shù)代表正面，偶數(shù)代表反面。第三步：模擬多次重復。生成10個隨機數(shù)字記錄開心的事件（至少連續(xù)三個正面或反面）是否發(fā)生，如果發(fā)生，記為1,否則為0重復1

5、0次（或者100，1000，1000000次），計算概率=開心事件發(fā)生/總重復次數(shù)。 真正的概率是 0.826。大部分的人認為連續(xù)正面或反面不太容易發(fā)生。但模擬結果足以修正我們直覺錯誤。*=begin=capt program drop seq3program seq3,rclass/rclass選項表示計算結果將由return返回到r()version 9drop _all/清空所有數(shù)據(jù)，不能用clearset obs10/將生成10個觀察值tempvar x y z/設定x,y,z 為臨時變量gen 'x' =int(10*uniform()/ 產(chǎn)生 10 個隨機變量，可能

6、為 0, 1,，9gen 'y' =(mod*' ,2)=0)/如果生成的隨機變量為奇數(shù)，則y=0;為偶數(shù)，y=1gen 'z' =0/ 生成 Z=0forvalues i=3/10 replace 'z' =1 ify' =y' _-1 & ' y' =y' _2 in、i'/ 連續(xù)三個變量相等時 z=1sum 'z'return scalar max=r(max)/z 取1表示高興的事發(fā)生(三連續(xù))，否則失敗endsimulate max=r(max),reps(

7、10000) nodots:seq3 / 重復 1 萬次，取平均結果sum=end=由于上述命令要不停生成變量，進行變量代換，所以計算速度較慢，如果使用矩陣，則計算速度會大大加快，命令也更簡捷。*MATA=begin= mataA=uniform(10000,10):>0.5/每行為一次試驗，每列為某次試驗中硬幣的正反面結果B=J(10000,1,0)/記錄每次試驗的結果,先記為0.若高興結果出現(xiàn)再改為 1for (j=1;jv=rows(A);j+) / 第 j 次試驗for (i=3;iv=cols(A);i+) /第j次試驗中第i次硬幣出現(xiàn)結果if ( Aj,(i-2,i-1,i)

8、=(0,0,0) | Aj,(i-2,i-1,i)=(1,1,1) )Bj,1=1 II若連續(xù)為正面或者連續(xù)為反面，則記為1breakmean(B) II求開心事件的次數(shù)end=end=10.3復雜模擬例2:我們要女兒任務：一對夫婦計劃生孩子生到有女兒才停，或者生了三個就停，他們擁有女兒的概率是多大？思考：理論上，概率是多少？第一步：概率模型每一個孩子是女孩的概率是0.49，且各個孩子的性別是互相獨立的。第二步：分配數(shù)字00，01， 02 ,。，49=女孩49，50，51 ,。，99=男孩第三步：模擬從隨機數(shù)表中生成一對一對的數(shù)字，直到有了女兒，或已有3個孩子。重復100次。6905 16 4

9、8 17 8717 648987男女 女 女 女 男女 男男男用數(shù)學可以計算出來，有女孩的真正概率是0.867*=begin= capt program drop girlprogram girl, rclassdrop _allset obs3gen x=i nt(100*u niform()gen y=(x<49)/ 生女孩 y=1,生男孩，y=0sum yreturn scalar max=r(max) II若有一個女孩，則max取1，若三個全為男孩，取0endsimulate max=r(max),reps(10000) nodots:girlsum *=end=更快的模擬方式

10、.*=begin=capt program drop girlprogram girlmat A=matuniform (1,3)scalar girl=1if A1,1>0.49 & A1,2>0.49 & A1,3>0.49 scalar girl=0 II 若三個全為男孩，生女孩數(shù)為 0 endsimulate girl,reps(10000) nodots :girltab _sim=end=begin= mataA=uniform(10000,3):<0.49B=J(10000,1,1)for (i=1;i<=rows(A);i+) if

11、 (Ai,.=(0,0,0) Bi,1=0 II 若三個全為男孩，生女孩數(shù)為 0 mean(B) end=end=10.4 多階段模擬例 3：腎臟移植 張三腎臟不行了，正在等待換腎，他的醫(yī)師提供了和他情況類似的病人資料。“ 換腎后的五年存活率：撐過手術的有 0.9 ，術后存活的人中有 0.6移植成功， 0.4 還得回去洗腎；五年后， 有新腎的人 70%仍然活著，而洗腎的只有一半仍活著?！睆埲Ｍ?，他能活過五年的概率， 然后再決定是否換腎。思考：計算理論概率 第一步：采用概率樹將信息組織起來。第二步：分配數(shù)字階段1：0=死亡1-9=存活 階段2 :0-5=移植成功6- 9=仍需洗腎 階段3，

12、成功0-6=存活五年7- 9=死亡 階段3 :洗腎0-4=存活五年5-9=死亡第三步：模擬數(shù)學計算結果為0.558 。*=begin=/*換腎后的五年存活率：撐過手術牛0.9，術后存活的人中有0.6移植成功，0.4還得回去洗腎；五年后，有新腎的人 70%仍然活著，而洗腎的只有一半仍活著。計算換腎的人活過五 年的概率。*/capt program drop surv program surv,rclass drop _all set obslgenz=1gen x1=i nt(10*u niform() if x1=0 replace z=0else gen x2=i nt(10*u nifor

13、m() if x2<6 gen x3=i nt(10*u niform() if x3>6 replace z=0 else gen x4=in t(10* uniform。)if x4>4 replace z=0sum zreturn scalar max=r(max)endsimulate max=r(max),reps(10000) nodots:surv sum*=end= 更簡捷的方式*=begin= capt program drop survprogram survscalar z=1if uniform()<0.1 scalar z=0else if u

14、niform()<0.6 if uniform()>=0.7 scalar z=0else if uniform()>=0.5 scalar z=0endsimulate z,reps(10000) nodots:survsum*=end=10.5商店案例任務：設某種商品每周的需求量X是10,30上均勻分布的隨機變量，而某店進貨數(shù)量為10,30中的某一整數(shù)。商店每銷售一單位商品可獲利500元；若供大于求則削價處理，每處理一單位商店虧損100元；若供不應求則可從外部調劑供應，此時每一單位商店僅獲利300元。為使商店所獲利潤期望值不少于9280元，試確定最少進貨量。解：由題設，隨

15、機變量 X的概率密度為1f x = 20 【010 Ex 乞30其它設商店進貨量為a則商店利潤Z為.500X -100(a-X )Z 、500a+300(X -a )_ 600X -100a10 乞 X a一 300X200a a 空 X 乞 30qQE Zx f x dx1030 1J；0dxZ x dx 0dx20I06°°x1°°a dx1 0300a 200x dx=丄300x2 - 100axF +20 - 10£ |150x2 200ax2-7.5a350a 5250 _ 92802二 7.5a 350a+4030 E03a-62

16、2.5a-65 <02 626520a263 32.5故知最小進貨量為21個單位可使期望利潤大于9280元。*=end=captu program drop goodsprogram goodsscalar x=10+i nt(20*u niform()if、1'>=x scalar z=500*x-100*('1'-x)else scalar z=500*'1'+300*(x-'1')endset more offquietly forvalues i=10/30 simulate z,rep( 1000) nodots:

17、goods 'i'quietly sumscalar z'i'=r(mea n)scalar list=end=10.6練習亞洲隨機甲蟲的命運（Asian stochastic beetle）,這種昆蟲的繁殖模式是：（1） 20%的蟲還沒有生雌幼蟲之前就死掉了，30%生1只雌蟲，50%生2只雌蟲。（2）個別雌蟲的繁殖情況互相獨立。問：亞洲隨機甲蟲的前途會：繁殖很快？勉強保持數(shù)目？逐漸滅絕？生日問題：只要一間屋里有23個人，則至少有兩人同一天生日的概率會超過1/2。試模擬兩個班60名學生同一天過生日的概率，并用你們班，或者前幾屆后幾屆班組驗證之。古羅馬時代的賭博:

18、擲4塊綿羊距骨是古羅馬最受歡迎的賭博,擲多次后骨頭的四個面挖概率如下:窄而平的一面 0.1寬而凹的一面 0.4寬而凸的一面 0.4窄而凹的一面 0.1擲4塊距骨最好的結果叫“維納斯”，這是朝上的四個面都不一樣的情況,估計擲出“維納斯“的概率。中國會有多少人成為可憐的單身漢？不少中國人（尤其是農(nóng)村）有重男輕女思想，他們總是想盡辦法生男孩，性別失調可能會導致越來越多的可憐的男性單身漢。假設所有夫婦都直到生出男孩，或者生夠3個孩子才停止生育，則一個家庭平均會有幾個孩子？多少個男孩？ 多少個女孩？商店的平均獲利： 一商店經(jīng)銷某種商品每周進貨的數(shù)量X與顧客對該種商品的需求量 Y是相互獨立的隨機變量，且都

19、服從區(qū)間10,20上的均勻分布。商店每售出一單位商品可獲得利 潤1000元，若需求量超過了進貨量，商店可從其它商店調劑供應，這時每單位商品可獲利 潤500元，試計算商店經(jīng)銷該種商品所獲利潤的期望值。f x, y = fx x 彳丫 y110 _x _20=210010 豈 y 空 20.0 其它設隨機變量Z表示商店每周所獲利潤，則100CY1000X 500 Y - X1000X Y 乞 X500 X Y Y XE Zx f x,y dxdy=1000y*D11100dxdy 亠 11500 x y dxdyD2202020 y= 1010dy y ydx 10dy105 x ydx20=10

20、 10 y 20-y dy20 3 25-y2 -10Y -50201010213=10 10y -§y5 1 y3 _5y2 _50 yQnnnn10.7 附錄亞洲隨機甲蟲的命運(Asian stochastic beetle),這種昆蟲的繁殖模式是：(1) 20%的蟲還沒有生雌幼蟲之前就死掉了， 30%生 1 只雌蟲， 50%生 2只雌蟲。(2)個別雌蟲的繁殖情況互 相獨立。問：亞洲隨機甲蟲的前途會：繁殖很快？勉強保持數(shù)目？逐漸滅絕？*=begin= capture prog drop bb prog bb,rclasslocal beetle='1'while

21、'beetle'<500 & 'beetle'。drop _allset obs'beetle'tempvar x ygen 'x'=uniform()gen 'y'=1replace 'y'=2 if 'x'<0.5replace 'y'=0 if 'x'>0.8quietly sum 'y'local beetle=r(sum)noisily display as error 'beetle'

22、;endset more offquietly bb 10*=end=生日問題： 只要一間屋里有 23 個人， 則至少有兩人同一天生日的概率會超過1/2。試模擬兩個班 60 名學生同一天過生日的概率，并用你們班，或者前幾屆后幾屆班組驗證之。*=begin=captu prog drop birthprog birthdrop _allset obs 60tempvar ygen 'y'= int(365* uniform。)sort 'y'scalar z=0forvalues i=1/59 if 'y''i'='y

23、9;'i'+1 scalar z=1continue, break endsimulate "birth" z,reps(100)sum=end=古羅馬時代的賭博：擲 4 塊綿羊距骨是古羅馬最受歡迎的賭博，擲多次后骨頭的四個面挖概率如下：窄而平的一面 0.1寬而凸的一面 0.4寬而凹的一面 0.4窄而凹的一面 0.1擲 4 塊距骨最好的結果叫維納斯”，這是朝上的四個面都不一樣的情況，估計擲出“維納斯“的概率。=begincaptu prog drop wnsprog wns drop _all set obs 4tempvar x ygen、y'=

24、uniform。recode、y' (0/0.1=1) (0.1/0.2=2) (0.2/0.6=3) (0.6/1=4), gen('x') sort 'x'scalar z=1forvalues i=1/4 if 'x''i'='x''i'+1 scalar z=0 continue, breakendsimulate "wns" z,reps(1000)sum*=end=中國會有多少人成為可憐的單身漢？ 不少中國人（尤其是農(nóng)村）有重男輕女思想，他們總 是想盡辦法生男

25、孩， 性別失調可能會導致越來越多的可憐的男性單身漢。 假設所有夫婦都直 到生出男孩， 或者生夠 3 個孩子才停止生育， 則一個家庭平均會有幾個孩子？多少個男孩？ 多少個女孩？*=begin=captu prog drop childprog childdrop _allset obs 3tempvar ygen、y'=int(100*uniform()scalar boy=0scalar girl=0if、y'1<49 scalar girl=1if、y'2<49 scalar girl=2 /第二個為女孩if、y'3<49 scalar gi

26、rl=3 /三女else scalar boy=1 /2女 1男else scalar boy=1 /1女1男else scalar boy=1 /0女 1男endsimulate "child" boy=boy girl=girl,reps(1000)sum*=end= 如果不是重男輕女，而是生到三個為止，不論是男是女，則*=begin= captu prog drop child prog child, rclass drop _allset obs 3tempvar y b ggen、y'=int(100*uniform()gen、b'=0gen、g&

27、#39;=0replace、b'=1 if 'y'>=49sum、b'scalar boy=r(sum)replace、g'=1 if 'y'<49sum、g'scalar girl=r(sum)endsimulate "child" boy girl,reps(1000)sum*=end=商店的平均獲利：一商店經(jīng)銷某種商品每周進貨的數(shù)量X與顧客對該種商品的需求量Y是相互獨立的隨機變量，且都服從區(qū)間10,20上的均勻分布。商店每售出一單位商品可獲得利 潤1000元，若需求量超過了進貨量，商店可從其它商店調劑供應，這時每單位商品可獲利 潤500元，試計算商店經(jīng)銷該種商品所獲利潤的期望值。當X和Y取連續(xù)變量時，X=10，20，Y=10，20 *=begin=capture program drop pfprog pfscalar x=10+10*uniform()scalar y=10+10*uniform()if x>=y scalar z=1000*yelse scalar z=500*(x+y)e