哈爾濱工程大學2024年隨機過程上機作業(yè)?一、作業(yè)目的本次隨機過程上機作業(yè)旨在幫助同學們更好地理解和掌握隨機過程的基本概念、理論和方法，通過實際編程操作，提高同學們運用隨機過程知識解決實際問題的能力，培養(yǎng)同學們的編程實踐能力和數(shù)據分析能力。二、作業(yè)要求（一）題目1：泊松過程的模擬1.生成強度為\(\lambda\)的泊松過程\(\{N(t),t\geq0\}\)在\([0,T]\)上的樣本路徑。2.計算并繪制該泊松過程在\([0,T]\)上的均值函數(shù)\(E[N(t)]\)和方差函數(shù)\(Var[N(t)]\)。（二）題目2：馬爾可夫鏈的模擬1.給定一個有限狀態(tài)的馬爾可夫鏈，其狀態(tài)空間為\(S=\{1,2,\cdots,n\}\)，轉移概率矩陣為\(P=(p_{ij})\)。2.生成該馬爾可夫鏈的一條長度為\(m\)的樣本路徑。3.計算并繪制該馬爾可夫鏈的狀態(tài)轉移頻率直方圖。（三）題目3：平穩(wěn)隨機過程的功率譜估計1.生成一個零均值的平穩(wěn)隨機過程\(\{X(t),t\inR\}\)，其自相關函數(shù)為\(R_X(\tau)\)。2.采用周期圖法和韋爾奇法估計該平穩(wěn)隨機過程的功率譜密度\(S_X(f)\)。3.繪制真實的自相關函數(shù)\(R_X(\tau)\)和估計的功率譜密度\(S_X(f)\)。三、相關理論知識（一）泊松過程1.定義：設\(\{N(t),t\geq0\}\)是一個隨機過程，如果它滿足以下條件：\(N(0)=0\)；\(N(t)\)是獨立增量過程；對于任意的\(s,t\geq0\)，\(N(t+s)N(s)\)服從參數(shù)為\(\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.均值函數(shù)和方差函數(shù)：對于泊松過程\(\{N(t),t\geq0\}\)，其均值函數(shù)\(E[N(t)]=\lambdat\)，方差函數(shù)\(Var[N(t)]=\lambdat\)。（二）馬爾可夫鏈1.定義：設\(\{X_n,n=0,1,2,\cdots\}\)是一個隨機過程，狀態(tài)空間為\(S\)，如果對于任意的\(n\geq0\)和\(i_0,i_1,\cdots,i_{n+1}\inS\)，有\(zhòng)(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\}\)的狀態(tài)空間為\(S=\{1,2,\cdots,n\}\)，則稱矩陣\(P=(p_{ij})\)為轉移概率矩陣，其中\(zhòng)(p_{ij}=P\{X_{n+1}=j|X_n=i\}\)。（三）平穩(wěn)隨機過程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))\)具有相同的聯(lián)合分布，則稱\(\{X(t),t\inR\}\)是平穩(wěn)隨機過程。2.自相關函數(shù)：平穩(wěn)隨機過程\(\{X(t),t\inR\}\)的自相關函數(shù)定義為\(R_X(\tau)=E[X(t)X(t+\tau)]\)，它是\(\tau\)的偶函數(shù)，且\(R_X(0)=E[X^2(t)]\)為隨機過程的平均功率。3.功率譜密度：平穩(wěn)隨機過程\(\{X(t),t\inR\}\)的功率譜密度\(S_X(f)\)與自相關函數(shù)\(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\)。四、具體實現(xiàn)（一）泊松過程的模擬1.算法思路：根據泊松過程的定義，利用泊松分布的隨機數(shù)生成方法生成\(N(t)\)在不同時刻的值。計算均值函數(shù)\(E[N(t)]=\lambdat\)和方差函數(shù)\(Var[N(t)]=\lambdat\)。2.代碼實現(xiàn)（Python）：```pythonimportnumpyasnpimportmatplotlib.pyplotasplt參數(shù)設置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]))計算均值函數(shù)和方差函數(shù)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.結果分析：從生成的樣本路徑可以直觀地看到泊松過程的隨機特性。均值函數(shù)和方差函數(shù)的理論值與繪制的曲線完全吻合，驗證了泊松過程的均值和方差公式。（二）馬爾可夫鏈的模擬1.算法思路：根據給定的轉移概率矩陣，利用隨機數(shù)生成方法確定馬爾可夫鏈的狀態(tài)轉移。統(tǒng)計狀態(tài)轉移頻率，繪制狀態(tài)轉移頻率直方圖。2.代碼實現(xiàn)（Python）：```pythonimportnumpyasnpimportmatplotlib.pyplotasplt狀態(tài)空間和轉移概率矩陣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]])計算狀態(tài)轉移頻率transition_count=np.zeros((n,n))foriinrange(m1):transition_count[X[i]][X[i+1]]+=1繪制狀態(tài)轉移頻率直方圖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.結果分析：生成的樣本路徑展示了馬爾可夫鏈的狀態(tài)轉移情況。狀態(tài)轉移頻率直方圖反映了不同狀態(tài)之間的轉移概率，與給定的轉移概率矩陣趨勢相符。（三）平穩(wěn)隨機過程的功率譜估計1.算法思路：生成零均值平穩(wěn)隨機過程的樣本數(shù)據。分別采用周期圖法和韋爾奇法估計功率譜密度。計算真實的自相關函數(shù)，并與估計的功率譜密度進行傅里葉變換對比。2.代碼實現(xiàn)（Python）：```pythonimportnumpyasnpimportmatplotlib.pyplotaspltfromscipy.signalimportwelch生成平穩(wěn)隨機過程樣本數(shù)據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]))計算真實的自相關函數(shù)的傅里葉變換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.結果分析：周期圖法和韋爾奇法都能對平穩(wěn)隨機過程的功率譜密度進行估計。真實的功率譜密度與兩種估計方法的結果在趨勢上基本一致，但周期圖法可能存在