1. 程式人生 > >高斯分佈估計子的效能與克拉默勞下界的討論

高斯分佈估計子的效能與克拉默勞下界的討論

本文演示估計子的定義以及評估方法

% 假設一個N長度的高斯隨機過程x,按照平均數u、方差a分佈,然後根據隨機過程估計他的均
% 值u1,並對評估的有效性進行檢驗

clc
clear 
close all

%% 生成一個隨機變數,並將統計結果與概率密度函式做對比
N = 10;
u = 2.5;
a = 0.8;
xrandom = random('norm',u,a,1,N);

xaxis =(1-a)*u:0.1:(1+a)*u;
fx = normpdf(xaxis,u,a);

figure(1)
subplot(2,1,1)
hist(xrandom,10)
xlim([min(xaxis),max(xaxis)])
title('隨機變數xrandom的統計結果')
subplot(2,1,2)
hold on
plot(xaxis,fx,'LineWidth',4)


%% 引數估計
% 我們根據xrandom過程估計uest(估計子),這個過程可以計作對映T:Xrandom^N -> U,(其中
% Xrandom^N 是N維樣本空間,U是一維引數空間,因為我們只對引數u做估計)
% 顯然,我們隨機過程xrandom(i)中的每一個數據都可以用來生成一個估計子引數uest(i),而且
% 生成的估計子的值uest(i)的值也是隨機的,因此有必要引入誤差B(uest)的概念,用來評估如機
% 子的有效性
% 定義 B(uest) = E{uest - u},如果B(uest) = 0,則認為uest是無偏的的
% 相應的,還有有偏估計 B(uest) ≠ 0,以及漸進無偏估計 limN->∞ B(uest)=0
% 當然這個有偏、無偏的判斷只能通過解析式來判斷,因為計算誤差的存在,數值方法上是無法
% 對其進行準確判斷的。

Euest1 = 1/N*(sum(xrandom,2));
Euest2 = 1/(N+1)*(sum(xrandom,2));
fxest1 = normpdf(xaxis,Euest1,a);
fxest2 = normpdf(xaxis,Euest2,a);
fprintf('u = %d \n E{uest1} = %d \n E{uest2} = %d \n',u,Euest1,Euest2)

%% 對比原始密度分佈函式以及估計出來的密度分佈函式
% 顯然的uest1是u的無偏估計,uest2是u的漸進無偏估計
% 無偏的一定是漸進無偏的,但是漸進無偏的卻不一定是無偏的

plot(xaxis,fxest1)
plot(xaxis,fxest2)
xlim([min(xaxis),max(xaxis)])
legend('原始分佈','估計分佈1','估計分佈2')
title(' xrandom的概率密度函式對比')
hold off


%% 估計子的有效性對比
% 估計子有效性的考察,依賴於估計子代價函式的定義
% 簡單來講,可以通過偏差B(uest)來確定誰更有效
% 但是當兩個估計子都是無偏無極的時候B(uest) = 0,這時候就要考察兩個估計子的方差var(uest)
% 方差更小的被認為更有效
% 當兩個估計子的偏差B(est) 與 aest^2都不相同的時候,就可以引入均方差的概念加以區別
% 定義 M^2(uest) = var(uest) + B^2(uest)

% 問題:如何計算估計子的方差呢?
% 答案:不斷重複隨機過程xrandom,對u不斷的重複進行估計,會得到以系列的uest,這個無偏估計
% uest系列是隨機分佈,且具有方差以及均值的。var(uest)能夠達到的最小值就是克拉默-勞下界
% 此處程式碼不再贅述

關於高斯分佈的Fisher資訊,以及克拉默勞下界是否能夠達到的求解、證明方法。

費舍爾資訊,以及克拉默-勞下界的求解程式碼如下

%% 高斯分佈的費舍爾資訊
% 費舍爾資訊是在被估計引數u已知的情況下,定義出來的一個均值
% 其表示式為 J(u) = E{(x-u)^2/(4*a^4)}
% 其中品質函式 V(x) = (x-u)/(2*a^2)
% 克拉默-勞下界 1/J(u) 是估計子 uest 均方誤差{M^2(uest) = var(uest) + B^2(uest)}能取得的最小值。

Vxrandom = (xrandom - u)/(2*a^2);
JVxsquare = mean(Vxrandom.^2);
fprintf('費舍爾資訊 JVxsquare = %d \n',JVxsquare)
fprintf('克拉默勞下界 1/JVxsquare = %d \n',1/JVxsquare)

以下為程式碼執行結果

u = 2.500000e+00 
 E{uest1} = 2.146888e+00 
 E{uest2} = 1.951716e+00 
費舍爾資訊 JVxsquare = 7.860628e-01 
克拉默勞下界 1/JVxsquare = 1.272163e+00