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*uniform())
任意均匀分布随机数（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 seed 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）、
一旦有了可靠的概率模型，模拟是出复杂事件发生概率的有效工具。一个事件在重复结果中发生的比例，迟早会接近它的概率，所以模拟可以对概率做适当的估计。
例1：如何执行模拟
掷一枚硬币10次，结果中会出现至少3个连续正面或者至少3个连续反面的概率是多少？
思考：
（1） 猜想这个概率大约是多少？
（2） 如何从理论上计算出这个概率？
（3） 如何模拟计算这个概率？
第一步：提出概率模型。
● 每一次掷，正面和反面的概率各为0.5
● 投掷之间，彼此是独立的。也就是说，知道某一次掷出的结果，不会改变任何其他次所掷结果的概率。
第二步：分配随机数字以代表不同的结果。
● 随机数字表中的0-9每个数字出现的概率都是0.1
● 每个数字模拟掷一次硬币的结果。
● 奇数代表正面，偶数代表反面。
第三步：模拟多次重复。
● 生成10个随机数字
● 记录开心的事件（至少连续三个正面或反面）是否发生，如果发生，记为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,y,z为临时变量
gen `x’=int(10*uniform())    //产生10个随机变量，可能为0，1，…scala不是内部或外部命令，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 `z’
return scalar max=r(max)        //z取1表示高兴的事发生（三连续），否则失败
end
simulate max=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次。
6905  16  48  17  8717  648987
男女    女    女  女  男女  男男男
用数学可以计算出来，有女孩的真正概率是0.867
*============================begin====================================
capt program drop girl
program girl, rclass
drop _all
set obs 3
gen x=int(100*uniform())
gen y=(x<49)          //生女孩y=1,生男孩，y=0
sum y
return scalar max=r(max)    //若有一个女孩，则max取1，若三个全为男孩，取0
end
simulate max=r(max),reps(10000) nodots:girl
sum
*============================end====================================
更快的模拟方式.
*============================begin====================================
capt program drop girl
program girl
mat A=matuniform(1,3)
scalar girl=1
if A[1,1]>0.49 & A[1,2]>0.49 & A[1,3]>0.49 {
scalar girl=0    //若三个全为男孩，生女孩数为0
}
end
simulate girl,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=死亡
第三步：模拟

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。