統(tǒng)計(jì)信號(hào)處理實(shí)驗(yàn)二

./《統(tǒng)計(jì)信號(hào)處理》實(shí)驗(yàn)二一、實(shí)驗(yàn)?zāi)康模?、掌握參數(shù)估計(jì)方法；2、掌握用計(jì)算機(jī)分析數(shù)據(jù)的方法。二、實(shí)驗(yàn)容：假設(shè)一個(gè)運(yùn)動(dòng)目標(biāo),在外力作用下作一維勻加速運(yùn)動(dòng)。其運(yùn)動(dòng)軌跡滿足的方程為：。其中為目標(biāo)的加速度,為t=0時(shí)目標(biāo)運(yùn)動(dòng)的速度〔初速度,為目標(biāo)在t=0時(shí)的初始位置。對(duì)目標(biāo)位置的觀測(cè)結(jié)果為：其中為觀測(cè)到的目標(biāo)位置,,為白色觀測(cè)噪聲。假設(shè)在t=0,1,2,…,99s時(shí)刻分別取得了100個(gè)觀測(cè)結(jié)果x<0>,x<1>,…,x<99>。1、分別用最大似然、最小二乘方法,根據(jù)觀測(cè)結(jié)果求出,和；2、用Monte_Carlo法,計(jì)算出上面兩種方法求出的參數(shù)的偏差和方差；3、利用估計(jì)出的參數(shù),得到目標(biāo)位置的估計(jì),并用Monte_Carlo法計(jì)算在t=0,1,2,…,99s等各個(gè)時(shí)間點(diǎn)上對(duì)目標(biāo)位置估計(jì)的方差和偏差；4、將噪聲的分布改為在〔-1,+1區(qū)間的均勻分布,應(yīng)用最大似然法對(duì)參數(shù)進(jìn)行估計(jì),自己推導(dǎo)該分布下的計(jì)算公式。三、實(shí)驗(yàn)要求：1、設(shè)計(jì)仿真計(jì)算的Matlab程序,給出軟件清單；2、完成實(shí)驗(yàn)報(bào)告,給出實(shí)驗(yàn)結(jié)果,并對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行描述,對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析。四、實(shí)驗(yàn)過(guò)程：假設(shè)s0=0,v0=0.1,a=0.01為實(shí)際值,將其代入,得到。,為白色觀測(cè)噪聲。假設(shè)在t=0,1,2,…,99s時(shí)刻分別取得了100個(gè)觀測(cè)結(jié)果x<0>,x<1>,…,x<99>。分別用最大似然、最小二乘方法,根據(jù)觀測(cè)結(jié)果求出,和；利用最大似然估計(jì)〔ML估計(jì)的目標(biāo)是尋找使得先驗(yàn)概率密度函數(shù)最大的條件作為估計(jì)的結(jié)果,即。這個(gè)數(shù)值可以用下面的公式計(jì)算出：或。即可解出讓似然函數(shù)取得最大值的。本實(shí)驗(yàn)中,其中,分別令,,,,,,,,,可以得到矩陣方程：,則。其中,,〔2線性最小二乘法限定觀測(cè)結(jié)果和待估計(jì)參量之間有下列線性關(guān)系：等式中的H是根據(jù)先驗(yàn)知識(shí)已知的參數(shù)矩陣,n是在觀測(cè)中附加的未知干擾。設(shè)計(jì)目標(biāo)就是尋找一個(gè)使得觀測(cè)誤差平方和最小的參數(shù)矢量作為估計(jì)結(jié)果。2、用Monte_Carlo法,計(jì)算出上面兩種方法求出的參數(shù)的偏差和方差；3、利用估計(jì)出的參數(shù),得到目標(biāo)位置的估計(jì),并用Monte_Carlo法計(jì)算在t=0,1,2,…,99s等各個(gè)時(shí)間點(diǎn)上對(duì)目標(biāo)位置估計(jì)的方差和偏差；4、將噪聲的分布改為在〔-1,+1區(qū)間的均勻分布,應(yīng)用最大似然法對(duì)參數(shù)進(jìn)行估計(jì),自己推導(dǎo)該分布下的計(jì)算公式。將噪聲的分布改為在〔-1,+1區(qū)間分布nt=unifrnd<-1,1>,最大似然估計(jì)的方法要進(jìn)行變化：。對(duì)于均勻分布的情況,無(wú)法寫出聯(lián)合密度函數(shù),計(jì)算最大似然比較困難,用正態(tài)分布結(jié)果進(jìn)行近似。五、實(shí)驗(yàn)結(jié)果：最大似然和最小二乘法得到的結(jié)果如下：可以發(fā)現(xiàn)對(duì)于估測(cè)結(jié)果a、v最好,估測(cè)結(jié)果s0最差。從信號(hào)變化的角度,或許可以這樣理解：隨著時(shí)間變化,信號(hào)發(fā)生變化。其中,提供的加速度和初始速度的信息最多,提供的初始位移信息量最少。最大似然估計(jì)和最小二乘法相比,結(jié)果相差很小,基本符合實(shí)際值。偏差和方差結(jié)果如下：兩種方法得到的結(jié)果基本是無(wú)偏的,并很類似。各個(gè)時(shí)間點(diǎn)上對(duì)目標(biāo)位置估計(jì)的方差和偏差如下：從圖中可以看出,兩種方法下的估計(jì)結(jié)果偏差很小,方差也不大,估計(jì)的效果很不錯(cuò)。4、將程序中estimation函數(shù)nt=randn<1,100>換成nt=unifrnd<-1,1>,重復(fù)〔1〔2〔3容?！?〔3比較可知,噪聲分布改為均勻分布后,參數(shù)估計(jì)的方差和偏差變小,估計(jì)原因是因?yàn)閇-1,1]區(qū)間均勻分布的噪聲方差小于白噪聲。六、程序清單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函數(shù)將當(dāng)前圖形<gca>的x軸坐標(biāo)刻度<xtick>標(biāo)志為:s0,v0,atitle<{'最大似然估計(jì)值'}>;subplot<1,2,2>bar<mean<ls'>>;set<gca,'XTickLabel',{'s0';'v0';'a'}>;title<{'最小二乘法估計(jì)值'}>;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<'最大似然各點(diǎn)偏差'>;subplot<2,2,2>plot<0:99,bias_every_ls>title<'最小二乘法各點(diǎn)偏差'>;subplot<2,2,3>plot<0:99,variance_every_ml>title<'最大似然各點(diǎn)方差'>;subplot<2,2,4>plot<0:99,variance_every_ls>title<'最小二乘法各點(diǎn)方差'>;%estimation函數(shù)function[theta_mltheta_lsbias_every_mlbias_every_lsvariance_every_mlvariance_every_ls]=estimation<s0,v0,a,N>forj=1:N%%%%%%%%%%%%--計(jì)算初始化--%%%%%%%%%%%%%%%%%%%%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>;%'表示求轉(zhuǎn)置,inv表示求逆矩陣,theta_ml<:,j>表示取出第j列%%%%%%%%%%%%--LS計(jì)算--%%%%%%%%%%%%%%%%%%%%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';end