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

《统计信号处理》实验二 实验目的：
1．掌握参数估计方法；
2．掌握用计算机分析数据的方法。

实验内容：
假设一个运动目标，在外力作用下作一维匀加速运动。

其运动轨迹满足的方程为：022
1)(s vt at t s ++=。

其中a 为目标的加速度，v 为t=0时目标运动的速度（初速度），0s 为目标在t=0时的初始位置。

对目标位置的观测结果为：
()()()
x t s t n t =+ 其中()x t 为观测到的目标位置，()(0,1)n t N ∈，为白色观测噪声。

假设在t=0,1,2,…,99s 时刻分别取得了100个观测结果x(0),x(1),…,x(99)。

1） 分别用最大似然，最小二乘方法，根据观测结果求出a ，v 和0s ；
2） 用Monte_Carlo 法，计算出上面两种方法求出的参数的偏差和方差；
3）利用估计出的参数，得到目标位置的时间参数()s t 的估计()s
t ，并用Monte_Carlo 法计算在t=0,1,2,…,99s 等各个时间点上对目标位置估计的方差和偏差；
4） 将噪声的分布改为在（-1，+1）区间分布，应用上面推导出的最大似然，最小二乘公式对参数进行估计，并计算估计的偏差和方差。

实验要求：
1)设计仿真计算的Matlab 程序，给出软件清单；
2)完成实验报告，对实验结果进行描述，并给出实验结果，对实验数据进行分析。

实验结果如下：
最小二乘法：
010********
60708090100010
20
30
40
5060
70
80
90
time(t)s (t )&x (t )The least square method——a=0.0098;v=0.35;s0=3.2;delta=1
-----------------------------最大似然估计方法----------------------
010********
60708090100
010
20
30
40
5060
70
80
90
time(t)s (t )&x (t )The most likelihood method——a=0.0098;v=0.35;s0=3.2;delta=1
实验结果分析：
可以发现对于估测结果a最好，估测结果v次之，估测结果s0最差。

从信号变化的角度，或许可以这样理解：随着时间变化，信号发生变化。

其中，提供的加速度的信息最多，其次是初始速度，提供的初始位移信息量最少。

最大似然估计和最小二乘法相比之下，最大似然估计结果较好一些，但相差很小，基本符合实际值。

（2）利用Monte_Carlo法，计算出上面两种方法求出的参数的偏差和方差
实验结果如下：
实验结果分析：
实际测得两种估计方法，得到的结果基本上是无偏的。

同等观测条件下，两种方法的性能是一致的。

（3）利用估计出的参数，得到目标位置的时间参数()s t 的估计()s
t ，并用Monte_Carlo 法计算在t=0,1,2,…,99s 等各个时间点上对目标位置估计的方差和偏差；
实验结果如下：
实验结果分析：
从图中可以看出，两种方法下的估计结果偏差很小，方差也不大，估计的效果很不错。

（4）将噪声的分布改为在（-1，+1）区间分布
w=2*delta*rand(1,100)-delta;
当噪声为均匀分布时，最小二乘公式不需要改变，但是最大似然估计的方法要进行变化：
)|(max ˆθθθ
x p ml = 对于均匀分布的情况，联合密度函数为二值函数，计算最大似然比较困难，近似用正态分布结果进行近似。

实验结果分析：
实际测得，这种情况下，两种估计方法，得到的结果也基本上是无偏的。

同等观测条件下，两种方法的性能趋于一致。

软件清单：
Estimate_show.m
s0=0;
v0=0.1;
a=0.01;
N=100;
%1
figure(1)
[ml ls]=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
[ml ls bias_every_ml bias_every_ls variance_every_ml variance_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 的方差'});
%3
figure(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.m
function [theta_ml theta_ls bias_every_ml bias_every_ls variance_every_ml variance_every_ls]=estimation(s0,v0,a,N)
for j=1:N
%%%%%%%%%%%%--计算初始化--%%%%%%%%%%%%%%%%%%%%
t = 0:99;
nt = randn(1,100);
st = s0 + v0*t + 0.5*a*t.^2;
xt = st + nt;
%%%%%%%%%%%%--ML计算--%%%%%%%%%%%%%%%%%%%%
%计算ML参数[s0 v0 a]
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);
%%%%%%%%%%%%--LS计算--%%%%%%%%%%%%%%%%%%%%
h_s=ones(100,1);
h_v=t';
h_a=(0.5*t.^2)';
H=[h_s h_v h_a];
theta_ls(:,j)=inv(H'*H)*H'*xt';
end
xt_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'))');。