东南大学统计信号处理实验二

《统计信号处理》实验二实验目的：1．掌握参数估计方法；2．掌握用计算机分析数据的方法。实验内容：假设一个运动目标，在外力作用下作一维匀加速运动。其运动轨迹满足的方程为：0221)(svtatts。其中a为目标的加速度，v为t=0时目标运动的速度（初速度），0s为目标在t=0时的初始位置。对目标位置的观测结果为：()()()xtstnt其中()xt为观测到的目标位置，()(0,1)ntN，为白色观测噪声。假设在t=0,1,2,…,99s时刻分别取得了100个观测结果x(0),x(1),…,x(99)。1）分别用最大似然，最小二乘方法，根据观测结果求出a，v和0s；2）用Monte_Carlo法，计算出上面两种方法求出的参数的偏差和方差；3）利用估计出的参数，得到目标位置的时间参数()st的估计()st，并用Monte_Carlo法计算在t=0,1,2,…,99s等各个时间点上对目标位置估计的方差和偏差；4）将噪声的分布改为在（-1，+1）区间分布，应用上面推导出的最大似然，最小二乘公式对参数进行估计，并计算估计的偏差和方差。实验要求：1)设计仿真计算的Matlab程序，给出软件清单；2)完成实验报告，对实验结果进行描述，并给出实验结果，对实验数据进行分析。实验结果如下：最小二乘法：01020304050607080901000102030405060708090time(t)s(t)&x(t)Theleastsquaremethod——a=0.0098;v=0.35;s0=3.2;delta=1-----------------------------最大似然估计方法----------------------01020304050607080901000102030405060708090time(t)s(t)&x(t)Themostlikelihoodmethod——a=0.0098;v=0.35;s0=3.2;delta=1实验结果分析：可以发现对于估测结果a最好，估测结果v次之，估测结果s0最差。从信号变化的角度，或许可以这样理解：随着时间变化，信号发生变化。其中，提供的加速度的信息最多，其次是初始速度，提供的初始位移信息量最少。最大似然估计和最小二乘法相比之下，最大似然估计结果较好一些，但相差很小，基本符合实际值。（2）利用Monte_Carlo法，计算出上面两种方法求出的参数的偏差和方差实验结果如下：实验结果分析：实际测得两种估计方法，得到的结果基本上是无偏的。同等观测条件下，两种方法的性能是一致的。（3）利用估计出的参数，得到目标位置的时间参数()st的估计()st，并用Monte_Carlo法计算在t=0,1,2,…,99s等各个时间点上对目标位置估计的方差和偏差；实验结果如下：实验结果分析：从图中可以看出，两种方法下的估计结果偏差很小，方差也不大，估计的效果很不错。（4）将噪声的分布改为在（-1，+1）区间分布w=2*delta*rand(1,100)-delta;当噪声为均匀分布时，最小二乘公式不需要改变，但是最大似然估计的方法要进行变化：)|(maxˆxpml对于均匀分布的情况，联合密度函数为二值函数，计算最大似然比较困难，近似用正态分布结果进行近似。实验结果分析：实际测得，这种情况下，两种估计方法，得到的结果也基本上是无偏的。同等观测条件下，两种方法的性能趋于一致。软件清单：Estimate_show.ms0=0;v0=0.1;a=0.01;N=100;%1figure(1)[mlls]=estimation(s0,v0,a,N);subplot(1,2,1)bar(mean(ml'));set(gca,'XTickLabel',{'s0';'v0';'a'});title({'最大似然估计值'});subplot(1,2,2)bar(mean(ls'));set(gca,'XTickLabel',{'s0';'v0';'a'});title({'最小二乘法估计值'});%2[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(2);subplot(2,3,1)bar([bias_ml(1,1)bias_ls(1,1)]);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的方差'});%3figure(3)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.mfunction[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);st=s0+v0*t+0.5*a*t.^2;xt=st+nt;%%%%%%%%%%%%--ML计算--%%%%%%%%%%%%%%%%%%%%%计算ML参数[s0v0a]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);%%%%%%%%%%%%--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'))');