随机过程 番外篇（随机拟合作业解答）

一晚上写了三道随机过程的随机模拟的代码，分享出来给大家做个参照。

---

1.如果一个随机变量服从的是期望为$\mu$，协方差矩阵为$\Sigma$的正态分布。那么如果这个多维正态分布是非退化的，那么这个协方差矩阵必然是正定的(why)，而且结合其本身的属性，可以知道$\Sigma$是一个对称正定矩阵。对于这样的矩阵，我们可以对其进行Cholesky分解（留“坑“”），将其分解为$\Sigma =L\cdot L^T$的形式，其中$L$是下三角矩阵。

这样的话，如果随机变量$X$服从$N(\text{0},I)$的分布，那么$Y=L\cdot X+\mu$服从$N(\mu ,\Sigma )$的正态分布。按照这个思路，我们可以对$X$进行线性变换得到我们希望得到的随机变量。

我用python实现的代码如下:

from scipy import linalg
import numpy as np
import random
import matplotlib.pyplot as plt
'''
Mu=np.array([[1,1]]).T
print(np.shape(Mu))
Delta=np.array([[1,1],[1,3]])
L = linalg.cholesky(Delta, lower=True)
print(np.dot(L,L.T))
# L \cdot L^T =Delta=E((x-mu) \cdot (x-mu)^T)(2*1\cdot 1*2)=E(L y \cdot y^T L^T)
#设y服从N(0,E),则 x=Ly+Mu
'''
def TwoD_Normalnum(Mu,Delta):
    '''
    Create random number of 2-Dimension Norm distribution on expection of Mu and Coverence matrix of Delta.
    :param Mu: expection
    :param Delta: Coverence matrix
    :return: random vector
    '''
    L = linalg.cholesky(Delta, lower=True)#获得下三角矩阵
    x=random.normalvariate(0,1)
    y=random.normalvariate(0,1)
    Y=np.array([[x],[y]])
    X=np.dot(L,Y)+Mu
    return X

if __name__=='__main__':
    Mu=np.array([[1,1]]).T
    Delta=np.array([[1,1],[1,3]])
    x=[]
    y=[]
    for i in range(2000):
        a=TwoD_Normalnum(Mu,Delta)
        x.append(a[0])
        y.append(a[1])
    plt.scatter(x,y)
    plt.grid(True)
    plt.savefig('2003.jpg')
    plt.show()

得到的以$\left(\begin{array}{c}1\\ 1\end{array}\right)$为期望，以$\left(\begin{array}{cc}1& 1\\ 1& 3\end{array}\right)$为协方差矩阵的二维正态分布的图像如下：

---

2.首先我们确定整体的思路，用计算机模拟大数定理，其实是看频率随n变化过程是否能够收敛到一个点附近。因此我们通过画出$n-Frequent$图像来对大数定律进行验证。

首先我们学过，正态分布符合大数定理，而Cauchy分布因为期望不存在而不满足大数定理。

正态分布在random库中有相应函数random.normalvariate()可以生成随机数，而Cauchy分布则没有现成的函数了。

这时候，我们联想到随机变量的存在性定理，对于给定的分布函数，我们可以通过对服从$U[0,1]$的随机变量$X$生成相应的随机变量$Y$。
$$P(F^{-1}(X)<x_0)=P(X<F(x_0))=F(x_0)$$
因而随机变量$Y=F^{-1}(X)$即为服从分布函数为$F$的随机变量。

我们用python代码对其进行实现：

import matplotlib.pyplot as plt
import random
import numpy as np
import math

if __name__=='__main__':
    plt.subplot(1,2,1)
    #正态分布(n,frequent)图像
    n=list(np.arange(1,1000000,1))
    F1=[]
    sum=0
    for i in n:
        sum+=random.normalvariate(0,1)
        F1.append(sum/i)
    plt.plot(n,F1)
    plt.grid(True)
    plt.subplot(1,2,2)
    #Cauchy分布(n,frequent)图像
    #Cauchy分布:f=(1/\pi)*(1/(1+x^2))
    #需要通过F^(-1)(X),X~U[0,1]构造复合Cauchy分布的随机变量
    #x=tan(\pi*(F-1/2))
    F2= []
    sum = 0
    for i in n:
        sum += math.tan(math.pi*(random.uniform(0, 1)-1/2))
        F2.append(sum / i)
    plt.plot(n, F2)
    plt.grid(True)
    plt.savefig('1608.jpg')
    plt.show()

得到的对应于正态分布与Cauchy分布的次数-频率图如下图所示：

---

3.第三题用计算机模拟，最容易想到的就是蒙特卡洛方法(MC),通过生成服从$U[0,1]\times U[0,1]$的二维随机变量，检验落点的位置：

如果落点在曲线下方，则记数变量$M$加一，否则不加。
那么根据大数定律，应该有$li{m}_{N\to \infty }\frac{M}{N}=\frac{Int(e^{-x})}{S_{[0,1]\times [0,1]}}$
我们将这个想法代码化，并将结果附后：

import random
import numpy as np
import math
import matplotlib.pyplot as plt

#Using MC mathod to culculate the internal of e^(-x) from 0 to 1.
n=list(np.arange(1,100000,1))
I=[]
sum=0
for i in n:
    X=random.uniform(0,1)
    Y=random.uniform(0,1)
    if Y<math.exp(-X):
        sum+=1
    I.append(sum/i)

plt.plot(n,I)
plt.text(70000,0.7,'stimulink result:'+str(round(I[-1],4)))
plt.text(70000,0.65,'culculate result:'+str(round(1-1/math.e,4)))
plt.savefig('1622.jpg')
plt.show()