哈爾濱工程大學2024年隨機過程上機作業?一、作業目的本次隨機過程上機作業旨在幫助同學們更好地理解和掌握隨機過程的基本概念、理論和方法，通過實際編程操作，提高同學們運用隨機過程知識解決實際問題的能力，培養同學們的編程實踐能力和數據分析能力。二、作業要求（一）題目1：泊松過程的模擬1.生成強度為\(\lambda\)的泊松過程\(\{N(t),t\geq0\}\)在\([0,T]\)上的樣本路徑。2.計算并繪制該泊松過程在\([0,T]\)上的均值函數\(E[N(t)]\)和方差函數\(Var[N(t)]\)。（二）題目2：馬爾可夫鏈的模擬1.給定一個有限狀態的馬爾可夫鏈，其狀態空間為\(S=\{1,2,\cdots,n\}\)，轉移概率矩陣為\(P=(p_{ij})\)。2.生成該馬爾可夫鏈的一條長度為\(m\)的樣本路徑。3.計算并繪制該馬爾可夫鏈的狀態轉移頻率直方圖。（三）題目3：平穩隨機過程的功率譜估計1.生成一個零均值的平穩隨機過程\(\{X(t),t\inR\}\)，其自相關函數為\(R_X(\tau)\)。2.采用周期圖法和韋爾奇法估計該平穩隨機過程的功率譜密度\(S_X(f)\)。3.繪制真實的自相關函數\(R_X(\tau)\)和估計的功率譜密度\(S_X(f)\)。三、相關理論知識（一）泊松過程1.定義：設\(\{N(t),t\geq0\}\)是一個隨機過程，如果它滿足以下條件：\(N(0)=0\)；\(N(t)\)是獨立增量過程；對于任意的\(s,t\geq0\)，\(N(t+s)N(s)\)服從參數為\(\lambdat\)的泊松分布，即\(P\{N(t+s)N(s)=k\}=\frac{(\lambdat)^k}{k!}e^{\lambdat},k=0,1,2,\cdots\)，則稱\(\{N(t),t\geq0\}\)是強度為\(\lambda\)的泊松過程。2.均值函數和方差函數：對于泊松過程\(\{N(t),t\geq0\}\)，其均值函數\(E[N(t)]=\lambdat\)，方差函數\(Var[N(t)]=\lambdat\)。（二）馬爾可夫鏈1.定義：設\(\{X_n,n=0,1,2,\cdots\}\)是一個隨機過程，狀態空間為\(S\)，如果對于任意的\(n\geq0\)和\(i_0,i_1,\cdots,i_{n+1}\inS\)，有\(P\{X_{n+1}=i_{n+1}|X_0=i_0,X_1=i_1,\cdots,X_n=i_n\}=P\{X_{n+1}=i_{n+1}|X_n=i_n\}\)，則稱\(\{X_n,n=0,1,2,\cdots\}\)是一個馬爾可夫鏈。2.轉移概率矩陣：設馬爾可夫鏈\(\{X_n,n=0,1,2,\cdots\}\)的狀態空間為\(S=\{1,2,\cdots,n\}\)，則稱矩陣\(P=(p_{ij})\)為轉移概率矩陣，其中\(p_{ij}=P\{X_{n+1}=j|X_n=i\}\)。（三）平穩隨機過程1.定義：設\(\{X(t),t\inR\}\)是一個隨機過程，如果對于任意的\(t_1,t_2,\cdots,t_n\inR\)和\(\tau\inR\)，\((X(t_1),X(t_2),\cdots,X(t_n))\)與\((X(t_1+\tau),X(t_2+\tau),\cdots,X(t_n+\tau))\)具有相同的聯合分布，則稱\(\{X(t),t\inR\}\)是平穩隨機過程。2.自相關函數：平穩隨機過程\(\{X(t),t\inR\}\)的自相關函數定義為\(R_X(\tau)=E[X(t)X(t+\tau)]\)，它是\(\tau\)的偶函數，且\(R_X(0)=E[X^2(t)]\)為隨機過程的平均功率。3.功率譜密度：平穩隨機過程\(\{X(t),t\inR\}\)的功率譜密度\(S_X(f)\)與自相關函數\(R_X(\tau)\)是一對傅里葉變換對，即\(S_X(f)=\int_{\infty}^{\infty}R_X(\tau)e^{j2\pif\tau}d\tau\)，\(R_X(\tau)=\int_{\infty}^{\infty}S_X(f)e^{j2\pif\tau}df\)。四、具體實現（一）泊松過程的模擬1.算法思路：根據泊松過程的定義，利用泊松分布的隨機數生成方法生成\(N(t)\)在不同時刻的值。計算均值函數\(E[N(t)]=\lambdat\)和方差函數\(Var[N(t)]=\lambdat\)。2.代碼實現（Python）：```pythonimportnumpyasnpimportmatplotlib.pyplotasplt參數設置lambda_=3T=10t=np.linspace(0,T,100)生成泊松過程樣本路徑N_t=np.zeros(len(t))foriinrange(1,len(t)):N_t[i]=N_t[i1]+np.random.poisson(lambda_*(t[i]t[i1]))計算均值函數和方差函數mean_N_t=lambda_*tvar_N_t=lambda_*t繪制樣本路徑plt.figure(figsize=(10,5))plt.plot(t,N_t,label='PoissonProcessSamplePath')plt.plot(t,mean_N_t,'r',label='MeanFunction')plt.plot(t,var_N_t,'g',label='VarianceFunction')plt.xlabel('t')plt.ylabel('N(t)')plt.title('PoissonProcessSimulation')plt.legend()plt.grid(True)plt.show()```3.結果分析：從生成的樣本路徑可以直觀地看到泊松過程的隨機特性。均值函數和方差函數的理論值與繪制的曲線完全吻合，驗證了泊松過程的均值和方差公式。（二）馬爾可夫鏈的模擬1.算法思路：根據給定的轉移概率矩陣，利用隨機數生成方法確定馬爾可夫鏈的狀態轉移。統計狀態轉移頻率，繪制狀態轉移頻率直方圖。2.代碼實現（Python）：```pythonimportnumpyasnpimportmatplotlib.pyplotasplt狀態空間和轉移概率矩陣n=5P=np.array([[0.1,0.2,0.3,0.2,0.2],[0.3,0.1,0.2,0.2,0.2],[0.2,0.2,0.1,0.3,0.2],[0.2,0.2,0.3,0.1,0.2],[0.2,0.2,0.2,0.2,0.1]])生成馬爾可夫鏈樣本路徑m=100X=np.zeros(m,dtype=int)X[0]=np.random.choice(n)foriinrange(1,m):X[i]=np.random.choice(n,p=P[X[i1]])計算狀態轉移頻率transition_count=np.zeros((n,n))foriinrange(m1):transition_count[X[i]][X[i+1]]+=1繪制狀態轉移頻率直方圖plt.figure(figsize=(10,5))foriinrange(n):forjinrange(n):plt.bar((i,j),transition_count[i][j],width=0.2,bottom=np.sum(transition_count[:i,j]),label=f'From{i}to{j}'ifi==0else"")plt.xlabel('FromState')plt.ylabel('TransitionCount')plt.title('MarkovChainStateTransitionFrequency')plt.legend()plt.grid(True)plt.show()```3.結果分析：生成的樣本路徑展示了馬爾可夫鏈的狀態轉移情況。狀態轉移頻率直方圖反映了不同狀態之間的轉移概率，與給定的轉移概率矩陣趨勢相符。（三）平穩隨機過程的功率譜估計1.算法思路：生成零均值平穩隨機過程的樣本數據。分別采用周期圖法和韋爾奇法估計功率譜密度。計算真實的自相關函數，并與估計的功率譜密度進行傅里葉變換對比。2.代碼實現（Python）：```pythonimportnumpyasnpimportmatplotlib.pyplotaspltfromscipy.signalimportwelch生成平穩隨機過程樣本數據N=1024t=np.linspace(0,1,N)X=np.random.randn(N)R_X=np.correlate(X,X,mode='full')/NR_X=R_X[R_X.size//2:]周期圖法估計功率譜密度f1,S1=plt.psd(X,NFFT=N,Fs=1/(t[1]t[0]))韋爾奇法估計功率譜密度f2,S2=welch(X,fs=1/(t[1]t[0]))計算真實的自相關函數的傅里葉變換S_true=np.fft.fft(R_X)/Nf_true=np.fft.fftfreq(N,t[1]t[0])繪制結果plt.figure(figsize=(10,5))plt.plot(f_true[:N//2],np.abs(S_true[:N//2]),label='TruePowerSpectrum')plt.plot(f1,S1,label='Periodogram')plt.plot(f2,S2,label='WelchMethod')plt.xlabel('Frequency(Hz)')plt.ylabel('PowerSpectrumDensity')plt.title('PowerSpectrumEstimation')plt.legend()plt.grid(True)plt.show()```3.結果分析：周期圖法和韋爾奇法都能對平穩隨機過程的功率譜密度進行估計。真實的功率譜密度與兩種估計方法的結果在趨勢上基本一致，但周期圖法可能存在