统计信号处理实验二

1、精品文档. 统计信号处理实验二一、实验目的：1、掌握参数估计方法；2、掌握用计算机分析数据的方法。二、实验内容：假设一个运动目标，在外力作用下作一维匀加速运动。其运动轨迹满足的方程为：2001( )2s tatv ts。其中a为目标的加速度，0v为 t=0 时目标运动的速度（初速度） ，0s为目标在t=0 时的初始位置。对目标位置的观测结果为：( )( )( )x ts tn t其中( )x t为观测到的目标位置，( )(0,1)n tn， 为白色观测噪声。 假设在 t=0， 1， 2， ，99s时刻分别取得了100个观测结果x(0)， x(1)， ，x(99)。1、分别用最大似然、最小二乘方

2、法，根据观测结果求出a，0v和0s；2、用 monte_carlo 法，计算出上面两种方法求出的参数的偏差和方差；3、 利用估计出的参数，得到目标位置( )s t的估计( )s t$， 并用 monte_carlo 法计算在t=0，1，2，99s 等各个时间点上对目标位置估计的方差和偏差；4、将噪声的分布改为在（-1，+1）区间的均匀分布，应用最大似然法对参数进行估计，自己推导该分布下的计算公式。三、实验要求：1、设计仿真计算的matlab 程序，给出软件清单；2、完成实验报告，给出实验结果，并对实验结果进行描述，对实验数据进行分析。四、实验过程：假 设s0=0,v0=0.1,a=0.01为

3、实 际 值 ， 将 其 代 入2001( )2s tatv ts， 得 到( )s t。( )( )( )x ts tn t，( )(0,1)n tn为白色观测噪声。假设在t=0，1，2，99s 时刻分别取得了 100 个观测结果x(0)，x(1)，x(99)。1、分别用最大似然、最小二乘方法，根据观测结果求出a，0v和0s；（1）利用最大似然估计（ml估计）的目标是寻找使得先验概率密度函数)|x(p最大的条件作为估计的结果，即?max( | )mlp x。这个数值可以用下面的公式计算出：精品文档. ?(| )0mlp x或?ln( |)0mlp x。即可解出让似然函数取得最大值的$ml。本实

4、验中，2)(211221)|,.,(niiixnnexxxp221121221niniiiniiixxneniininiiiininxxxxxp12112212121)21ln()|,.,(ln其中，2005.0ativsi分别令0ln0sp，0ln0vp，0lnap，niit11，niit122，niit133，niit144，niixtx10_，niiixtx11_，niiixtx122_可以得到矩阵方程：ba，则bat。其中432321215 .02222ttttttttna，avs00，210_2_2_2txtxtxb（2）线性最小二乘法限定 观测结果和待估计参量之间有下列线性 关系：

5、xh n等式中的 h 是根据先验知识已知的参数矩阵，n 是在观测中附加的未知干扰。设计目标就是寻找一个使得观测误差平方和最小的参数矢量作为估计结果。$1ttlsh hh x2、用 monte_carlo 法，计算出上面两种方法求出的参数的偏差和方差；3、利用估计出的参数，得到目标位置( )s t的估计( )s t$，并用 monte_carlo 法计算在t=0，1，2， 99s 等各个时间点上对目标位置估计的方差和偏差；4、将噪声的分布改为在（-1，+1）区间的均匀分布，应用最大似然法对参数进行估计，自己推导该分布下的计算公式。将噪声的分布改为在（-1，+1）区间分布nt= unifrnd(-

6、1,1)，最大似然估计的方法要进行变化：?max(| )mlp x。对于均匀分布的情况，无法写出联合密度函数，计算最大似然比较困难，用正态分布结果进行近似。五、实验结果：精品文档. 1、最大似然和最小二乘法得到的结果如下：可以发现对于估测结果a、v 最好，估测结果s0 最差。从信号变化的角度，或许可以这样理解：随着时间变化，信号发生变化。其中，提供的加速度和初始速度的信息最多，提供的初始位移信息量最少。最大似然估计和最小二乘法相比，结果相差很小，基本符合实际值。精品文档. 2、偏差和方差结果如下：两种方法得到的结果基本是无偏的，并很类似。3、各个时间点上对目标位置估计的方差和偏差如下：从图中可

7、以看出，两种方法下的估计结果偏差很小，方差也不大，估计的效果很不错。4、将程序中estimation函数nt=randn(1,100) 换成nt= unifrnd(-1,1)，重复（ 1）（ 2）精品文档. （3）内容。（1）（2）（3）精品文档. 比较可知，噪声分布改为均匀分布后，参数估计的方差和偏差变小，估计原因是因为-1,1区间均匀分布的噪声方差小于白噪声。六、程序清单clear;%1s0=0;v0=0.1;a=0.01;n=100;figure(1)ml ls=estimation(s0,v0,a,n);subplot(1,2,1);bar(mean(ml);set(gca,xtick

8、label,s0; v0; a);%set 函数将当前图形 (gca)的 x 轴坐标刻度 (xtick)标志为 :s0,v0,atitle( 最大似然估计值 );subplot(1,2,2)bar(mean(ls);set(gca,xticklabel,s0; v0; a);title( 最小二乘法估计值 ); clear;%2精品文档. s0=0;v0=0.1;a=0.01;n=100;ml ls bias_every_ml bias_every_ls variance_every_ml variance_every_ls=estimation(s0,v0,a,n);ml=ml;ls=ls;

9、bias_ml=sum(ml)./n-0,0.1,0.01;variance_ml=var(ml);bias_ls=sum(ls)./n-0,0.1,0.01;variance_ls=var(ls);figure(1);subplot(2,3,1)bar(bias_ml(1,1) bias_ls(1,1); % 取矩阵的第一行第一列，即s0 的偏差set(gca,xticklabel, 最大似然 ; 最小二乘法 );title(s0的偏差 );subplot(2,3,2)bar(bias_ml(1,2) bias_ls(1,2);set(gca,xticklabel, 最大似然 ; 最小二乘法

10、 );title(v0的偏差 );subplot(2,3,3)bar(bias_ml(1,3) bias_ls(1,3);set(gca,xticklabel, 最大似然 ; 最小二乘法 );title(a的偏差 );subplot(2,3,4)bar(variance_ml(1,1) variance_ls(1,1);set(gca,xticklabel, 最大似然 ; 最小二乘法 );title(s0的方差 );subplot(2,3,5)bar(variance_ml(1,2) variance_ls(1,2);set(gca,xticklabel, 最大似然 ; 最小二乘法 );tit

11、le(v0的方差 );subplot(2,3,6)bar(variance_ml(1,3) variance_ls(1,3);set(gca,xticklabel, 最大似然 ; 最小二乘法 );title(a的方差 );clear;%3s0=0;v0=0.1;a=0.01;n=100;精品文档. ml ls bias_every_ml bias_every_ls variance_every_ml variance_every_ls=estimation(s0,v0,a,n);figure(1)subplot(2,2,1)plot(0:99,bias_every_ml)title( 最大似然

12、各点偏差 );subplot(2,2,2)plot(0:99,bias_every_ls)title( 最小二乘法各点偏差 );subplot(2,2,3)plot(0:99,variance_every_ml)title( 最大似然各点方差 );subplot(2,2,4)plot(0:99,variance_every_ls)title( 最小二乘法各点方差 );%estimation函数function theta_ml theta_ls bias_every_ml bias_every_ls variance_every_ml variance_every_ls=estimation(

13、s0,v0,a,n)for j=1:n%-计算初始化 -%t=0:99;nt=randn(1,100);%generate values from a normal distribution with mean 1 and standard deviation 2.r = 1 + 2.*randn(100,1);%nt 均值为 0，方差为1st=s0+v0*t+0.5*a*t.2;xt=st+nt;t1=sum(t);t2=t*t;t3=t.2*t;t4=(t.2)*(t.2);x_t0=sum(xt);x_t1=sum(xt.*t);x_t2=sum(xt.*t.2);a=200 2*t1 t2; 2*t1 2*t2 t3; t2 t3 0.5*t4;b=(2*x_t0 2*x_t1 x_t2);theta_ml(:,j)=(inv(a)*b); % 表示求转置， inv表示求逆矩阵， theta_ml(:,j)表示取出第j 列%-ls计算 -%h_s=ones(100,1);h_v=t;h_a=(0.5*t.2);精品文档. h=h_s