【统计信号处理kay 第七章 两道习题python仿真】

---

前言

本次仿真参考教材Kay Fundamentals Of Statistical Signal Processing - Estimation Theory第七章7.13和7.14题
学艺不精，希望大家多多包含，提供合理的指点

---

一、蒙特卡洛计算法

内容摘抄自书中7A的附录，介绍了一种近似估计量pdf的方式，作为后续仿真的参考方法

二、7.13题

1.题目内容

关于WGN中的DC电平信号，书中有详细的讨论。此外，这里涉及到样本均值的分布问题，如何求样本均值的期望和方差。这里我们无需使用7A提供的fortran程序，使用python进行仿真

2.流程图

3.示例代码

代码如下：

import numpy as np
import matplotlib.pyplot as plt
import random
########################################
# name:w_gen
# function：生成正态分布的w向量
# parameter:
#   N:样本长度
#   mu:均值
#   sigma:标准差
# return: w向量
def w_gen(N, mu, sigma):
    w = []
    for i in range(N):
        w.append(random.gauss(mu, sigma))
    return w
    
#########################################
# name:guass
# function：生成正态分布的pdf
# parameter:
#   n:坐标轴向量或列表
#   mu:均值
#   sigma2:方差
# return: 正态分布pdf向量
def guass(n, mu, sigma2):
    return (1 / np.sqrt(2 * np.pi * sigma2)) * np.exp((-(n - mu) ** 2) / (2 * sigma2))

if __name__ == '__main__':
    # #########################################数据定义##################################
    N = 50                                                   # 样本长度
    mu = 0                                                   # w(n)均值
    sigma2 = 0.1                                             # w(n)方差
    A = 1                                                    # A
    M = 5000                                                 # 蒙特卡洛计算次数
    max_x = 1.2                                              # 坐标轴上界
    min_x = 0.8                                              # 坐标轴下界
    n = np.arange(min_x, max_x, (max_x - min_x) / M)         # 理论pdf坐标轴生成
    A_ba = []                                                # 样本均值列表

    pdf_label = []                                           # 近似pdf坐标轴
    x_mean = 0
    # #########################################生成理论pdf###############################

    pdf = guass(n, A, sigma2 / N)
    # #########################################蒙特卡洛计算###############################
    for j in range(M):  # 循环M次
        # 生成样本x向量
        x = w_gen(N, mu, np.sqrt(sigma2)) + A * np.ones(N)
        # 计算样本均值
        x_mean = np.mean(x)
        # 加入样本均值列表
        A_ba.append(x_mean)

    # ##########################################计算样本均值的期望和方差####################
    EA_ba = np.mean(A_ba)
    varA_ba = 0
    # 样本方差计算
    for i in range(M):
        varA_ba = varA_ba + ((A_ba[i] - EA_ba) ** 2)
    varA_ba /= (M - 1)

    # ###########################################直方图数据统计###########################
    # bins为range=[min_x,max_x]内的区间数量，根据附录7A，注意控制区间数和M使得近似pdf达到收敛
    bins = 200
    # A_ba_pdf为直方图数据，bin_edges为直方图每个区间上下边界的坐标，因此bin_edges.size==bins+1
    A_ba_pdf, bin_edges = np.histogram(A_ba, bins=bins, range=[min_x, max_x])
    A_ba_pdf = A_ba_pdf.astype(np.float64)

    # ############################################近似pdf###############################
    for i in range(len(bin_edges) - 1):
        deltax = (bin_edges[i + 1] - bin_edges[i])
        # 附录7A提到的PDF近似方法
        A_ba_pdf[i] = (A_ba_pdf[i] / M) / deltax
        # 坐标取区间的中点
        pdf_label.append((bin_edges[i + 1] + bin_edges[i]) / 2)

    # ############################################对比图################################
    plt.scatter(pdf_label, A_ba_pdf, label='PDF of sample mean', s=10, c='r')
    plt.plot(n, pdf, label='PDF of normal')
    plt.xlabel('n')
    plt.ylabel('p')
    plt.title('M=' + str(M) + ' E=' + str(EA_ba) + ' var=' + str(varA_ba))
    plt.legend()
    plt.show()

4.实验结果

理论分析

摘抄自我的实验报告

三、7.14题

1.题目内容

本题的重点是分析估计量的PDF，作为初学者很难对均值除以标准差这样形式的估计量的PDF有概念，虽然我无法从计算费雪信息，去直接证明它的渐近性，但是可以联想到考研数学里的chi分布和t分布，从而确定理论PDF，能够证明N足够大时服从标准正态分布。

2.参考代码

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


########################################
# name:w_gen
# function：生成正态分布的w向量
# parameter:
#   N:样本长度
#   mu:均值
#   sigma:标准差
# return: w向量
def w_gen(N, mu, sigma):
    w = []
    for j in range(N):
        w.append(random.gauss(mu, sigma))
    return w


#########################################
# name:guass
# function：生成正态分布的pdf
# parameter:
#   n:坐标轴向量或列表
#   mu:均值
#   sigma2:方差
# return: 正态分布pdf向量
def guass(n, mu, sigma2):
    return (1 / np.sqrt(2 * np.pi * sigma2)) * np.exp((-(n - mu) ** 2) / (2 * sigma2))


if __name__ == '__main__':
    # #########################################数据定义##################################
    N = 10                                                     # 样本长度
    M = 5000                                                   # 蒙特卡洛计算次数
    mu = 0                                                     # w(n)均值
    sigma2 = 1                                                 # 方差
    max_x = 5                                                  # 坐标轴上界
    min_x = -5                                                 # 坐标轴下界
    n = np.arange(min_x, max_x, (max_x - min_x) / M)           # 理论pdf坐标轴生成
    A_ba = []                                                  # 估计量列表
    pdf = np.empty(n.size)                                     # 理论pdf数据列表
    pdf_label = []                                             # 近似pdf坐标轴
    sigma2_ba = 0                                              # 估计方差
    # #####################################生成理论pdf标准高斯#############################
    pdf = guass(n, 0, 1)
    # #########################################蒙特卡洛计算###############################
    for i in range(M):
        # 生成样本x向量
        x = w_gen(N, mu, sigma2)
        # 计算样本均值
        x_mean = np.mean(x)
        # 计算统计量
        for j in range(N):
            sigma2_ba = sigma2_ba + ((x[j] - x_mean) ** 2)
        sigma2_ba = sigma2_ba / N
        # 统计量存入
        A_ba.append(x_mean * np.sqrt(N) / np.sqrt(sigma2_ba))
    # ##########################################计算统计量的均值和方差#####################
    EA_ba = np.mean(A_ba)
    varA_ba = 0
    # 样本方差计算
    for i in range(M):
        varA_ba = varA_ba + ((A_ba[i] - EA_ba) ** 2)
    varA_ba /= (M - 1)
    # ###########################################直方图数据统计###########################
    # bins为range=[min_x,max_x]内的区间数量，根据附录7A，注意控制区间数和M使得近似pdf达到收敛
    bins = 200
    A_ba_pdf, bin_edges = np.histogram(A_ba, bins=bins, range=[min_x, max_x])
    A_ba_pdf = A_ba_pdf.astype(np.float64)

    # ############################################近似pdf###############################
    for i in range(len(bin_edges) - 1):
        deltax = (bin_edges[i + 1] - bin_edges[i])
        # 附录7A提到的PDF近似方法
        A_ba_pdf[i] = A_ba_pdf[i] / M / deltax
        # 坐标取区间的中点
        pdf_label.append((bin_edges[i + 1] + bin_edges[i]) / 2)

    plt.plot(pdf_label, A_ba_pdf, label='PDF of estimate')
    plt.plot(n, pdf,label='PDF of normal')
    plt.xlabel('n')
    plt.ylabel('p')
    plt.title('N=' + str(N) + ' E=' + str(EA_ba) + ' var=' + str(varA_ba))
    plt.legend()
    plt.show()

3.理论分析