R语言蒙特卡洛计算和快速傅立叶变换计算矩生成函数

发布于

原文链接：http://tecdat.cn/?p=13734

概率论中，矩生成函数（Moment-generating Function）和特征函数（Characteristic Function）是定义 概率分布函数的另一种形式。

特征函数能够唯一确定随机变量的概率分布，如果随机变量的概率密度函数f(x)存在，特征函数相当于 f(x)的傅里叶变换。

如果随机变量分布的矩母函数存在，那么矩母函数和特征函数之间存在关系。

  • 蒙特卡洛计算

可以使用蒙特卡洛模拟来计算矩生成函数函数，

代码语言：javascript
复制

> F=function(x) ifelse(x<0,0,1-exp(-x)/3)
> Finv=function(u) uniroot(function(x) F(x)-u,c(-1e-9,1e4))$root

或

代码语言：javascript
复制


> Finv=function(u) ifelse(3*u>1,0,uniroot(function(x)


    

  • F(x)-u,c(-1e-9,1e4))$root))
代码语言：javascript
复制



> plot(u,v,type="b",col='blue')

> lines(u,Mtheo(u),col="red")

可以计算

代码语言：javascript
复制


> M(3)

[1] 5748134

有限总和始终可以通过数字计算。就算在这里 不存在。就像Cauhy样本的平均值一样，即使期望值不存在，我也总是可以计算出来

代码语言：javascript
复制


> mean(rcauchy(1000000))

[1] 0.006069028

这些生成函数在存在时会很有趣。也许使用特征函数是一个更好的主意。

  • 特征函数

当我们处理独立随机变量的总和时，特征函数很有趣，因为总和的特征函数是特征函数的乘积。考虑计算Gamma随机变量复合和的99.5％分位数的问题，即

策略是分散损失金额，

然后，要计算的代码 ， 我们用

99.5％分位数

代码语言：javascript
复制
代码语言：javascript
复制
> sum(cumsum(f)<.995)
代码语言：javascript
复制
考虑以下损失金额
代码语言：javascript
复制


> print(X[1:5])

[1] 75.51818 118.16428 14.57067 13.97953 43.60686
代码语言：javascript
复制
让我们拟合一个伽玛分布。我们可以用
代码语言：javascript
复制


shape         rate

1.309020256   0.013090411

(0.117430137) (0.001419982)


> alpha

[1] 1.308995

> beta

[1] 0.01309016
代码语言：javascript
复制
无论如何，我们都有个人损失的Gamma分布参数。并假设泊松计数变量的均值为
代码语言：javascript
复制
代码语言：javascript
复制
> lambda <- 100
代码语言：javascript
复制
同样，可以使用蒙特卡洛模拟。我们可以使用以下通用代码：首先，我们需要函数来生成两种感兴趣的变量，

如果我们生成一百万个变量，我们可以得到分位数的估算，

代码语言：javascript
复制


> set.seed(1)

> quantile(rcpd4(1e6),.995)

99.5%

13651.64

另一个想法是记住Gamma分布的比例：独立Gamma分布的总和仍然是Gamma（在参数上有附加假设，但在此我们考虑相同的Gamma分布）。因此，可以计算复合和的累积分布函数，

如果我们求解那个函数，我们得到分位数

代码语言：javascript
复制

> uniroot()$root
[1] 13654.43

这与我们的蒙特卡洛计算一致。现在，我们也可以在此处使用快速傅立叶变换，

代码语言：javascript
复制


> sum(cumsum(f)<.995)

[1] 13654

让我们比较获得这三个输出的计算时间

代码语言：javascript
复制


> system.time

user      system     elapsed

2.453       0.106       2.611

> system.time

user      system     elapsed

0.041       0.012       0.361

> system.time

user      system     elapsed

0.527       0.020       0.560

#function