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

对精算科学来说，当我们处理独立随机变量的总和时，特征函数很有趣，因为总和的特征函数是特征函数的乘积。

除了表示中心、离散程序、斜度这些特性外，更高阶的矩可以描述分布的其它特性。矩统计中有重要的地位，比如参数估计的一种重要方法就是利用了矩。然而，根据矩的定义，我们需要对不同阶的X幂求期望，这个过程包含复杂的积分过程，并不容易。矩同样催生了矩生成函数(moment generating function)，它是求解矩的一样有力武器。

可下载资源

完整代码、数据和文档（word）

在概率论中，让 $F(x)=1-e^{-x}/3$ 对于 $x\geq 0$ 和 $F(x)=0$ 对于 $x<0$ 是一些随机变量的累积分布函数 $X$ ，即 $F(x)=\mathbb{P}(X\leq x)$ 。什么是矩生成函数 $X$ ，即 $M(t)=\mathbb{E}(e^{tX})$ ？

如何编写 $\mathbb{E}(g(X))$ ？

在概率教科书中，标准答案是

如果 $X$ 是离散的

$\mathbb{E}(g(X))=\sum_x g(x)\cdot \mathbb{P}(X=x)$

如果 $X$ （绝对）连续，

$\mathbb{E}(g(X))=\int g(x)\cdot f(x)dx$

$f(\cdot)$ 是的密度 $X$ 。这里， $X$ 显然不是离散变量。但是是连续的。需要绘制该分布函数以查看， $x\mapsto F(X)$ ，对所有 $x\in\mathbb{R}$

视频

马尔可夫链蒙特卡罗方法MCMC原理与R语言实现

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

我们有一个不连续的0。因此，我们在这里必须谨慎一些： $X$ 既不是连续的也不是离散的。让我们使用公式，

$\mathbb{E}(Y)=\mathbb{E}(\mathbb{E}(Y\vert Z))$

如果也可以写 $Z\in\{A,B\}$ ，

$\mathbb{E}(Y)=\mathbb{P}(Z=A)\cdot \mathbb{E}(Y\vert Z=A)+\mathbb{P}(Z=B)\cdot \mathbb{E}(Y\vert Z=B)$

这只是说总体平均值是每个子组平均值的重心。 $Y=g(X)$ 然后让 $A=\{X=0\}$ 而 $B=\{X>0\}$ $\mathbb{P}(X\in \{A\cup B\})=1$ ）。

$\mathbb{E}(g(X))=\mathbb{P}(X=0)\cdot \mathbb{E}(g(X)\vert X=0)+\mathbb{P}(X>0)\cdot \mathbb{E}(g(X)\vert X>0)$

让我们考虑三个不同的组成部分。

$\mathbb{P}(X=0)=F(0)=1-1/3=2/3$

$latex i\hbar\frac{\partial}{\partial }\left|\Psi(t)\right>=H\left|\Psi(t)\right>$

$\mathbb{P}(X>0)=1-\mathbb{P}(X=0)=1/3$

$\mathbb{E}(g(X)\vert X=0)=g(0)$

（因为它是一个实值常量），在这里 $g(0)=e^{0}=1$ 。

所以最后，我们计算 $\mathbb{E}(g(X)\vert X>0)$ 。观察一下 $X$ 给定 $X>0$ 是具有密度的（绝对）连续随机变量。观察所有 $x>0$ ，

$\overline{F}_\star(x)=\mathbb{P}(X>x\vert X>0)=\frac{\mathbb{P}(X>x)}{\mathbb{P}( X>0)}=\frac{e^{-x}/3}{e^{-0}/3}=e^{-x}$

和 $f_\star(x)=e^{-x}$ ，即 $X$ 给定 $X>0$ 是指数分布。

因此， $X$ 是指数变量和Dirac质量之间的混合 $0$ 。这实际上是问题的棘手部分，因为当我们看到上面的公式时，它并不明显。

从现在开始，这是高中阶段的计算，

$\mathbb{E}(g(X)\vert X>0)=\int_0^\infty g(x) f_\star (x)dx=\int_0^\infty e^{(t-1)\ x} dx=\frac{1}{1-t}$

如果 $t\leq 1$ 。如果把所有的放在一起

$M(t)=\frac{2}{3}\cdot 1 + \frac{1}{3}\cdot \frac{1}{1-t}=\frac{3-2t}{3-3t}$

蒙特卡洛计算

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

> 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

或（以避免不连续的问题）

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

在这里，逆很容易获得，因此我们可以使用

最受欢迎的见解

1.R语言中实现马尔可夫链蒙特卡罗MCMC模型

2.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

3.R语言中的Stan概率编程MCMC采样的贝叶斯模型

4.R语言Rstan概率编程规划MCMC采样的贝叶斯模型简介

5.matlab实现MCMC的马尔可夫切换ARMA – GARCH模型估计

6.matlab对MCMC贝叶斯方法用于加筋复合板的冲击载荷识别

7.R语言用Backfitting MCMC抽样算法进行贝叶斯推理案例

8.使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析

9.R语言使用马尔可夫链Markov Chain, MC来模拟抵押违约

然后，我们使用


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

蒙特卡洛模拟的问题在于，仅当它们有效时才应使用它们。我可以计算


> M(3)
[1] 5748134

有限总和始终可以通过数字计算。就算在这里 $https://latex.codecogs.com/gif.latex?\mathbb{E}(e^{3X})$ 不存在。就像Cauhy样本的平均值一样，即使期望值不存在，我也总是可以计算出来

> mean(rcauchy(1000000))
[1] 0.006069028

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

生成函数

首先，让我们定义那些函数。

$\left\{\begin{array}{l}\text{if }x>x_0, \ \vert x\vert^k \leq \exp(\vert t x\vert)\\\text{if }x\leq x_0, \ \vert x\vert^k \leq K \cdot \exp(\vert t x\vert)\end{array}\right.$

$\mathbb{E}\left\vert X\vert^k\right) \leq \mathbb{E}\left( e^{\vert tX \vert}\right)\leq \mathbb{E}\left( e^{- tX}\right)+\mathbb{E}\left( e^{ tX}\right)\infty$

如果 $https://latex.codecogs.com/gif.latex?t$ 足够小。

现在，如果我们使用泰勒展开式

$M_X(t)=\mathbb{E}\left( e^{ tX}\right)=\mathbb{E}\left(\sum_{k=0}^\infty \frac{(tX)^k}{k!}\right)=\sum_{k=0}^\infty\frac{t^k}{k!}\mathbb{E}[X^k]$

和

$\frac{\partial^k M_X(t)}{\partial t^k} =\mathbb{E}\left( X^k e^{tX} \right)$

如果我们看一下该函数在0点的导数的值，那么

$\left. \frac{\partial^k M_X(t)}{\partial t^k}\right\vert_{0} =\mathbb{E}\left( X^k \right)$

可以为某些随机矢量在更高维度上定义一个矩生成函数 $https://latex.codecogs.com/gif.latex?\boldsymbol{X}=(X_1,\cdots,X_d)$ ，

$\left. \frac{\partial^2 M_X(\boldsymbol t)}{\partial t_i \partial t_j}\right\vert_{\boldsymbol{0}} =\mathbb{E}\left( X_iX_j \right)$ 如果要导出给定分布的矩，则一些矩生成函数很有趣。另一个有趣的特征是，在某些情况下，此矩生成函数（在某些条件下）完全表征了随机变量的分布。 $https://latex.codecogs.com/gif.latex?%20h%3E0$ ，
$https://latex.codecogs.com/gif.latex?%20M_X(t)=M_Y(t)$ 对所有人 $https://latex.codecogs.com/gif.latex?%20t\in(-h,+h)$ ，然后 $https://latex.codecogs.com/gif.latex?X\overset{\mathcal{L}}{=}Y$ 。

快速傅立叶变换

回想一下欧拉公式，

$e^{it}=\cos(t)+i \ \sin(t)$

因此，看到傅立叶变换就不会感到惊讶。从这个公式，我们可以写

$\phi_X(t)=\mathbb{E}\left( e^{i tX} \right) =\mathbb{E}\left( \cos[tX] \right)+i\ \mathbb{E}\left( \sin[tX] \right)$

使用傅立叶分析中的一些结果，我们可以证明概率函数满足

$\mathbb{P}(X=x)=\lim_{T\rightarrow \infty} \frac{1}{2T}\int_{-T}^{+T}e^{-itx}\phi_X(t)dt$

也可以写成

$f_X(x)=\frac{1}{2\pi}\int_{-\pi}^{+\pi}e^{-itx}\phi_X(t)dt$

如果在点处的分布是绝对连续的，则可以获得类似的关系 $x$ ，

$f_X(x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-itx}\phi_X(t)dt$

实际上，我们可以证明，

$f_X(x)=\frac{1}{2\pi}\int_{-\pi}^{+\pi} \text{Re}\left(e^{-itx}\phi_X(t)\right)dt$

然后可以使用1951年获得的吉尔-佩莱阿兹（Gil-Peleaz）的反演公式来获得累积分布函数，

$F_X(x)=\frac{1}{2} +\frac{1}{2\pi}\int_{0}^{\infty} \frac{e^{itx}\phi_X(-t)+e^{-itx}\phi_X(t)}{it}dt$

这意味着，在金融市场上工作的任何人都知道用于定价期权的公式（例如，参见 Carr＆Madan（1999））。好处是，可以使用任何数学或统计软件来计算这些公式。

特征函数和精算科学

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

$https://latex.codecogs.com/gif.latex?%20S=\sum_{n=1}^N%20X_i$

$https://latex.codecogs.com/gif.latex?%20X_i\sim\mathcal{G}(\alpha,\beta)$ 和 $https://latex.codecogs.com/gif.latex?%20N\sim\mathcal{P}(\lambda)$ 。策略是分散损失金额，

然后，要计算的代码 $https://latex.codecogs.com/gif.latex?%20\tilde%20f(s)=\mathbb{P}(S\in[s\pm1/2])$ ，我们用

99.5％分位数

> sum(cumsum(f)<.995)

考虑以下损失金额


> print(X[1:5])
[1] 75.51818 118.16428 14.57067 13.97953 43.60686

让我们拟合一个伽玛分布。我们可以用

      shape         rate    
  1.309020256   0.013090411 
 (0.117430137) (0.001419982)


> alpha
[1] 1.308995
> beta
[1] 0.01309016

无论如何，我们都有个人损失的Gamma分布参数。并假设泊松计数变量的均值为

> lambda <- 100

同样，可以使用蒙特卡洛模拟。我们可以使用以下通用代码：首先，我们需要函数来生成两种感兴趣的变量，

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

> set.seed(1)
> quantile(rcpd4(1e6),.995)
   99.5% 
13651.64

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

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

> uniroot()$root
[1] 13654.43

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

> sum(cumsum(f)<.995)
[1] 13654

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

> 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

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！