R语言泊松过程及在随机模拟应用可视化

简介: R语言泊松过程及在随机模拟应用可视化

泊松分布是概率论中最重要的分布之一,在历史上泊松分布是由法国数学家泊松引人的。近数十年来,泊松分布日益显现了其重要性而将泊松随机变量的概念加以推广就得到了泊松过程的概念点击文末“阅读原文”获取完整代码数据


相关视频

image.png 泊松过程是被研究得最早和最简单的一类点过程,它在点过程的理论和应用中占有重要的地位。泊松过程在现实生活的许多应用中是一个相当适合的模型,它在物理学、天文学、生物学、医学、通讯技术、交通运输和管理科学等领域都有成功运用的例子。

考虑一个具有非均匀强度的泊松过程。

假设顾客到达服务站的人数服从强度为4的泊松过程,到达的顾客很快就可以接受服务,并且假设服务时间是独立的并且服从一个普通的分布,记为G。

为了计算在时刻t已完成服务和正在接受服务的顾客的联合分布,把在时刻t<=3完成服务的顾客称为第一类,在时刻t未完成服务的顾客称为第二类顾客,现在,如果第一个顾客到来的时间为tSS,,如果他的服务时间少于t - s,那么 他就是第一类顾客,并且因为服务时间服从G分布,所以服务时间少于t - s的概率为G(t - s)因而,P(s) = G(t -s); S≤ t。利用定理2我们得到的)(1tN的分布。到时间t为止,已完成服务的顾客的数目服从泊松分布 在这里,我们考虑一个确定性函数,而不是随机强度。定义累积强度:

9919267fafb873761958fcefe1e26524.png

发生的事件的数量是随参数分布的泊松的随机变量。

lambda=function(x) 100*(sin(x*pi)+1)

这个想法是在有限的时间间隔上生成泊松过程

1.            开始 7d19e9c0ac3f73c59e29929a1f5cb7b2.png

2.            生成 7a069a52cef2d201e6d866b15d243cec.png

3.            设置 a58e49ba5e18e2db4a93e893e5bffd6f.png

4.            设置 6da379d05d8b45e3eae404e330b74c28.png 属于 e6559de3334c37ce36ba44c91ce4d9d0.png

5.            更新t

6.            返回第 2步.

为了得到最小值 9d03e825f5dc678614700027536d48f3.png ,考虑代码

这里,生成泊松过程的代码是

X= 0  
   while(X[length(X)]<=Tmax){  
     u=runif(1)

在这里,我们得到以下直方图,

hist(X,breaks=seq(0,max(X)+1,by=.1),col="yellow")  
   u=seq(0,max(X),by=.02)

9a2ed065fd917b69eb13e3b1aa47ff52.png

点击标题查阅往期内容


R语言和Python用泊松过程扩展:霍克斯过程Hawkes Processes分析比特币交易数据订单到达自激过程时间序列



左右滑动查看更多



01


现在考虑另一个策略。这个想法是在下一个事件之前使用条件分布,假设一个事件发生在时间t,

df2aef5d7415697e037c95408d13e133.png

1.            开始 e5ecd4cefe796bf5579b0aa7fa727e6e.png


2.            生成 a472a3e78f892ed149ebdea6f286b4ef.png

3.            设置 2d0ea481570d4db52756f4c4f52bca6b.png

4.            更新t

5.            返回第 2步.

我们可以使用二分法算法,

for(j in 1:20){  
       if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}  
       if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
       
       
       a=0  
      b=Tmax  
      for(j in 1:20){  
        if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}  
        if(Ft((a+b)/2)>=u)

在这里,我们得到以下直方图,

lines(u,lambda(u)/10,lwd=2,col="red")

678a22e32938c24220dbf45c48133cd3.png

第三个代码是基于经典算法在有限间隔上生成均匀泊松过程:首先,我们生成事件数,然后,绘制均匀变量,然后对它们进行排序。

1.            生成时间间隔

4b0f0659a6c952baba1f70d8c175260e.png

2.        生成 bba664c685542e2bbf145f33cdb408b7.png 其中   bf8d92fec319fdebf671abd7c3589b32.png  

3.            设置 6ec99049554b88f1986232525216b6ca.png i.e. 8362b40edb3aafcf72b12ac0ab56dacc.png

4.            更新 dd820338b70345bce6f71d424152956e.png ‘s

这个算法非常简单,而且速度也很快。这是一个反函数的函数,它不在循环中,

n=rpois(1,Lambda(Tmax))  
   Ft=function(x) Lambda(x)/Lambda(Tmax)  
   Ftinv=function(u){  
     a=0  
     b=Tmax  
     for(j in 1:20){

在这里,我们得到以下直方图

u=seq(0,max(X),by=.02)

9154c14dcf1a77d0ebb81c5012a28175.png

一种替代方案基于拒绝技术 。这里,我们需要一个强度的上限,这样计算可能会快得多。

1.            开始 4a2cedbd0cd24cf0ff8edb00c428e7e2.png

2.            生成 7a069a52cef2d201e6d866b15d243cec.png

3.            设置 00d8abdd50dc58db2f8e6ba47caef729.png

4.            生成 cd118b018f89146149b482657b31f907.png

5.            如果 c7c5d30faf8ee66850afa3723aa1c0a9.png  然后 更新 1c0fe6679196a083b733eefeeeed15da.png

6.            返回第 2步.

这里,考虑一个恒定的上界,

t=0  
   X=  0  
   while(X[length(X)]<=Tmax){  
     u=runif(1)

在这里,我们得到以下直方图

hist(X,breaks=seq(0,max(X)+1,by=.1),col="yellow")  
   u=seq(0,max(X),by=.02)

f1c9dd93a331253a1a61ce34c0310d07.png

最后,一个也是基于拒绝技术,与第二个混合。也就是说 定义

91a169669ce0c1309b871e27b3ebe18c.png

这个函数可以很容易

9271fc3eae0d63e543b0abd691adbe61.png

1.            开始

4a2cedbd0cd24cf0ff8edb00c428e7e2.png

2.            生成

7b4609fef93d96e227cbb868422ac773.png

3.            设置

041e7492cf0e0d2929cfdced20cb73c2.png

4.            生成

7a069a52cef2d201e6d866b15d243cec.png

5.            如果 0d2452fa43cb2e4f995fae7b8613fb65.png 然后 更新

6.            返回第二步.

Ftinvu=function(u) -log(1-x)/lambdau  
     x=Ftinvu(runif(1))

在这里,我们得到以下直方图

f1c9dd93a331253a1a61ce34c0310d07.png

相关文章
|
4天前
【R语言实战】——fGARCH包在金融时序上的模拟应用
【R语言实战】——fGARCH包在金融时序上的模拟应用
|
4天前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
4天前
|
Web App开发 数据可视化 数据挖掘
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
利用R语言进行聚类分析实战(数据+代码+可视化+详细分析)
|
4天前
|
数据采集 数据可视化
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
利用R语言进行因子分析实战(数据+代码+可视化+详细分析)
|
4天前
|
移动开发 数据可视化
广义线性模型beta二项分布的淋巴结疾病风险预测可视化R语言2实例合集|附数据代码
广义线性模型beta二项分布的淋巴结疾病风险预测可视化R语言2实例合集|附数据代码
|
4天前
|
机器学习/深度学习 数据可视化
R语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
R语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
|
4天前
|
数据可视化
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码2
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码
|
4天前
|
数据可视化 数据挖掘
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码1
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码
|
4天前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
4天前
利用R语言进行典型相关分析实战
利用R语言进行典型相关分析实战

热门文章

最新文章