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

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

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

这样的话，如果随机变量 X X X服从 N ( 0 , I ) N(\textbf{0},I) N(0,I)的分布，那么 Y = L ⋅ X + μ Y=L\cdot X+\mu Y=L⋅X+μ服从 N ( μ , Σ ) N(\mu,\Sigma) N(μ,Σ)的正态分布。按照这个思路，我们可以对 X X 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()

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

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

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

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

这时候，我们联想到随机变量的存在性定理，对于给定的分布函数，我们可以通过对服从 U [ 0 , 1 ] U[0,1] U[0,1]的随机变量 X X X生成相应的随机变量 Y Y Y。
P ( F − 1 ( X ) < x 0 ) = P ( X < F ( x 0 ) ) = F ( x 0 ) P(F^{-1}(X)<x_0)=P(X<F(x_0))=F(x_0) P(F−1(X)<x0​)=P(X<F(x0​))=F(x0​)
因而随机变量 Y = F − 1 ( X ) Y=F^{-1}(X) Y=F−1(X)即为服从分布函数为 F F 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 ] × U [ 0 , 1 ] U[0,1]×U[0,1] U[0,1]×U[0,1]的二维随机变量，检验落点的位置：

如果落点在曲线下方，则记数变量 M M M加一，否则不加。
那么根据大数定律，应该有 l i m N → ∞ M N = I n t ( e − x ) S [ 0 , 1 ] × [ 0 , 1 ] lim_{N\rightarrow\infty}\frac{M}{N}=\frac{Int(e^{-x})}{S_{[0,1]×[0,1]}} limN→∞​NM​=S[0,1]×[0,1]​Int(e−x)​
我们将这个想法代码化，并将结果附后：

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()