STATA入门10 随机模拟.docx
《STATA入门10 随机模拟.docx》由会员分享，可在线阅读，更多相关《STATA入门10 随机模拟.docx（22页珍藏版）》请在冰点文库上搜索。
STATA入门10随机模拟
10随机模拟
只要你自己试试模拟随机现象几次，就会加强对概率的了解，比读很多页的数理统计和概率论的文章还有用。
学习模拟，不仅是为了解模拟本身，也是为更了解概率而了解模拟。
10.1伪随机数
生成（0，1）之间均匀分布的伪随机数的函数为uniform（）
diuniform（）
diuniform（）
diuniform（）
每次都得到一个大于零小于1的随机数。
如果要生成一位数的随机数（即0，1，2，3，4，5，6，7，8，9），可以取小数点后第一位数，通常用下面的命令
diint（10*uniform（））
两位随机数（0-99）则取小数点后两位小数，即
diint（100*uniform（））
任意均匀分布随机数（a,b）由下述函数得到
a+（b-a）*uniform（）
任意均匀分布整数随机数（a,b）由下述函数得到
a+int（（b-a）*uniform（））
也可以同时生成多个随机数，然后将该随机数赋给某个变量。
要注意的是，电脑中给出的随机数不是真正的随机数，而是伪随机数，因为它是按照一定的规律生成的。
如果给定基于生成伪随机数的初始数值（即setseed#），则对相同的初始数值，生成的伪随机数序列完全一样。
*============================begin====================================
clear
setobs10
genx1=uniform（）
genx2=uniform（）//注意到x1与x2不一样
setseed1234
geny1=uniform（）
setseed1234
geny2=uniform（）
geny3=uniform（）//注意到y1与y2一样，但均与y3不同
setseed5634
genz1=uniform（）
setseed1234
genz2=uniform（）//注意到z2与y1,y2一样，但z1与z2不同
list
*============================end====================================
10.2简单模拟
利用随机数字表或者电脑软件中的随机数字，来模仿机遇现象，叫模拟（simulation）、
一旦有了可靠的概率模型，模拟是找出复杂事件发生概率的有效工具。
一个事件在重复结果中发生的比例，迟早会接近它的概率，所以模拟可以对概率做适当的估计。
例1：
如何执行模拟
掷一枚硬币10次，结果中会出现至少3个连续正面或者至少3个连续反面的概率是多少？
思考：
（1）猜想这个概率大约是多少？
（2）如何从理论上计算出这个概率？
（3）如何模拟计算这个概率？
第一步：
提出概率模型。
●每一次掷，正面和反面的概率各为0.5
●投掷之间，彼此是独立的。
也就是说，知道某一次掷出的结果，不会改变任何其他次所掷结果的概率。
第二步：
分配随机数字以代表不同的结果。
●随机数字表中的0-9每个数字出现的概率都是0.1
●每个数字模拟掷一次硬币的结果。
●奇数代表正面，偶数代表反面。
第三步：
模拟多次重复。
●生成10个随机数字
●记录开心的事件（至少连续三个正面或反面）是否发生，如果发生，记为1，否则为0
●重复10次（或者100，1000，1000000次），计算概率=开心事件发生/总重复次数。
真正的概率是0.826。
大部分的人认为连续正面或反面不太容易发生。
但模拟结果足以修正我们直觉错误。
*============================begin====================================
captprogramdropseq3
programseq3,rclass//rclass选项表示计算结果将由return返回到r（）
version9
drop_all//清空所有数据，不能用clear
setobs10//将生成10个观察值
tempvarxyz//设定x,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
forvaluesi=3/10{
replace`z’=1if`y’==`y’[_n-1]&`y’==`y’[_n-2]in`i'//连续三个变量相等时z=1
}
sum`z’
returnscalarmax=r（max）//z取1表示高兴的事发生（三连续），否则失败
end
simulatemax=r（max）,reps（10000）nodots:
seq3//重复1万次，取平均结果
sum
*============================end====================================
由于上述命令要不停生成变量，进行变量代换，所以计算速度较慢，如果使用矩阵，则计算速度会大大加快，命令也更简捷。
*MATA
*============================begin====================================
mata
A=uniform（10000,10）:
>0.5//每行为一次试验,每列为某次试验中硬币的正反面结果
B=J（10000,1,0）//记录每次试验的结果,先记为0.若高兴结果出现再改为1
for（j=1;j<=rows（A）;j++）{//第j次试验
for（i=3;i<=cols（A）;i++）{//第j次试验中第i次硬币出现结果
if（A[j,（i-2,i-1,i）]==（0,0,0）|A[j,（i-2,i-1,i）]==（1,1,1））{
B[j,1]=1//若连续为正面或者连续为反面,则记为1
break
}
}
}
mean（B）//求开心事件的次数
end
*============================end====================================
10.3复杂模拟
例2：
我们要女儿
任务：
一对夫妇计划生孩子生到有女儿才停，或者生了三个就停，他们拥有女儿的概率是多大？
思考：
理论上，概率是多少？
第一步：
概率模型
每一个孩子是女孩的概率是0.49，且各个孩子的性别是互相独立的。
第二步：
分配数字
00，01，02，。
。
。
，49=女孩
49，50，51，。
。
。
，99=男孩
第三步：
模拟
从随机数表中生成一对一对的数字，直到有了女儿，或已有3个孩子。
重复100次。
69051648178717648987
男女女女女男女男男男
用数学可以计算出来，有女孩的真正概率是0.867
*============================begin====================================
captprogramdropgirl
programgirl,rclass
drop_all
setobs3
genx=int（100*uniform（））
geny=（x<49）//生女孩y=1,生男孩，y=0
sumy
returnscalarmax=r（max）//若有一个女孩，则max取1，若三个全为男孩，取0
end
simulatemax=r（max）,reps（10000）nodots:
girl
sum
*============================end====================================
更快的模拟方式.
*============================begin====================================
captprogramdropgirl
programgirl
matA=matuniform（1,3）
scalargirl=1
ifA[1,1]>0.49&A[1,2]>0.49&A[1,3]>0.49{
scalargirl=0//若三个全为男孩，生女孩数为0
}
end
simulategirl,reps（10000）nodots:
girl
tab_sim
*============================end====================================
*============================begin====================================
mata
A=uniform（10000,3）:
<0.49
B=J（10000,1,1）
for（i=1;i<=rows（A）;i++）{
if（A[i,.]==（0,0,0））{
B[i,1]=0//若三个全为男孩，生女孩数为0
}
}
mean（B）
end
*============================end====================================
10.4多阶段模拟
例3：
肾脏移植
张三肾脏不行了，正在等待换肾，他的医师提供了和他情况类似的病人资料。
“换肾后的五年存活率：
撑过手术的有0.9，术后存活的人中有0.6移植成功，0.4还得回去洗肾；五年后，有新肾的人70%仍然活着，而洗肾的只有一半仍活着。
”张三希望知道，他能活过五年的概率，然后再决定是否换肾。
思考：
计算理论概率
第一步：
采用概率树将信息组织起来。
第二步：
分配数字
阶段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%仍然活着，而洗肾的只有一半仍活着。
计算换肾的人活过五年的概率。
*/
captprogramdropsurv
programsurv,rclass
drop_all
setobs1
genz=1
genx1=int（10*uniform（））
ifx1==0{
replacez=0
}
else{
genx2=int（10*uniform（））
ifx2<6{
genx3=int（10*uniform（））
ifx3>6{
replacez=0
}
}
else{
genx4=int（10*uniform（））
ifx4>4{
replacez=0
}
}
}
sumz
returnscalarmax=r（max）
end
simulatemax=r（max）,reps（10000）nodots:
surv
sum
*============================end====================================
更简捷的方式
*============================begin====================================
captprogramdropsurv
programsurv
scalarz=1
ifuniform（）<0.1{
scalarz=0
}
else{
ifuniform（）<0.6{
ifuniform（）>=0.7{
scalarz=0
}
}
else{
ifuniform（）>=0.5{
scalarz=0
}
}
}
end
simulatez,reps（10000）nodots:
surv
sum
*============================end====================================
10.5商店案例
任务：
设某种商品每周的需求量X是[10,30]上均匀分布的随机变量，而某店进货数量为[10,30]中的某一整数。
商店每销售一单位商品可获利500元；若供大于求则削价处理，每处理一单位商店亏损100元；若供不应求则可从外部调剂供应，此时每一单位商店仅获利300元。
为使商店所获利润期望值不少于9280元，试确定最少进货量。
解：
由题设，随机变量X的概率密度为
设商店进货量为a则商店利润Z为
故知最小进货量为21个单位可使期望利润大于9280元。
*============================end====================================
captuprogramdropgoods
programgoods
scalarx=10+int（20*uniform（））
if`1'>=x{
scalarz=500*x-100*（`1'-x）
}
else{
scalarz=500*`1'+300*（x-`1'）
}
end
setmoreoff
quietlyforvaluesi=10/30{
simulatez,rep（1000）nodots:
goods`i'
quietlysum
scalarz`i'=r（mean）
}
scalarlist
*============================end====================================
10.6练习
亚洲随机甲虫的命运（Asianstochasticbeetle），这种昆虫的繁殖模式是：
（1）20%的虫还没有生雌幼虫之前就死掉了，30%生1只雌虫，50%生2只雌虫。
（2）个别雌虫的繁殖情况互相独立。
问：
亚洲随机甲虫的前途会：
繁殖很快？
勉强保持数目？
逐渐灭绝？
生日问题：
只要一间屋里有23个人，则至少有两人同一天生日的概率会超过1/2。
试模拟两个班60名学生同一天过生日的概率，并用你们班，或者前几届后几届班组验证之。
古罗马时代的赌博：
掷4块绵羊距骨是古罗马最受欢迎的赌博，掷多次后骨头的四个面挖概率如下：
窄而平的一面0.1宽而凹的一面0.4
宽而凸的一面0.4窄而凹的一面0.1
掷4块距骨最好的结果叫“维纳斯”，这是朝上的四个面都不一样的情况，估计掷出“维纳斯“的概率。
中国会有多少人成为可怜的单身汉？
不少中国人（尤其是农村）有重男轻女思想，他们总是想尽办法生男孩，性别失调可能会导致越来越多的可怜的男性单身汉。
假设所有夫妇都直到生出男孩，或者生够3个孩子才停止生育，则一个家庭平均会有几个孩子？
多少个男孩？
多少个女孩？
商店的平均获利：
一商店经销某种商品每周进货的数量X与顾客对该种商品的需求量Y是相互独立的随机变量，且都服从区间[10,20]上的均匀分布。
商店每售出一单位商品可获得利润1000元，若需求量超过了进货量，商店可从其它商店调剂供应，这时每单位商品可获利润500元，试计算商店经销该种商品所获利润的期望值。
设随机变量Z表示商店每周所获利润，则
10.7附录
亚洲随机甲虫的命运（Asianstochasticbeetle），这种昆虫的繁殖模式是：
（1）20%的虫还没有生雌幼虫之前就死掉了，30%生1只雌虫，50%生2只雌虫。
（2）个别雌虫的繁殖情况互相独立。
问：
亚洲随机甲虫的前途会：
繁殖很快？
勉强保持数目？
逐渐灭绝？
*===========================begin====================================
captureprogdropbb
progbb,rclass
localbeetle=`1'
while`beetle'<500&`beetle'>0{
drop_all
setobs`beetle'
tempvarxy
gen`x'=uniform（）
gen`y'=1
replace`y'=2if`x'<0.5
replace`y'=0if`x'>0.8
quietlysum`y'
localbeetle=r（sum）
}
noisilydisplayaserror`beetle'
end
setmoreoff
quietlybb10
*===========================end====================================
生日问题：
只要一间屋里有23个人，则至少有两人同一天生日的概率会超过1/2。
试模拟两个班60名学生同一天过生日的概率，并用你们班，或者前几届后几届班组验证之。
*===========================begin===================================
captuprogdropbirth
progbirth
drop_all
setobs60
tempvary
gen`y'=int（365*uniform（））
sort`y’
scalarz=0
forvaluesi=1/59{
if`y'[`i']==`y'[`i'+1]{
scalarz=1
continue,break
}
}
end
simulate"birth"z,reps（100）
sum
*===========================end====================================
古罗马时代的赌博：
掷4块绵羊距骨是古罗马最受欢迎的赌博，掷多次后骨头的四个面挖概率如下：
窄而平的一面0.1宽而凹的一面0.4
宽而凸的一面0.4窄而凹的一面0.1
掷4块距骨最好的结果叫“维纳斯”，这是朝上的四个面都不一样的情况，估计掷出“维纳斯“的概率。
*===========================begin===================================
captuprogdropwns
progwns
drop_all
setobs4
tempvarxy
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`x'
scalarz=1
forvaluesi=1/4{
if`x'[`i']==`x'[`i'+1]{
scalarz=0
continue,break
}
}
end
simulate"wns"z,reps（1000）
sum
*===========================end====================================
中国会有多少人成为可怜的单身汉？
不少中国人（尤其是农村）有重男轻女思想，他们总是想尽办法生男孩，性别失调可能会导致越来越多的可怜的男性单身汉。
假设所有夫妇都直到生出男孩，或者生够3个孩子才停止生育，则一个家庭平均会有几个孩子？
多少个男孩？
多少个女孩？
*===========================begin===================================
captuprogdropchild
progchild
drop_all
setobs3
tempvary
gen`y'=int（100*uniform（））
scalarboy=0
scalargirl=0
if`y'[1]<49{
scalargirl=1
if`y'[2]<49{
scalargirl=2//第二个为女孩
if`y'[3]<49{
scalargirl=3//三女
}
else{
scalarboy=1//2女1男
}
}
else{
scalarboy=1//1女1男
}
}
else{
scalarboy=1//0女1男
}
end
simulate"child"boy=boygirl=girl,reps（1000）
sum
*===========================end====================================
如果不是重男轻女，而是生到三个为止，不论是男是女，则
*===========================begin===================================
captuprogdropchild
progchild,rclass
drop_all
setobs3
tempvarybg
gen`y'=int（100*uniform（））
gen`b'=0
gen`g'=0
replace`b'=1if`y'>=49
sum`b'
scalarboy=r（sum）
replace`g'=1if`y'<49
sum`g'
scalargirl=r（sum）
end
simulate"child"boygirl,reps（1000）
sum
*===========================end====================================
商店的平均获利：
一商店经销某种商品每周进货的数量X与顾客对该种商品的需求量Y是相互独立的随机变量，且都服从区间[10,20]上的均匀分布。
商店每售出一单位商品可获得利润1000元，若需求量超过了进货量，商店可从其它商店调剂供应，这时每单位商品可获利润500元，试计算商店经销该种商品所获利润的期望值。
当X和Y取连续变量时，X=[10，20]，Y=[10，20]
*============================begin====================================
captureprogramdroppf
progpf
scalarx=10+10*uniform（）
scalary=10+10*uniform（）
ifx>=y{
scalarz=1000*y
}
else{
scalarz=500*（x+y）
升级会员 