202211011039+梁子俊+第4周作业

题5.3

(反函数法)

\[ U=\frac{x^2+x}{2} \] \[ 2U+\frac{1}{4}=(x+\frac{1}{2})^2 \] \[ x=\sqrt{2U+\frac{1}{4}}-\frac{1}{2} \]

U=runif(1)
X=sqrt(2*U+1/4)-1/2
X

## [1] 0.129529

题5.8

复合抽样法

#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.6616721

#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.634159

#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.9004155

题5.10

prob=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.1061

题5.14

F是已知a<=X<=b时X的条件分布因为 \[ F(x|a\leq X\leq b)=\frac{P(X\leq x,a\leq X\leq b)}{P(a\leq X\leq b)}=\frac{G(x)-G(a)}{G(b)-G(a)},a\leq x\leq b \]
按照舍选抽样法生成的Y，其分布函数为 \[ F(y)=P(Y\leq y)=P(X\leq y|a\leq X\leq b)=\frac{G(x)-G(a)}{G(b)-G(a)},a\leq y\leq b \]

题5.20

舍选抽样法

标准正态分布的绝对值的密度函数与均值为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] -1.174913

题5.22

#由指数分布产生泊松分布
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] 9987

S[1:100]

##   [1]   0.06539254   1.55802734   1.87167353   2.56573573   3.50584603
##   [6]   3.54888404   4.01710213   4.03636324   5.22353291   5.91070440
##  [11]   9.16854512  10.99332152  11.92933415  12.08837301  14.94470060
##  [16]  15.54247935  16.68937114  18.02198784  18.97397911  19.87514087
##  [21]  20.61851680  21.78530211  22.13568685  22.41193657  22.44860410
##  [26]  23.36118626  24.55216116  28.20254576  29.76581327  30.77601713
##  [31]  31.97215715  32.70828804  33.66005302  37.01387314  37.40033561
##  [36]  37.49097240  38.07186469  39.98189591  40.74800705  41.61959386
##  [41]  43.01218305  43.72935796  45.13666559  45.57243396  48.13876094
##  [46]  48.55406048  51.09972304  51.16977551  51.97633253  52.06484284
##  [51]  52.49725704  52.93231874  53.88095697  55.57361510  57.01787488
##  [56]  57.84002745  60.88704858  62.41217580  63.77778398  65.55262985
##  [61]  66.02115973  68.41753555  68.78915542  71.51100276  79.36057032
##  [66]  79.62507685  80.92719743  81.62023525  81.88174182  83.46170796
##  [71]  83.76536352  84.24359852  86.24536388  89.02567672  89.93168986
##  [76]  90.68663481  90.96414303  92.94582491  93.81328478  94.24363524
##  [81]  95.13902904  96.88124518  97.19802653  98.94605511  99.45853250
##  [86]  99.90700443 100.55841097 100.96203504 101.45389939 102.15322310
##  [91] 102.49681381 105.69431689 106.49993828 106.85305904 111.67878476
##  [96] 112.57641382 113.08560420 114.58616388 114.65536700 116.07170459

题5.25

#瘦身方法
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.5207

S #某一次泊松过程的事件发生时刻

## [1] 0

end_time=Sys.time() #记录程序的结束时刻
end_time-start_time #程序的运行时间

## Time difference of 1.703132 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.5971

S #某一次泊松过程的事件发生时刻

## [1] 0

end_time=Sys.time() #记录程序的结束时刻
end_time-start_time #程序的运行时间

## Time difference of 1.286681 secs

#复合方法效率略优于瘦身方法

题5.26

#复合方法结合分段瘦身方法
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.0988

S #某一次泊松过程的事件发生时刻

## [1] 0

end_time=Sys.time()
end_time-start_time

## Time difference of 1.516857 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] 70.0748

S #某一次泊松过程的事件发生时刻

## [1] 0

end_time=Sys.time() #记录程序的结束时刻
end_time-start_time #程序的运行时间

## Time difference of 0.903168 secs

#分布函数法略优