统计模拟和R第讲.ppt

f=function(n){ X=0;Z=0 for(i in 1:n){ U1=runif(1); Y=-log(U1);U2=runif(1) while(U2exp(-(Y-1)^2/2)){ U1=runif(1);U2=runif(1); Y=-log(U1)} X[i]=Y U=runif(1) if(U=1/2)Z[i]=X[i] else Z[i]=-X[i] } Z } 注:若想生成 的随机变量,只需取 正态随机变量的极坐标法 设X,Y是两个独立的标准正态随机变量,记 为(X,Y) 的极坐标. 由X和Y独立,它们的联合密度为 由极坐标变换可计算出 的联合密度 正态随机变量的极坐标法 易看出上述密度为均值为2的指数密度与 上的均匀 密度的乘积. 产生两个独立的标准正态随机变量算法: Step1:生成随机数 Step2:令 Step3:令 此式称为Box-Muller变换. 正态随机变量的极坐标法 由于使用上述算法效率不高,因为需要计算正余弦函数 改进算法如下: 产生两个独立的标准正态随机变量算法: Step1:生成随机数 Step2:令 Step3:如果S1,则转至Step1 Step4:产生独立的标准正态随机变量 f=function(n){ X=0;Y=0 for(i in 1:n){ repeat{ U1=runif(1); U2=runif(1) V1=2*U1-1;V2=2*U2-1;S=V1^2+V2^2 if(S=1)break } X[i]=sqrt(-2*log(S)/S)*V1 Y[i]=sqrt(-2*log(S)/S)*V2 } data.frame(X=X,Y=Y) } 泊松过程的生成 希望产生一个参数为 的泊松过程的前n个事件 发生的时间,利用如下的结论: 该过程的两个相继事件发生的时间间隔为参数为 的指数随机变量. 先生成n个随机数 ,令 表示此过程第i-1个和第i个事件发生的间隔时间,第j个 事件发生的时刻等于前j个间隔时间之和,则前n个事件 发生的时间分别为 泊松过程的生成 想模拟此泊松过程在时刻T前的状态,由上述分析,相继 产生其间隔时间,直至其和超过T. 一些记号:T:模拟的总时间, t:表示时间,I表示时间t内发生 的事件数,S(I):最近事件发生的时间. Step1:t=0,I=0 Step2:生成一个随机数 Step3: ,如果tT,停止. Step4:I=I+1,S(I)=t 一个参数为 的泊松过程在时刻T前的状态模拟 Step5:转至Step2 f=function(T,lamb){ t=0;I=0;S=0 repeat{ U=runif(1);t=t-1/lamb*log(U) if(tT)break I=I+1;S[I]=t } list(S=S,I=I) } 另外一种模拟参数为 的泊松过程时刻T前状态的方法: 先模拟时刻T发生的事件总数N(T),由于N(T)是一均值 的泊松随机变量.若N(T)的模拟值为n,则生成n个随机数 取 为时刻T时此泊松过程的 事件发生的时间集合. 其对应的程序为: f1=function(T,lamb){ N=0; S=0 N=rpois(1,lamb*T) U=runif(N) S=sort(T*U) list(S=S,I=N) } 比较其效率 system.time(f(10000,5)) 用户 系统 流逝 14.92 0.00 15.38 system.time(f1(10000,5)) 用户 系统 流逝 0.03 0.00 0.03 例:设公共汽车到达某比赛场地的时间服从参数5(即每 小时有5个事件发生)的泊松过程,每辆公共汽车等可能 地有20,21,22,…,40个球迷乘坐.且不同公共汽车上的 球迷数是相互独立的.编写模拟t=1时刻抵达此比赛场地 的球迷数的程序.并把模拟的球迷数与理论球迷数作比较 X=20:40 sum(sample(X,f1(1,5)$I,replace=TRUE)) X=20:40 for(i in 1:10000){ Y[i]=sum(sample(X,f1(1,5)$I,replace=TRUE)) } mean(Y) 分析：首先要模拟出在此段时间内到达的公共汽车数，再模拟 每辆车上的球迷数，则球迷数就是它们的和 R中