\[ U=\frac{x^2+x}{2} \] \[ 2U+\frac{1}{4}=(x+\frac{1}{2})^2 \] \[ x=\sqrt{2U+\frac{1}{4}}-\frac{1}{2} \]
## [1] 0.8211206#F_1(x)=x,0<x<1;p_i=1/3
#F_2(x)=x^3,0<x<1;p_i=1/3
#F_3(x)=x^5,0<x<1;p_i=1/3
U=runif(3)
U0=runif(1)
if(U0<=1/3){
X=U[1]
} else if(U0<=2/3){
X=U[2]^(1/3)
} else{
X=U[3]^(1/5)
}
X## [1] 0.5109733#F_1(x)=1-exp(-2x),0<x;p_i=1/3
#F_2(x)=F_3(x)=x,0<x<1;p_i=1/3
U=runif(3)
U0=runif(1)
if(U0<=1/3){
X=-log(1-U[1])/2
} else if(U0<=2/3){
X=U[2]
} else{
X=U[3]
}
X## [1] 0.8370865#F_i(x)=x^i,0<x<1;p_i=α_i
n=10;alpha=runif(10);alpha=alpha/sum(alpha)
alpha=cumsum(alpha)
U=runif(n)
U0=runif(1)
X=0;i=0
while(X==0){
i=i+1
if(U0<=alpha[i]){
X=U[i]^(1/i)
}
}
X## [1] 0.6725448prob=0 #初始化
for(k in 1:1e4){#模拟1e4次
U=runif(1000) #模拟1000次二项分布
n=sum(U<=0.05) #成功概率为0.05，n为成功次数
U=runif(n)
t=-800*log(prod(U)) #n个指数分布随机变量的和
V=rep(0,n+1)
V[n+1]=1
V[2:n]=sort(runif(n-1)) #n-1个均匀分布的次序统计量
X=rep(0,n)
for(i in 1:n){
X[i]=t*(V[i+1]-V[i])
#前i个指数分布随机变量的和的条件分布与n-1个均匀分布的次序统计量同分布
}
prob=prob+(sum(X)>50000) #超过50000则+1
}
prob=prob/1e4
prob## [1] 0.1101标准正态分布的绝对值的密度函数与均值为1的指数分布的密度函数的比值为 \[ \frac{f(x)}{g(x)}=\frac{2}{\sqrt{2\pi}}e^{x-\frac{x^2}{2}} \] 最大值为 \[ c=\frac{2}{\sqrt{2\pi}}e^{1-\frac{1^2}{2}}=\sqrt{\frac{2e}{\pi}} \] 此时 \[ \frac{f(x)}{cg(x)}=e^{x-\frac{x^2}{2}-\frac{1}{2}}=e^{-\frac{(x-1)^2}{2}} \] 用舍选抽样法选出X的绝对值后，再以1/2的概率随机取正负号
X=Inf
mu=1;sigma=2
while(X==Inf){#当X有取值时退出循环
U=runif(3)
Y=-log(U[1]) #取指数分布的随机变量
if(U[2]<=exp(-(Y-1)^2)/2){ #取X的绝对值
if(U[3]<=1/2) X=Y #取正负
else X=-Y
}
}
X=mu+X*sigma #均值为μ，标准差为σ
X## [1] 3.506146#由指数分布产生泊松分布
T=1e4;lambda=1
#为了运行程序给各变量设了初值，实际应用中初值可以任选
t=0;i=0;S=rep(0,1)
#t为当前时刻，i为泊松过程的计数，S储存事件发生时刻
while(t<T){
U=runif(1)
t=t-log(U)/lambda #模拟事件发生的时间间隔
if(t>=T) break
i=i+1
S[i]=t
}
i## [1] 9966## [1] 1.411444 3.856369 5.528962 6.124154 6.745523 7.008187 7.878562
## [8] 8.319570 8.798036 8.879112 9.155390 11.211243 12.874869 14.967165
## [15] 15.676646 16.120426 16.562775 16.638961 17.823894 17.911738 18.088964
## [22] 19.587882 20.453731 21.223846 21.908984 23.228827 23.709385 23.964406
## [29] 25.129828 25.842041 26.171998 26.547975 28.631460 30.216562 31.551500
## [36] 32.758725 32.897087 33.459191 34.734870 35.405234 36.034224 38.591840
## [43] 39.198833 40.623053 42.360686 43.900657 43.924671 44.564689 44.967087
## [50] 46.250147 46.950782 50.147997 50.462425 52.186215 54.057015 54.204769
## [57] 54.434045 56.234000 56.247618 56.549786 56.572033 60.226080 60.258429
## [64] 60.336780 61.207674 62.808385 65.048705 65.380524 66.752230 68.506118
## [71] 69.987092 69.989217 70.399563 70.579251 72.078758 72.648985 74.841777
## [78] 76.338289 76.353363 77.686829 78.549210 79.051202 79.533249 80.224072
## [85] 81.757547 82.048195 82.691386 84.282481 84.605806 85.398703 85.875943
## [92] 87.658994 88.795617 90.784486 92.655183 95.184342 97.785427 97.851001
## [99] 98.953108 99.997489#瘦身方法
start_time=Sys.time() #记录程序的开始时刻
i=0;lambda=7
#i为泊松过程的计数
for(k in 1:1e4){
T=10;lambda=7
t=0;S=rep(0,1)
#t为当前时刻，i为泊松过程的计数，S储存事件发生时刻
while(t<T){
U=runif(1)
t=t-log(U)/lambda #模拟事件发生的时间间隔
if(t>=T) break
U=runif(1)
if(U*lambda<3+4/(1+t)){#舍选抽样
i=i+1
if(k==1) S[i]=t
}
}
}
i/1e4 #事件的平均发生次数## [1] 39.5752## [1] 0## Time difference of 1.789313 secs#复合方法
start_time=Sys.time() #记录程序的开始时刻
i=0;lambda1=3+4/11;lambda2=4-4/11
#i为泊松过程的计数，lambda1是齐次泊松过程，lambda2是非齐次泊松过程
for(k in 1:1e4){
T=10
t=0;S=rep(0,1)
#t为当前时刻，S储存事件发生时刻
while(t<T){
U=runif(1)
t=t-log(U)/lambda1 #模拟事件发生的时间间隔
if(t>=T) break
i=i+1
if(k==1) S[i]=t
}
t=0
while(t<T){
U=runif(1)
t=t-log(U)/lambda2 #模拟事件发生的时间间隔
if(t>=T) break
U=runif(1)
if(U*lambda2<4/(1+t)-4/11){#舍选抽样
i=i+1
if(k==1) S[i]=t
}
}
}
i/1e4 #事件的平均发生次数## [1] 39.6564## [1] 0## Time difference of 1.456223 secs#复合方法结合分段瘦身方法
start_time=Sys.time() #记录程序的开始时刻
i=0
#i为泊松过程的计数
for(k in 1:1e4){
t=0;T=5;S=rep(0,1);lambda=1
#t为当前时刻，S储存事件发生时刻
#前5个单位时间，用瘦身方法处理
while(t<T){
U=runif(1)
t=t-log(U)/lambda #模拟事件发生的时间间隔
if(t>=T) break
U=runif(1)
if(U<t/5){#舍选抽样
i=i+1
if(k==1) S[i]=t
}
}
t=5;T=10;J=6;lambda=1
#t为当前时刻，S储存事件发生时刻
#后5个单位时间，用复合方法结合分段瘦身方法处理
#首先模拟齐次泊松过程，对于长度为1的时间区段，lambda固定为1+5*floor(t-5)，即1,6,11,16……
while(t<T){
U=runif(1)
delta_t=-log(U)/lambda #模拟事件发生的时间间隔
while(t+delta_t>J){
#判断下一次事件的发生时间是否超出了当前时间区段
#若超出，则依次平移当前时刻t，直到找到下一次事件发生的时间区段
delta_t=(t+delta_t-J)*lambda/(lambda+5)
#不同时间区段的lambda不同，因此要更新时间间隔
lambda=lambda+5
#更新lambda
t=J
#将t移动到该时间区段的起始位置
J=J+1
#平移当前时间区段
}
t=t+delta_t
if(t>=T) break
i=i+1
if(k==1) S[i]=t
}
t=5;T=10;J=6;lambda=5
#再模拟非齐次泊松过程，对于长度为1的时间区段，lambda为5*(t-floor(t))，其最大值为5
while(t<T){
U=runif(1)
delta_t=-log(U)/lambda #模拟事件发生的时间间隔
while(t+delta_t>J){
#判断下一次事件的发生时间是否超出了当前时间区段
#若超出，则依次平移当前时刻t，直到找到下一次事件发生的时间区段
delta_t=t+delta_t-J
#在不同的时间区段，lambda的最大值总是5，因此无需更新lambda
t=J
J=J+1
}
t=t+delta_t
if(t>=T) break
U=runif(1)
if(U*lambda<(t-floor(t))*5){#舍选抽样
i=i+1
if(k==1) S[i]=t
}
}
i
}
i/1e4 #事件的平均发生次数## [1] 70.0387## [1] 0## Time difference of 1.478707 secs在前5个单位时间内，由于 \[ \int^x_0(\frac{t+y}{5})\,dy=\frac{t}{5}x+\frac{1}{10}x^2 \] 而 \[ 1-e^{-(\frac{t}{5}x+\frac{1}{10}x^2)}=U \] 的反函数为 \[ {\frac{t}{5}x+\frac{1}{10}x^2}=-\log (1-U) \] \[ x^2+2tx=-10\log U(因U与1-U同分布) \] \[ (x+t)^2=-10\log U+t^2 \] \[ x=\sqrt{-10\log U+t^2}-t \] 而在后5个单位时间内，由于 \[ \int^x_0(1+5(t+y-5))\,dy=-24x+5tx+\frac{5}{2}x^2 \] 而 \[ 1-e^{-(-24x+5tx+\frac{5}{2}x^2)}=U \] 的反函数为 \[ -24x+5tx+\frac{5}{2}x^2=\log (1-U) \] \[ (\sqrt{5}x+\frac{5t-24}{\sqrt{5}})^2=-2\log U+\frac{(5t-24)^2}{5}(因U与1-U同分布) \] \[ x=\frac{\sqrt{-2\log U+\frac{(5t-24)^2}{5}}-\frac{5t-24}{\sqrt{5}}}{\sqrt{5}} \]
start_time=Sys.time() #记录程序的开始时刻
i=0
#i为泊松过程的计数
for(k in 1:1e4){
t=0;T=5;S=rep(0,1)
#t为当前时刻，S储存事件发生时刻
while(t<T){
U=runif(1)
t=t+sqrt(-10*log(U)+t^2)-t #模拟事件发生的时间间隔
if(t>T) break
i=i+1
if(k==1) S[i]=t
}
t=5;T=10
while(t<T){
U=runif(1)
t=t+sqrt(-2*log(U)+(5*t-24)^2/5)/sqrt(5)-(5*t-24)/5 #模拟事件发生的事件间隔
if(t>T) break
i=i+1
if(k==1) S[i]=t
}
}
i/1e4 #事件的平均发生次数## [1] 69.9721## [1] 0## Time difference of 1.004304 secs