STATA入门10 随机模拟.docx

上传人:b****5 文档编号:12596148 上传时间:2023-04-20 格式:DOCX 页数:22 大小:74KB
下载 相关 举报
STATA入门10 随机模拟.docx_第1页
第1页 / 共22页
STATA入门10 随机模拟.docx_第2页
第2页 / 共22页
STATA入门10 随机模拟.docx_第3页
第3页 / 共22页
STATA入门10 随机模拟.docx_第4页
第4页 / 共22页
STATA入门10 随机模拟.docx_第5页
第5页 / 共22页
点击查看更多>>
下载资源
资源描述

STATA入门10 随机模拟.docx

《STATA入门10 随机模拟.docx》由会员分享,可在线阅读,更多相关《STATA入门10 随机模拟.docx(22页珍藏版)》请在冰豆网上搜索。

STATA入门10 随机模拟.docx

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)

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 成人教育 > 专升本

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1