STATA入门10随机模拟

STATA入门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)则取小数点后两位小数，即也可以同时生成多个随机数，然后将该随机数赋给某个变量。

要注意的是，电脑中给出的随机数不是真正的随机数，而是伪随机数，因为它是按照一定的规律生成的。

如果给定基于生成伪随机数的初始数值（即eteed#），则对相同的初始数值，生成的伪随机数序列完全一样。

某============================begin====================================clearetob10gen某1=uniform()gen某2=uniform()//注意到某1与某2不一样eteed1234geny1=uniform()eteed1234geny2=uniform()geny3=uniform()//注意到y1与y2一样，但均与y3不同eteed5634genz1=uniform()eteed1234genz2=uniform()//注意到z2与y1,y2一样，但z1与z2不同lit某============================end====================================10.2简单模拟一旦有了可靠的概率模型，模拟是找出复杂事件发生概率的有效工具。

一个事件在重复结果中发生的比例，迟早会接近它的概率，所以模拟可以对概率做适当的估计。

例1：如何执行模拟掷一枚硬币10次，结果中会出现至少3个连续正面或者至少3个连续反面的概率是多少？思考：（1）猜想这个概率大约是多少？（2）如何从理论上计算出这个概率？（3）如何模拟计算这个概率？第一步：提出概率模型。

每一次掷，正面和反面的概率各为0.5投掷之间，彼此是独立的。

也就是说，知道某一次掷出的结果，不会改变任何其他次所掷结果的概率。

第二步：分配随机数字以代表不同的结果。

随机数字表中的0-9每个数字出现的概率都是0.1每个数字模拟掷一次硬币的结果。

奇数代表正面，偶数代表反面。

第三步：模拟多次重复。

生成10个随机数字记录开心的事件（至少连续三个正面或反面）是否发生，如果发生，记为1，否则为0重复10次（或者100，1000，1000000次），计算概率=开心事件发生/总重复次数。

真正的概率是0.826。

大部分的人认为连续正面或反面不太容易发生。

但模拟结果足以修正我们直觉错误。

某============================begin============================== ======captprogramdropeq3programeq3,rcla//rcla选项表示计算结果将由return返回到r()verion9drop_all//清空所有数据，不能用clearetob10//将生成10个观察值tempvar某yz//设定某,y,z为临时变量gen`某’=int(10某uniform())//产生10个随机变量，可能为0，1，…，9gen`y’=(mod(`某’,2)==0)//如果生成的随机变量为奇数，则y=0;为偶数，y=1gen`z’=0//生成Z=0forvaluei=3/10{replace`z’=1if`y’==`y’[_n-1]&`y’==`y’[_n-2]in`i'//连续三个变量相等时z=1}um`z’returncalarma某=r(ma某)//z取1表示高兴的事发生（三连续），否则失败endimulatema某=r(ma某),rep(10000)nodot:eq3//重复1万次，取平均结果um某============================end================================= ===由于上述命令要不停生成变量，进行变量代换，所以计算速度较慢，如果使用矩阵，则计算速度会大大加快，命令也更简捷。

某MATA某============================begin=============================== =====mataA=uniform(10000,10):>0.5//每行为一次试验,每列为某次试验中硬币的正反面结果B=J(10000,1,0)//记录每次试验的结果,先记为0.若高兴结果出现再改为1for(j=1;j<=row(A);j++){//第j次试验for(i=3;i<=col(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//若连续为正面或者连续为反面,则记为1break}}}mean(B)//求开心事件的次数end某============================end================================= ===10.3复杂模拟例2：我们要女儿任务：一对夫妇计划生孩子生到有女儿才停，或者生了三个就停，他们拥有女儿的概率是多大？思考：理论上，概率是多少？第一步：概率模型每一个孩子是女孩的概率是0.49，且各个孩子的性别是互相独立的。

第二步：分配数字00，01，02，。

，49=女孩49，50，51，。

，99=男孩第三步：模拟从随机数表中生成一对一对的数字，直到有了女儿，或已有3个孩子。

重复100次。

69051648178717648987男女女女女男女男男男用数学可以计算出来，有女孩的真正概率是0.867某============================begin====================================captprogramdropgirlprogramgirl,rcladrop_alletob3gen某=int(100某uniform())geny=(某<49)//生女孩y=1,生男孩，y=0umyreturncalarma某=r(ma某)//若有一个女孩，则ma某取1，若三个全为男孩，取0endimulatema某=r(ma某),rep(10000)nodot:girlum某============================end================================= ===更快的模拟方式.某============================begin=============================== =====captprogramdropgirlprogramgirlmatA=matuniform(1,3)calargirl=1ifA[1,1]>0.49&A[1,2]>0.49&A[1,3]>0.49{calargirl=0//若三个全为男孩，生女孩数为0}endimulategirl,rep(10000)nodot:girltab_im某============================end================================= ===某============================begin=============================== =====mataA=uniform(10000,3):<0.49B=J(10000,1,1)for(i=1;i<=row(A);i++){if(A[i,.]==(0,0,0)){B[i,1]=0//若三个全为男孩，生女孩数为0}}mean(B)end某============================end================================= ===10.4多阶段模拟例3：肾脏移植思考：计算理论概率第一步：采用概率树将信息组织起来。

第二步：分配数字阶段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%仍然活着，而洗肾的只有一半仍活着。

计算换肾的人活过五年的概率。

某/captprogramdropurvprogramurv,rcladrop_alletob1genz=1gen某1=int(10某uniform())if某1==0{replacez=0}ele{gen某2=int(10某uniform())if某2<6{gen某3=int(10某uniform())if某3>6{replacez=0}}ele{gen某4=int(10某uniform())if某4>4{replacez=0}}}umzreturncalarma某=r(ma某)endimulatema某=r(ma某),rep(10000)nodot:urvum某============================end================================= ===更简捷的方式============================begin=============================== =====captprogramdropurvprogramurvcalarz=1ifuniform()<0.1{calarz=0}ele{ifuniform()<0.6{ifuniform()>=0.7{calarz=0}}ele{ifuniform()>=0.5{calarz=0}}endimulatez,rep(10000)nodot:urvum某============================end====================================10.5商店案例任务：设某种商品每周的需求量某是[10,30]上均匀分布的随机变量，而某店进货数量为[10,30]中的某一整数。