统计信号处理实验二

《统计信号处理》实验二一、实验目的：1、掌握参数估计方法；2、掌握用计算机分析数据的方法。二、实验内容：假设一个运动目标，在外力作用下作一维匀加速运动。其运动轨迹满足的方程为：2001()2statvts。其中a为目标的加速度，0v为t=0时目标运动的速度（初速度），0s为目标在t=0时的初始位置。对目标位置的观测结果为：()()()xtstnt其中()xt为观测到的目标位置，()(0,1)ntN，为白色观测噪声。假设在t=0，1，2，…，99s时刻分别取得了100个观测结果x(0)，x(1)，…，x(99)。1、分别用最大似然、最小二乘方法，根据观测结果求出a，0v和0s；2、用Monte_Carlo法，计算出上面两种方法求出的参数的偏差和方差；3、利用估计出的参数，得到目标位置()st的估计()st，并用Monte_Carlo法计算在t=0，1，2，…，99s等各个时间点上对目标位置估计的方差和偏差；4、将噪声的分布改为在（-1，+1）区间的均匀分布，应用最大似然法对参数进行估计，自己推导该分布下的计算公式。三、实验要求：1、设计仿真计算的Matlab程序，给出软件清单；2、完成实验报告，给出实验结果，并对实验结果进行描述，对实验数据进行分析。四、实验过程：假设s0=0,v0=0.1,a=0.01为实际值，将其代入2001()2statvts，得到()st。()()()xtstnt，()(0,1)ntN为白色观测噪声。假设在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。本实验中，2)(211221)|,...,,(NiiixNNexxxP221121221NiNiiiNiiixxNeNiiNiNiiiiNiNxxxxxP12112212121)21ln()|,...,,(ln其中，2005.0ativsi分别令0ln0sP，0ln0vP，0lnaP，Niit11，Niit122，Niit133，Niit144，Niixtx10_，Niiixtx11_，Niiixtx122_可以得到矩阵方程：bA，则bAT。其中432321215.02222ttttttttNA，avs00，210_2_2_2txtxtxb（2）线性最小二乘法限定观测结果和待估计参量之间有下列线性关系：xHθn等式中的H是根据先验知识已知的参数矩阵，n是在观测中附加的未知干扰。设计目标就是寻找一个使得观测误差平方和最小的参数矢量作为估计结果。1TTLSθHHHx2、用Monte_Carlo法，计算出上面两种方法求出的参数的偏差和方差；3、利用估计出的参数，得到目标位置()st的估计()st，并用Monte_Carlo法计算在t=0，1，2，…，99s等各个时间点上对目标位置估计的方差和偏差；4、将噪声的分布改为在（-1，+1）区间的均匀分布，应用最大似然法对参数进行估计，自己推导该分布下的计算公式。将噪声的分布改为在（-1，+1）区间分布nt=unifrnd(-1,1)，最大似然估计的方法要进行变化：ˆmax(|)mlpx。对于均匀分布的情况，无法写出联合密度函数，计算最大似然比较困难，用正态分布结果进行近似。五、实验结果：1、最大似然和最小二乘法得到的结果如下：可以发现对于估测结果a、v最好，估测结果s0最差。从信号变化的角度，或许可以这样理解：随着时间变化，信号发生变化。其中，提供的加速度和初始速度的信息最多，提供的初始位移信息量最少。最大似然估计和最小二乘法相比，结果相差很小，基本符合实际值。2、偏差和方差结果如下：两种方法得到的结果基本是无偏的，并很类似。3、各个时间点上对目标位置估计的方差和偏差如下：从图中可以看出，两种方法下的估计结果偏差很小，方差也不大，估计的效果很不错。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)[mlls]=estimation(s0,v0,a,N);subplot(1,2,1);bar(mean(ml'));set(gca,'XTickLabel',{'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;%2s0=0;v0=0.1;a=0.01;N=100;[mllsbias_every_mlbias_every_lsvariance_every_mlvariance_every_ls]=estimation(s0,v0,a,N);ml=ml';ls=ls';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',{'最大似然';'最小二乘法'});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',{'最大似然';'最小二乘法'});title({'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;[mllsbias_every_mlbias_every_lsvariance_every_mlvariance_every_ls]=estimation(s0,v0,a,N);figure(1)subplot(2,2,1)plot(0:99,bias_every_ml)title('最大似然各点偏差');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_mltheta_lsbias_every_mlbias_every_lsvariance_every_mlvariance_every_ls]=estimation(s0,v0,a,N)forj=1:N%%%%%%%%%%%%--计算初始化--%%%%%%%%%%%%%%%%%%%%t=0:99;nt=randn(1,100);%Generatevaluesfromanormaldistributionwithmean1andstandarddeviation2.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=[2002*t1t2;2*t12*t2t3;t2t30.5*t4;];b=([2*x_t02*x_t1x_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_sh_vh_a];theta_ls(:,j)=inv(H'*H)*H'*xt';endxt_ba_ml=H*theta_ml;xt_ba_ls=H*theta_ls;bias_every_ml=((sum(xt_ba_ml'))')/N-st';bias_every_ls=((sum(xt_ba_ls'))')/N-st';variance_every_ml=((var(xt_ba_ml'))');variance_every_ls=((var(xt_ba_ls'))');