allenlu2007

A great WordPress.com site

Category: math

Forney Style Factor Graph

by allenlu2007

Which Graph to Use

以下的四種 graphs 都可用來表示 p(u, w, x, y, z), 但並非等價。可以找到一些例子說明。

例如 Baysian network (directed graph) 的 V structure 無直接對應的 undirected graph.  Undirected graph 的 loop structure 也沒有對應的 (acyclic) directed graph.

即使 undirected graphic 的 factor graph 和 markov random field (MRF) 也並非完全等價。Factor graph 以 factor 為單位,MRF 似乎以 clique 為單位 (?).   一般認為 factor graph 比較嚴謹且 general,但是太繁瑣因為同時要畫 variable and factor node (如下圖 original factor graph).   Forney-style factor graph 用 arc 來取代 variable node, 只留下 factor node, 看來比較精簡。不過需要解決 shared variable (下圖的 X 連接兩個以上 factors) 的情況,因而引入 “=” 的 node.  

一般有清楚因果關係的大多用 Baysian network.  各 variable 平等的應用大多用 MRF, 如圖像處理。Factor graph 在 error correction code 中常常使用。以下使用 Forney style factor graph.

 

NewImage

 

Message Passing in Factor Graph

Message passing MP (or belief propagation BP) 分為兩類:  sum product and max product.

在 factor graph 的 MP for sum product 如下圖所示。基本上以 variable 為單位

 

NewImage

 

or with message passing:

NewImage

可以 reuse message:

NewImage

 

 

In general sum product and max product 可以用以下的 summary message.  針對每個相連的 factor 可以得到 summary message ( summation or maximisation of products of factor and all other incoming messages except for the outgoing one.  最後的結果是 summary message 的 products.

 

Step 1: compute all messages (starting at the leaves!)

NewImage

 NewImage

在 loop-free 的情況下,sum-product algorithm:NewImage

Message Passing in MRF

 

MRF 的 message passing 兩個例子

1: Markov Chain (沒有 shared variable)

NewImage 

Node: variable X1, .. X5.  

Factor: arc (每一個 factor 只有兩個 variable, no shared variable)

剛好和 Forney style factor graph 相反。

P(X2, X5) = P(X5) P(X2|X5) 

P(X2|X5) :  X5 是 evident node; X2 是 unknown node; 其他是 hidden node

所以 X2 看到的 message 為左右的 product.  Format 和上文一樣: 

sum of products:  sig{ P(X3|X2) * delta(4) } and sig { P(X2|X1) * delta(1)}

 

2.  Clique tree example (簡化後和以上相同)

NewImage

 

 

Message Passing in WiKi

NewImage

NewImage

Wiki 的說法和之前的公式一致,只是分為兩部份。 messages: (variable v to factor u) and (factor u to variable v)

 

General Formula

Find P(X | Y) where X and Y are vectors

X: x1, x2, x3… are unknows

Y: y1, y2, y3… are evidences

Z: z1, z2, z3… are hiddens

 

 

 

 

 

 

 

Hidden Markov Model of Filtering, Smoothing – Probabilistic Interpretation

by allenlu2007

前文 Kalman filter 是從 (linear) state space model 出發:

NewImage

 

 

我們可以從 Hidden Markov Model (HMM) 出發,更 general 的 probabilistic model:

Sometimes it is called Bayes filter or recursive Bayesian filter.

Kalman filter 是一個 linear space + Gaussian distribution 的特例。

好處:

1. 更 general

2. 更好的 nonlinear filter or large dimensional filter: Particle filter, Unscented KF, EnKF

 

The following is reference from xxx

NewImage

 

上述的 joint PDF 其實是 Bayesian Network (BN) 的直接結果。Joint PDF 等於所有的 conditional PDF product!

 

The key to solve the estimation problem is divided into filtering and smoothing!! 

The first equation is directly from Markov property (since t1:n is given, there is only one proportional parameter!)

NewImage

同樣的結論 from another article:

NewImage 

(1) state is X;  output is e

(1a) combine prediction and update into one equation

(2) 多了 smoothing (backward pass).  Smoothing 只有 one step; 不像 filter: prediction+update

(3) summation vs. integration

 

 

如何把上述 equation 轉為常見的 Kalman filter and smoother?

 

Filtering

Filtering 分為兩步: prediction and update (or innovation)

Prediction:

(y is state, t is observed data)

NewImage

or ( x is state, z is observed data)

NewImage

or ( x is state, e is observed data)

NewImage

 

Update (innovation)

NewImage

or more precisely

NewImage

where alpha > 1

 

 

Case 1: Linear model : Kalman filter

NewImage

After prediction  …

After update:

NewImage

一個重點是 Kt and Sigma_t 和 observation data 無關,可以 offline compute.  問題是如果 Sigma_x, Sigma_z, 和初始的 Sigma_t 如果不準,所算出的 Kt 自然也不準。可能需要用 EM algorithm 事先 training. 

Or use EnKF to compute the K!

 

Case 2: Nonlinear model:  Extended Kalman filter (1st order approximation)

如果 nonlinear (therefore non-Gaussian) but well define function with reasonable dimension.

We can use extended Kalman filter.  It is essentially a linear approximation of the nonlinear function.

 NewImage

接下來就和 KF 是一樣的。

 

Case 3: Nonlinear model: Unscented Kalman filter (2nd order approximation)

NewImage

比較 Taylor expansion and statistical linearisation 對 input Gaussian PDF to output 的差異如下圖

NewImage

可以看到 statistical linearisation 所產生的 output PDF 的 1st and 2nd order (mean and variance) 都比 Taylor expansion linearisation 的近似好。這是 unscented KF 的基本概念。

NewImage

NewImage

D 是 Gaussian 的 dimension.  1-D 3 points; 2-D 5 points, etc.

NewImage

Sigma point KF 的好處:

* 不需要計算 derivative

* Based on deterministic sampling (非 Monte Carlo)

* Gaussian approximation filter with exact nonlinear models (非 Taylor expansion approximation).  

缺點: 收歛速度比較慢

 

NewImage

 

 

Case 4: particle filter (or sequential Monte Carlo filter)

For nonlinear and non-Gaussian and high dimension case.  There is no well defined nonlinear function.

NewImage

 

Case 5: Ensemble Kalman filter (EnKF)

A special case of particle filter but everything is Gaussian.  The only problem is high dimensional so it is difficult to compute the covariance matrix.

The EnKF is a monte carlo approximation of the Kalman filter, which avoids evolving the covariance matrix of the state vector x.   K = CH’ (HCH’ + R ) ^ -1

Recall in Kalman filter, Kt and Sigma_t are independent of the observed data (only depends on Q and R)!!  In the EnKF, this is not the case, but it is better because it takes the observed data into account? and avoid the inaccuracy prior assumption of the state covariance?

 

 

 

 

 

Q: what’s the link of HMM and conventional Markov chain?

   (a) conventional Markov chain (如下), 沒有 show time step, only the state transition and transition probability.

       BTW, the state transition and probability is time invariant.

    (b) conventional Markov chain, the state is directly observable.

    © HMM 的圖如最上。each state circle includes all possible states; and shows the time step transition.

 

    Answer:  Daphne Kroller 的解釋很清楚 (PGM in coursera).  Conventional Markov chain 定義了 state transition model and transition probability, 而 hidden markov model 代表的是 temporal unroll 的圖。 Markov property and time invariant 是一般 HMM 的 default 假設。  當然 HMM, state 是 unobservable.  

NewImage

NewImage

Kalman Filter, Kalman Smoother, and Expectation Maximization

by allenlu2007

 

前文提到 Kalman filter 的好處 (i) 探索 data 的內在結構 (state equation and noise covariance matrix), 以及 (ii) 發展出 recursive equation 來 estimate 內在的結構 (state or sometimes noise).

實際上遇到的困難是 (a) 這樣的 state equation 是否存在? (b) 如果存在如何找到正確的 model 參數? and (c) 如果 model or 參數不準確,Kalman filter estimate 的結果是否 robust?

(a) State equation 是否存在

有兩種解釋: (i) state equation bases on physics  如力學、電子學、etc. (ii) dynamic system 的線性近似。

(b) 如何找到 model 的參數

(i) 由 equation 推導出. (ii) 是由 training 出來的。如果是 (ii) 一般不會差太多。

(c) 如果 model or 參數不準確,Kalman filter estimate 的結果是否 robust?

Again, 可能要用 training 方式來找 Kalman filter 的 Q, R, A, C.  避免 Kalman filter output 和現實差太多。



Kalman filter Training (parameter estimation) 的方式可以用 EM algorithm 或其他 system identification algorithm (ARMA, .. etc.)

Kalman filter is typically used for on-line state estimation and a minimum-variance smoother may be employed for off-line or batch state estimation. However, these minimum-variance solutions require estimates of the state-space model parameters. EM algorithms can be used for solving joint state and parameter estimation problems.

1. 似乎是 Kalman smoother 才可以。Kalman filter 不能?  是因為必須 batch mode?


Step 0: Kalman filter and Kalman smoother 的差異如下:

(a) Kalman filter 主要是”當下”的最佳 estimation.  依照 Markov state 的定義,只要有前一 state 的 pdf/information (mean and covariance matrix; mean 是 x_hat; covariance matrix 是 P), 就可以 ignore 更早的 information.  

(b) Kalman smoother 則是用未來的 observed data 來回推最佳的 estimation.  By the same token, 似乎只要”下一個” sample 就足夠。實則不然,因為需要的是下一 state 的 pdf/information, 並非只是 observed data.  下一 state 的 pdf/information 可以被 estimate 更精確如果知道下下 sample (because they are correlated), 依此類推。當然也可以只看下一個 (或下幾個) sample 來做 Kalman smoother, 稱為 fixed-length or fixed-interval smoothing.  只是所得到的 variance reduction 不是最佳。實務上很少在 real time application 用 Kalman smoother, 主要是 complexity 以及額外的 latency for smoothing.

NewImage

Step 1: Kalman smoother State Estimation:

因此完整的 Kalman smoother 需要所有的 data 來做 post-processing (Viterbi 例外?)

 

NewImage

對 Kalman filter 來說,xt 是一個 Gaussian random variable.  Mean 是 x_hat; covariance matrix 是 P.  P 在 prediction 時增加 (if Q>0),在 update 時減小。

Kalman smoother 多了 backward pass 可以給更好的 estimation, covariance matrix 應該要減小.  此時並非 real time.  理論上 Pt|T < Pt|t (how to prove?)

Find an example to show the difference of Kalman filter and Kalman smoother!

 

 

Step 3: EM algorithm of State Estimation

綜上所述 Kalman smoother 的 state estimation (mean and covariance) 優於 Kalman filter.  Trade-off 是 complexity and latency.  對於 real-time applications 如 tracking, navigation, etc., Kalman filter 比較適合。但若是一些非 real-time 的影像或聲音處理,Kalman smoother 可以多得到一些 performance.  

Kalman smoother 其實還有一個重要的應用,就是配合 EM algorithm, 用來 estimate A, C, Q, R (當然同時 estimate state 的 mean and covariance)! Magic!

因為在實務上,Q and R 大多非 prior information (如何事先知道 state noise and measurement noise 的大小?) 因此可以靠 training (with or without prior state information).  更進一步,有時 A and C 也是 unknown, 也需要由 training 而來。

結合 EM algorithm 和 Kalman smoother 可以做為 training (without prior state information!!) 而得到 A, C, Q, R 的 estimation.  再 performance real-time Kalman filter.   

EM 的 E-step 就是 Kalman smoother.  M-step 就是 update A, C, Q, R based on ML estimation.  EM algorithm 需要解決兩個問題: (i) 如何決定 A, C, Q, R 的 initial guess.  (ii) 如何避免 local minimum.

NewImage


Example: 


 

 

 

 

More EM Examples

by allenlu2007

Based on Andrew Ng’s lecture, there are more applications for EM algorithm.

 

1.  GMM: done in the previous article

2.  Naive Bayesian Mixture:  Two types:   multivariate binomial mixture;  multinomial mixture: commonly used in machine learning.

3.  Factor analysis: to reduce the dimension of observed data  (http://newgenerationresearcher.blogspot.tw/2009/02/factor-analysis.html)

 

以及其他的包含

* Kalman filter 的 A, C, Q, R estimation.

* Hidden markov model (HMM) 的 estimation? Baum Welch algorithm  (和以上應該是同一件事?)

* Possion random process: wavefront sensing

 

 

Machine Learning Lecture

by allenlu2007

這幾天把 Stanford CS229 Andrew Ng 的 machine learning lecture 看了一半,有幾點感想。

1. Ng 的理論背景扎實,lecture handout 都很清楚。每個証明的策略和細節他都一一走過,非常不容易。至少我做不到。

2. 不過有些部份以乎還不夠深入。例如牛頓法在最小值為什麼收歛比 gradient descent 快。object function with regulation 的 trade-off. Boyd 的 convex optimisation lecture 說的更深入。兩者有互補作用。

2a. supervise learning 的 ML 的 parameter estimation 是 maximize joint likelihood function 而非 conditional likelihood functions.  這點很重要!!  結論是 ML parameter estimation 就是要 maximise 所有 observed data 的 probability.  In supervise training 的 training phase,  both input and output are observed, 所以當然是 maximise joint likelihood.  但在一般的 parameter estimation, internal state or unsupervised learning case, 只有 output 是 observed, 所以 maximise conditional likelihood.

3. Ng 解釋很好的部份有: logistic regression, support vector machine,  machine learning debugging (high bias or high variance, algorithm convergence problem or objective function problem, etc.)  特別他解釋 expectation maximisation 部份讓我之前對於整個 logic 更清楚的了解。之後再整理這部份。

4. 交互比較 Ng and Boyd 的 lecture and note 應該可以收獲更多。 

Convex Optimization

by allenlu2007

 

Convex optimisation 的問題

找到最小值

NewImage

 

NewImage

 

通常上述的 equations 沒有 close form, 需要用迭代法逐步找到最小值

NewImage 

事實上應該是 a * -delta f(x).  a is a step parameter.  

這個方法稱為 gradient descent method 或 steepest descent method.  梯度下降法或最陡下降法。理論上如果 a 很小,就是1st order 近似法的最佳解。但是收歛的時間會很久。如果 a 太大,在一些 local ellipsoid 很扁平的情況下會 overshoot, 同樣會有收歛慢的情況。

2nd order 近似法的最佳解就是牛頓法。最大的缺點是 Hessian 可能很難算。

從另一角度來看,也可以說 steepest descent method 是基於 unit circle (Euclidean norm) 的最佳解。但實際情況需要用 quadratic norm (ellipsoid).  另外也有 l1 norm.

這部份可以參考 Boyd 的 notes.

NewImage

Steepest descent   –>  Euclidean norm  –> 1st order approximation

NewImage

 

牛頓法原理

原始的牛頓法是用來解 ( f(x) = 0 )

NewImage

如上圖所示,可以得到:

NewImage 

以收歛速度而言,只是一階的近似。

其實牛頓法最 powerful 是在找最大或最小值。因為 f'(x) = 0.

x(n+1) = x(n) – f'(x(n)) / f”(x(n))   

以上圖求最小值來說,一階導數是直線。一次就可以找到最小值。

因此牛頓法找最大或最小值,是二階近似。

Two dimension or higher dimension, use Hessian to replace f” and gradient to replace f’.


 

 

 

 

 

Hidden Markov Chain – HMM

by allenlu2007

 

和 markov chain 的差異是非直接觀察 markov chain.  而是另外的 observables.

HMM  有各種不同的應用 Applications:

communication: convolution code, trellis code modulation

speech recognition

Kalman filter?

 

How to solve HMM?

1. Viterbi algorithm

2. Expectation Maximisation

 

Expectation Maximization EM Examples

by allenlu2007

 

 

 

 

本文主要參考  

 

Reference1 是 Stanford 的 tutorial paper.  

 

Reference2 是 Washington 的 tutorial paper. 

 

 

EM algorithm 主要用 iterative method 來計算 ML 或 MAP parameter estimate. EM 一個主要的好處是用在 missing data or hidden data or latent variables. 

不過至少對我而言,不容易直覺了解 EM algorithm.  Too good to be true?  以下是幾個 examples 

Missing Data Estimation

Example1a: Normal distribution with unknown mean (but known variance)

NewImage  

observables

NewImage  

NewImage 

equals 

NewImage

NewImage

 

For missing data

NewImage

NewImage

 

x is complete data, y is observables, z is missing or hidden data

x = [y, z]   where y = (x1, .. , xm), z = (xm+1… xn)

NewImage 

 

NewImage

 NewImage

 

依照前文的現在使用 EM algorithm, 第一步是 Expectation step:

 NewImage

NewImage

推導的過程如下

NewImage

when t–> infinity   u = 1/m (sum 1 to m Xi), same as ML estimation without missing data!

 

Example1b: Normal distribution with unknown mean and variance

其實 1a 是 1b 的特例。比較完整的推導如下

NewImage 

當然 1a 和 1b 是一個 trivial case, 因為 mean and variance 最後收歛到 ML estimator without missing data.  引入 missing data 無法榨出更多的 information, 也沒有什麼意義。因為  missing data (x(m+1) … x(n)) 和 observables (x(1) .. x(m)) 之間是 independent.  In general, can we prove missing data of i.i.d. always converge to MLE without missing data?

下面的例子是另一類的 hidden information.


Example2: Fisher’s Multinomial Distribution

The history is cited here.

NewImage

pmf  (y1, y2, y3, y4) are integers.

NewImage

where y = (y1, y2, y3, y4).    = (125, 18, 20, 34)

NewImage

 

Let  the complete data x = (y11, y12, y2, y3, y4) where y11+y12 = y1 

NewImage

 NewImage

 where Bin = Binomial  n = y1  p = 1/2 / (1/2 + phi/4)!   E(Y11) = np as follows:

   NewImage

y12(0) = y1 – y11(0).   This complete the E part.

 

E(log Lc | y ) = (y12(0)) + y4 )  log (phi) + (y2 + y3) log (1-phi)

Set derivative to 0  ==>  phi(1) = (y12(0) + y4 )  / (y12(0) + y2 + y3 + y4) = (y12(0) + y4 ) /  (n- y11(0))

ans = 0.6268214


Mixture Estimation 混合估計

混合估計在統計上有實際的應用如 pattern recognition, data mining, etc.  甚至數位通信都可視為混合估計.

混合估計的問題可以分為不同的困難度: (i) 已知有幾個 pdf 以及參數,估計未知的權重; (ii) 已知有幾個 partially specified pdf, 估計未知的參數和權重; (iii) 最困難的是未知幾個 partially specified pdf, 估計未知的參數和權重。

本文 (i) and (ii) by applying EM and Gibbs sampling schemes respectively. We will address task (iii) as an illustration of reversible jump MCMC later.

我們考慮幾種 distribution, Binomial, Poisson, and Gaussian.  

Binomial 和 Poisson mixtures 是 discrete distributions.  Poisson 常常用在光的 sensor 應用上如 communication or image. 

Gaussian 是 continuous distribution mixture, 有廣泛的應用如 data mining, big data, etc.  

Assume that a sample X1, X2, .. Xn comes from the mixture

NewImage

where the sum of wi equals to 1, and wi > = 0.   fi(x) could be Binomial, Poisson, or Gaussian.

 

Case 1: 

先考慮最簡單的 case, only estimate the weights, the log-likelihood function

NewImage

單純的 ML estimator 計算對 wi 的導數,沒有 closed form. 

如果用 EM 的 latent variable 方式,除了 (x1, x2, .., xn) observables 之外,引入 unobservables matrix Z = (zij, i and j from 1.. n).  zij = 1, observation xi comes from distribution fj;  0, else.

Therefore, the new likelihood and log-likelihood of Lc(w) is shown below:

NewImage

NewImage
E-Step

NewImage  

  NewImage

 

  The result is regardless of the distribution. 

 NewImage

   

Example3:  Gaussian Mixture

 

Case 2:  unknown mean and variance

First a common variance

NewImage

Independent variance

NewImage

 NewImage

NewImage

 NewImage

NewImage

 

  

  


Special case of two Gaussians

NewImage

NewImage

Given hidden variable sigma_i = (0, 1),

 

NewImage 

 

E – step:
NewImage  
NewImage
NewImage

A matlab code is listed below:

MU1 = [1 2];

SIGMA1 = [2 0; 0 .5];

MU2 = [-3 -5];

SIGMA2 = [1 0; 0 1];

X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];

 

scatter(X(:,1),X(:,2),10,‘.’)

hold on

 

options = statset(‘Display’,‘final’);

obj = gmdistribution.fit(X,2,‘Options’,options);

 

h = ezcontour(@(x,y)pdf(obj,[x y]),[-8 6],[-8 6]);

 

ComponentMeans = obj.mu

ComponentCovariances = obj.Sigma

MixtureProportions = obj.PComponents

 

  
 NewImage


 

 

 

  

 

  

 

 

 

 

 

 

 

 

 

 

常見估計 Estimation 的簡結

by allenlu2007

Least square estimation (LSE), maximum likelihood estimation (MLE), minimum mean square error estimation (MMSE), maximum a posterior esimation (MAP), minimum mean absolute value of error estimation (MAVE), 相當多的 estimation 方式。差異何在,以及如何歸類?

最常見也基本的是 least square estimation (LSE), 常用在 data fitting 或是 linear regression, 主要是 minimize sum of squares error.  對於 observed data or experimental data 沒有假設分佈或是假設高斯分佈.

Maximum likelihood estimation (MLE) 常用在 parameter estimation and inference, 基本上需要 oberved data 的 pdf 分佈.  但是對 estimate 的參數 ( 如 mean, variance, frequency, phase, etc.) 沒有任何假設或 prior information, 屬於 classical estimation.  對於沒有任何 prior information 的情況,MLE 一般是最好的 estimation.

在高斯分佈的特例中,MLE 和 LSE 會得到一樣的結果。

MMSE, MAP, and MAVE 則是屬於 Bayesian estimation, 基本上需要 oberved data 的 pdf 分佈.  同時對 estimate 的參數 ( 如 mean, variance, frequency, phase, etc.) 也有 prior information (Bayesian estimation).  MMSE, MAVE, and MAP 主要的差別是在 cost function: MMSE (mean square error), MAVE (mean absolute value of error), MAP (hit or miss).

Parameter Estimation 參數估計

by allenlu2007

What is Parameter Estimation

參數估計是指用樣本指標(稱為統計量)估計總體指標(稱為參數)。例如用樣本 mean 估計總體 mean 以及用樣本 variance 估計總體 variance。

Parameter Estimation 分類

Parameter estimation 有點估計 (point estimation) 和區間估計 (interval estimation) 兩種。

點估計是依據樣本估計總體分佈中所含的未知參數或未知參數的函數。通常它們是總體的某個特征值,如數學期望、方差相關係數等。點估計問題就是要構造一個只依賴於樣本的量,作為未知參數或未知參數的函數的估計值。例如,設一批產品的廢品率為θ。為估計θ,從這批產品中隨機地抽出N個作檢查,以X記其中的廢品個數,用X/N估計θ,這就是一個點估計。

點估計常用的方法是:

  (1) Classic estimation:   特征值參數θ是一個定值 (non-random), 同時沒有任何先驗資訊 (a priori information).  例如 Maximum likelihood 最大似然估計法。於1912年由英國統計學家 R.A. Fisher 提出,利用樣本分佈密度構造似然函數來求出參數的最大似然估計 (Maximum likelihood estimation)。

  (2) Bayesian estimation: 基於貝依氏學派的觀點而提出的估計法。特征值參數θ是一個 random variable, 同時允許 a priori PDF on θ.  例如 Maximum a posteriori (MAP) estimation, 或是 Minimum mean square error (MMSE) estimation.

區間估計是依據抽取的樣本,根據一定的正確度與精確度的要求,構造出適當的區間,作為總體分佈的未知參數或參數的函數的真值所在範圍的估計。例如人們常說的有百分之多少的把握保證某值在某個範圍內,即是區間估計的最簡單的應用。1934年統計學家 J. Neyman 創立了一種嚴格的區間估計理論。本文只討論點估計的 parameter estimation.

Bayesian Parameter Estimation

Given the conditional pdf

p({y}_{1}, ..., {y}_{N}|x) = p(\mathbf{y}|x)

以上 {y}_{1}, ..., {y}_{N} 是實驗值。以及 x 是被估計的參數 (unknown),根據貝依氏的理論是一個 random variable.

再來定義 \hat{x} 是 x 的 estimator, 以及估計的誤差

\epsilon_x = x - \hat{x}(\mathbf{y})

Cost Function and Bayes Risk

貝依氏參數估計的重點就是引入一個 cost function c(\epsilon_x), 或 loss function.   注意 cost function 仍是一個 random variable.

我們稱 cost function 的期望值為 Bayes risk, R relative to the joint pdf p(\mathbf{y}, x) .

R = E[c(x-\hat{x}(\mathbf{y}))] = \int \int c(x-\hat{x}(\mathbf{y})) p(x,\mathbf{y}) dx d\mathbf{y}

貝依氏的參數估計就是 minimize Bayes risk R.

以下舉幾個 cost function 的例子以及 Bayes estimation 如何 link 到常見的 estimator.

Quadratic cost function (如下圖一):

{c}_{q}({\epsilon}_{x}) = {{\epsilon}_{x}}^{2}

Uniform (binary) cost function (如下圖二):

{c}_{u}({\epsilon }_{x}) = \begin{cases} &0 \text{ if } |{\epsilon }_{x}| \leq \frac{\delta }{2} \\ &1 \text{ if } |{\epsilon }_{x}| > \frac{\delta }{2} \end{cases}

NewImage

Minimum Mean Square Error (MMSE) Estimator

MMSE estimator 對應為 quadratic cost function.  Bayes risk R:

R_q = E[c(x-\hat{x}(\mathbf{y}))] = \int \int (x-\hat{x}(\mathbf{y}))^2 p(x,\mathbf{y}) dx d\mathbf{y}

可以用 Bayes formula 簡化上述公式

R_q = \int d\mathbf{y} p(\mathbf{y}) \int (x-\hat{x}(\mathbf{y}))^2 p(x|\mathbf{y}) dx

Since both integrals are not negative, minimize risk R_q 相當於 minimize

\int (x-\hat{x}(\mathbf{y}))^2 p(x|\mathbf{y}) dx

可以証明

\hat{x}_{MMSE} = \int x p(x|\mathbf{y}) dx = E[x|\mathbf{y}] 


Maximum A Posteriori (MAP) Estimator

 NewImage

To minimize Ru 相當於 maximize the following integral.

NewImage

\hat{x}_{MAP} = \arg \max_{x} \log p(x|\mathbf{y})

可以觀察到  MMSE 和 MAP estimators 都用到 posterior pdf p(x|\mathbf{y}).

MMSE 是 center of mass, MAP 是 mode (peak) of pdf.  如果 conditional pdf is symmetric posterior pdf the MMSE and MAP estimators are equal.


What’s the difference between Cramer Rao bound and MMSE estimator??

NewImage

Maximum Likelihood (ML) Estimator

考慮一個特例, x 是未知但定值 (unknown but deterministic). 其實是 classic estimation, 是否能用以上 Bayesian estimation 解釋?

NewImage

NewImage

\hat{x}_{MAP} = \arg \max_{x} \log p(x|\mathbf{y})

可以簡化為

\hat{x}_{MAP} = \arg \max_{x} \log p(\mathbf{y} | x)

這是 Maximum likelihood estimator. 

當然從 classic estimation 角度,出發點完全不同。通常我們先把 x (random variable) 換成 \theta (non-random variable).  再定義 likelihood function L, ML estimator 即是找 Likelihood function maximum.   或是 Log-likelihood function maximum.

L(\theta | \mathbf{y}) = p(\mathbf{y} | \theta )

\frac{\partial \log L(\theta | \mathbf{y})}{\partial \theta } = 0

\hat{x}_{ML} = \arg \max_{x} \log p(\mathbf{y}|x)

所以 MAP estimator 根據 Bayes formula 可以改寫為

\hat{x}_{MAP} = \arg \max_{x} \log p(\mathbf{y}|x) p(x)

Remark

如果 p(x) is uniform distribution (此為充分條件,非必要條件), MAP 和 ML estimators 相同。

 NewImage

其實這正反應了 ML estimator 完全 no prior information of x, 所以才會是 uniform distribution (不管 x 在那的機會均等).  如果我們事先知道 x 的 pdf (非 uniform distribution), 就變成 MAP estimator, 理論上可以最更好的 estimation (假設 prior information 是準確的)。  這也是 Bayesian estimator 的精神所在。


 

Maximum Mean Absolute Value of Error (MAVE) Estimator

NewImage

 

 NewImage

From the above equation MAVE estimator is the median of the posterior density.