# 補充：Bayesian Estimation

我們直接從一個例子開始，假設我們要從一個 Poisson distribution 中挑出一個 observation，其中這個 distribution 的 parameter（同時也是 mean $\lambda$）的值，是 $2$ 或 $4$。
> 關於 Poisson distribution 的介紹，有興趣可參考筆記「[補充：Poisson Distribution](https://hackmd.io/@pipibear/SkrUqjGIA)」，但在此不熟悉這部分也不影響。

除此之外，假設在執行 experiment 之前，我們就相信 $\lambda = 2$ 的機率是 $\lambda = 4$ 的四倍，也就是說：

\begin{equation}
\begin{split}
P(\lambda = 2) = 0.8 \\
P(\lambda = 4) = 0.2
\end{split}
\end{equation}

接著，我們就真的去執行這個 experiment，並得到結果 $x=6$。

根據查表可以算出來的資訊，我們得到下方結果：


如果按照直覺，我們會覺得：

即使這整個 distribution 中，$\lambda =2$ 的機率是 $0.8$，但實驗解果出來 $\lambda=4$ 的情況下更有可能得到結果 $6$（$\because 0.104 > 0.012$），那或許 $\lambda$ 其實更有可能是 $4$。

為了實際驗證看看，我們就會去計算 posterior density of $\lambda$，意思是「在得到 outcome $=6$ 的情況下，$\lambda=2$ 和 $\lambda=4$ 的機率各自是多少。」

方法如同一般的 Bayes' 計算 posterior probability 的方式，過程如下：


從這樣的結果我們發現，在做出 observation $x=6$ 以後，$\lambda=2$ 的 probability 從 $0.8$（prior probability）降到 $0.316$（posterior probability）。

## prior pdf

雖然上面的例子中，我們的 parameter 只有 $2$ 和 $4$ 兩種可能，但在現實應用中，我們的 parameter （以 $\theta$ 表示）可能會有非常多種可能的值。

因此，我們可以對 parameter space 中的每一個可能的 $\theta$ assign 一個 prior probability。
> - parameter space：所有可能的 parameter 值所成的集合。
>
> $\rightarrow$ 意思就是假設 $\theta$ 可能的值是 $\{1,2,...,10\}$，那麼我們就 assign $p(1)=$某個值、⋯⋯、$p(10)=$某個值，在真正做實驗前先假設 parameter 分佈在各種可能的值的機率。
>> 至於要怎麼在真正做實驗前就知道這些值，那是統計學家要處理的事，我們只要使用就好。

既然要去表示 $\theta$ 分佈在各個值的機率，那麼也就類似於我們之前描述 random variable 各個 outcome 的 probability 時，用 pmf / pdf 來表示一樣：

$\rightarrow$ 我們用 <font color = "snake">prior pdf</font> ==$h(\theta)$== 來表示 $\theta$ 在各個可能的值的 probability distribution。
> 舉例來說，我們有一個 $\theta_1 = 1$，將它作為 prior pdf 的 input 後，$h()$ 這個 function 就會給出 parameter $\theta$ 真的是 $\theta_1$ 的機率，如： $h(\theta_1) = 0.1 = P(\theta = \theta_1)$。

### noninformative prior

如果 $h(\theta)$ 是 constant，則 $\theta$ 具 uniform distribution，那我們就說這是 <font color = "snake">noninformative prior</font>。
> 例如 $h(\theta_i)=0.1 \quad \forall i \in \{1,2,...,10\}$
> 代表我們的 $\theta$ 有十種可能值的情況下，每一種的機率都是 $0.1$。
>
> $\rightarrow$ 說每一種可能都一樣有機會，這樣的資訊基本上沒什麼幫助，所以才說是 "noninformative"。因此，如果 $\theta$ 有什麼資訊是我們在做實驗前就知道的，就要盡量避免 noninformative prior。

---

除了 $\theta$ 可以不只是兩個值，我們在實際應用中通常也不會只做一次 experiment，只用一個 observation，我們會取++許多的 observations++，也就是一個 ++random sample++。

有了 random sample 中那麼多的 observations 後，我們通常就能計算出一個還不錯的 ==$Y$==，代表 ++statistic for parameter++ $\theta$。
> 任何從 sample 計算出的值都可以稱為 statistic，在此的 $Y$ 這個 statistic 是由 sample 的資訊，經過計算取得的估計 $\theta$ 的值。
>> 關於 statistic 的介紹可回顧筆記「[4.1 Introduction](https://hackmd.io/@pipibear/S1oHBgxSA)」。

假設我們現在討論的 statistic $Y$ 是 continuous 的，那 $Y$ 也會有它的 distribution，因此也有 ++$Y$ 的 pdf++，我們用 ==$g(y;\theta)$== 表示，代表 ++given parameter 是 $\theta$ 的條件下，$Y$ 的 pdf++。
> 用一個例子回顧整個脈絡：
>
> 假設我們希望預測某個地區在七月的降雨量，根據一些知識背景或觀察，我們覺得這件事會符合某種特定的 distribution，也就是我們知道了預測的 function 大概長什麼樣子，但是在這個 function 之中有一個或多個 parameter 的值我們不清楚，於是我們找了一些過去幾年的資料來幫助我們找到這個（或這些）未知變數的值。
>
> 有了這些資訊之後，我們就用這些 observations 去做一些計算，進而猜測我們未知的 parameter $\theta$ 實際上是什麼（假設只有一個未知的 parameter。）
>
> 舉例來說，假設真正的 $\theta = 5$，而我們算出來的猜測值可能是 $6, \ 5.5, \ 4$⋯⋯，這些值也會各自有自己的機率，好比說我們取十年的資料，算出 $6$ 三次，也就是 $6$ 的機率是 $0.3$、算出 $5.5$ 四次，也就是是 $5.5$ 的機率是 $0.4$⋯⋯。
>
> 因此，這些猜測的值也會有一個 distribution，我們就把描繪這個 distribution 的 function 訂為 $g(y;\theta)$。
>
> 在這個例子中，$g(Y = 6;\theta) = 0.3$

- $g(y|\theta) = g(y;\theta)$ 是通用的。

## joint pdf of statistic and parameter

接著，我們可以將下方式子視為 ++statistic $Y$ 和 parameter $\theta$ 的 joint pdf++：

:::info
\begin{equation}
g(y|\theta)h(\theta) = k(y,\theta)
\end{equation}
:::
> 我們把 $k(y,\theta)$ 令為：
>
> parameter 是 $\theta$ 的機率（$h(\theta)$），乘上 parameter 是 $\theta$ 的情況下， statistic $Y=y$ 的機率。

$\rightarrow k(y,\theta)$ 也就是 statistic $Y$ 和 parameter $\theta$ 的 joint pdf。
> 關於 joint pdf 的介紹可參考筆記「[補充： Joint distribution functions](https://hackmd.io/@pipibear/rkrvIg2I0)」。

## marginal pdf of statistic

上面這個式子，也可以用來轉換成表示 ++marginal pdf of $Y$++：
> Recall： marginal pdf of $Y$ 的意義是，不管 $\theta$ 是什麼值，$Y$ 自己的 distribution。

:::info
\begin{equation}
k_1(y) = \int_{-\infty}^{\infty}g(y|\theta)h(\theta) \,d\theta
\end{equation}
:::
> 正因為是 $Y$ 自己的 distribution，所以我們將 joint pdf $k(y,\theta)$ sum over all possible $\theta$。

## posterior pdf

接著，如果我們把這兩個東西相除，再利用 Bayes' Rule，我們就可以得到下方的結果：


我們把 ==$k(\theta|y)$== 稱作 <font color = "snake">posterior pdf of $\theta$</font> (given that $Y=y$)。

對自稱 Bayesians 的統計學家們來說，他們相信所有關於我們未知的 parameter $\theta$ 的資訊都可以被總結在這個 posterior pdf $k(\theta|y)$ 中。

舉例來說：

:::warning
在對 parameter $\theta$ 做 point estimate 時，就等同於我們在 pdf 為 $k(\theta|y)$ 的情況下，去猜測 random variable $\theta$ 的值。
:::
> Recall：point estimate 的意思是，我們從 parameter space $\Omega$（包含所有 $\theta$ 可能的值所成的集合）中挑出某個 $\theta$ 值。

## estimating parameter：決定 w(y)

那要怎麼透過這個 pdf 去猜 $\theta$ 的值呢？

$\rightarrow$ 有很多方式，像是取平均、求中位數，或是看這個 distribution 是什麼樣的模式⋯⋯，但是其實最好的方式是：

:::warning
根據猜錯以後所產生的 penalties，來決定 best guess。
:::

舉例來說：

假設我們令++猜測的值++為 ==$w(y)$==，parameter 實際的值為 $\theta$。

猜錯的 penalty 訂為猜測的值和真正的值之間的差距平方，也就是 $(\theta - w(y))^2$，那我們就會發現，其實最好的（penalty 最小的）情況會發生在 $w(y)$ 訂為如下的 conditional mean：

\begin{equation}
w(y) = \int_{-\infty}^{\infty}\theta k(\theta \ | \ y) \,d\theta
\end{equation}

為什麼呢？

假設我們現在有一個 random variable $Z$，我們希望 mean square error $E[(Z-b)^2]$ 這個值越小越好，那麼其實最小的情況會發生在當 $b = E(Z)$ 時。
> - 關於 mean square error 可參考筆記「[4.3 Evaluating an Estimator: Bias and Variance](https://hackmd.io/@pipibear/BkS_CMuS0)」。
>
> 這裡我就不再證明為什麼當 $b = E(Z)$ 時 $E[(Z-b)^2]$ 會最小，但是直覺來想，當我們要將許多點和某個點的差距取平方再相加時，相加的結果要最小，理應要發生在那個要和這些點算距離的點是他們的平均值時。

把這個想法套用到我們的 penalty $(\theta - w(y))^2$，就可以看出 $w(y)$ 應該要是 $E[\theta]$，也就是要把猜測的值訂為 $\theta$ 的 expected value。

再回過頭來看上方 $w(y)$ 的那個式子：

對 $\theta$ 從 $-\infty$ 積到 $\infty$ 的意思是，我們把所有可能的 $\theta$ 都考慮一遍，其中函數圖形為 $\theta k(\theta \ | \ y)$，也就是每個 $\theta$ 的值我們都去乘上 estimate $Y=y$ 的情況下，parameter 是 $\theta$ 的機率。

這樣的想法就和我們取 mean 相同，只不過是因為 $Y=y$ 的條件而變成了 conditional 而已。

同樣的道理，如果我們的 penalty (loss) function 訂為 error 取絕對值，也就是 $|\theta - w(y)|$，那我們的 $w(y)$ 就應該要訂為 distribution 的 median。
> distribution 的 median 也就可以由我們的 posterior pdf $k(\theta|y)$ 來算出。
>
> 至於為什麼取中位數會讓 penalty 最小，我也沒證，但像上面一樣，直覺也會想，如果要讓許多點和某個點之間的相對距離最小，那個點應該要在他們的正中間。

### 例一

舉個課本上篇幅很長的例子，假設我們有一個 random variable $Y$ 具 binomial distribution，我們的目標是根據如同上面所講到的步驟，最後求出 $w(y)$ 應該要令為什麼樣的 function 來估計 binomial distribution 的 parameter $\theta$。
> 關於 binomial distribution，很簡略的介紹可參考筆記「[A.3.2 Binomial Distribution](https://hackmd.io/@pipibear/SkWtX39NR)」。

因為過程太長了，很容易不知道自己到底在幹嘛，所以我們先把整個流程講一次：

---
1. 首先我們有一個 binomial distribution 的 pmf，pmf 中有一個未知的 parameter $\theta$ 是我們要估計的。
> 記得 binomial distribution 是在講 $n$ 次 Bernoulli trials 中有幾次 success，因此 sample space 為 $\{0,1,...,n\}$ 為 discrete，所以是 pmf 而非 pdf。
2. $\theta$ 具有 prior pdf $h(\theta)$。
> 在我們進行實驗前就預設的 $\theta$ 的分佈情形，在這個例子裡我們設定為一個 beta distribution。
>> $\rightarrow$ 關於 beta distribution 的介紹我還沒寫，課本其實也沒有多講，之後有機會再補。這裡會用到的我會簡單說明，但不會證明。
3. 透過上面的兩個式子，我們能夠去計算 joint pdf。
4. 我們前面說過，可以用 joint pdf 來定義 $Y$ 的 marginal pdf，於是我們求出只和 $Y$ 有關的 distribution。
5. 接著，根據 posterior pdf 的定義，將 joint pdf 除以 marginal pdf。
6. 有了 posterior pdf 以後，我們希望能夠讓 posterior expected loss 為最小。
> 因為此處假設我們要求 Bayes Estimator，而 Bayes Estimator 的定義就是做出讓 posterior expected loss 最小的選擇。
>
> 簡單一點來講，posterior expected loss 的意義是：
>
> 在我們預設了一種 $w(y)$ 的定義方式之後，$w(y)$ 和它所要估計的 $\theta$ 之間會有誤差，以前我們講過，我們有不同種方式去定義這個 loss （或是稱作 penalty），例如取兩者差距的絕對值，或是取兩者差距的平方，這裡題目假設用的是後者。
>
> 接著，我們就加總每種不同的 $\theta$ 值所產生的這個 loss，且因為我們已經知道 posterior pdf，所以我們知道各個 $\theta$ 值發生的機率，因此我們加總的是 loss 乘上它發生的機率。
>
> 那我們希望求出的 $w(y)$ 能讓這個 posterior expected loss 值最小。
7. 利用一些觀念求出 $w(y)$ 的定義以後，我們最後把它寫成另一種形式（利用一些簡單的代數把原本的東西拆成兩個），再去解讀它的意義。
8. 最後，我們簡單交代其實我們不用一開始就求 marginal pdf of $Y$，只需要知道我們要求的 posterior pdf 等同一個可以被視為常數的 $y$ 的函數，乘上一串東西；因此去求 posterior pdf，我們就等同去找那個常數 $c(y)$。
---

1. pmf 的定義

2. prior pdf $h(\theta)$
> 這個定義就是 beta distribution 的標準 pdf 定義。
>
> 透過兩個 parameters $\alpha$ 和 $\beta$，只要我們好好去選擇它們的值，這個 pdf 就夠有彈性。

3. 上面兩式相乘求 joint pdf

4. 求 $Y$ 的 marginal pdf

5. 由上面兩式相除求 posterior pdf

6. minimize posterior expected loss

7. 拆開 $w(y)$ 解釋意義
> 這裡我省略了兩件事的證明，分別是：
>
> 1. binomial distribution pdf 中的 parameter $\theta$，它的 maximum likelihood estimate 是 success 次數佔總次數的比例。
> 2. beta distribution 的 expected value 為 $\frac{\alpha}{\alpha + \beta}$
>
> 另外，$\alpha + \beta$ 在 beta distribution 中代表的意義是 ++sample size++，有機會寫 beta distribution 時這些再一起講。

不考慮這些問題的話，我們能夠看到 $w(y)$ 變成由 maximum likelihood estimate $\frac{y}{n}$ 乘上一個 weight $\frac{n}{\alpha + \beta + n}$，加上它的 mean $\frac{\alpha}{\alpha + \beta}$ 乘上一個 weight $\frac{\alpha + \beta}{\alpha + \beta + n}$。

$\rightarrow$ 因此我們發現，在選 $\alpha,\beta$ 的值時我們除了需要考慮讓 $\frac{\alpha}{\alpha + \beta}$ 會是一個適合的 prior mean，$\alpha + \beta$ 也要是 sample size。

舉個例子來看在這樣的情況下 $\alpha,\beta$ 的值會受到什麼樣的限制、有什麼樣的影響：

> 此處 posterior 為什麼也具 beta distribution、又為什麼具那樣的 $\alpha,\beta$ 值，之後有寫相關筆記再補上。

這個例子的圖：


8. 求 posterior 只需找出 constant $c(y)$ 使得 probability 總和為 $1$
> 我們發現其實不需要先去決定 marginal pdf of $Y$，因為它（$k_1(y)$）只和 $y$ 有關，和 $\theta$ 無關，所以我們可以將它視為常數。
>
> 這樣一來，$k(y,\theta)$ 去除它就等同於乘上一個常數，所以我們說 posterior $k(\theta|y)$ 正比於 $k(y,\theta)$，並且根據定義也就正比於 $g(y|\theta)h(\theta)$
>
> 至於為什麼我們的常數 $c(y)$ 要滿足 posterior pdf 積分會積成 $1$ 是因為必須要滿足 probability 的定義，pdf 底下的面積總和為 $1$。

### 例二

以下這個例子主要是在展示我們要怎麼利用剛剛得到的結果：

\begin{equation}
k(\theta|y) \propto g(y|\theta)h(\theta)
\end{equation}

假設我們一開始的條件如下：


首先，我們先將已知的條件代入剛剛結果的式子：


上面我們簡化到一半，如果我們把所有 constant 去除（反正我們寫成正比，所以並不影響），會得到有著一串指數的 $e$：

> 到這裡，我們發現如果能把 $e$ 上面那串指數配成 normal distribution pdf 中 $e$ 的指數 $\frac{(x-\mu)^2}{2\sigma^2}$ 這樣的形式，這樣我們就能找出 posterior mean, variance。

於是就接著一些很醜的計算：


最後去對照 normal distribution pdf 的樣子，我們順利求得 posterior mean, variance。

我們可以利用這個結果來知道要怎麼估計 $\theta$ 的值。

假設我們計算 error（估計的 $\theta$ 值和實際的 $\theta$ 之間的差距）的方式是用 squared-error，那前面有講過，$w(y) =$ posterior mean 時會有最小的 error。

這樣一來，我們就得到一個以 $y$ 為變數的公式，讓我們可以去算當 $y=$ 某個數字代入時，$\theta$ 應該要估為多少。

詳細如下：

> 其中如果我們把 posterior mean 拆成圖中的形式，會發現 posterior mean 就是 $y$ 乘上一個權重，再加上 $\theta_0$ 乘上另一個權重。因此，我們可以透過把 $n$ 取大一點或取小一點，來決定 prior 和 mle 之間哪個的影響要大一點。
>> 如果我們的 prior information 有告訴我們哪個東西的影響應該要比較大，那我們就能去決定 $n$ 的值該設大一點還是小一點。

:::info
Note：我們將 Bayes estimator 設為 posterior mean ，是因為例子裡用的 loss function 是 squared error；但如果我們的 loss function 是 error 取絕對值，也就是 $|w(y) - \theta|$，那就像我們前面說過的，Bayes estimator 就要改為 $k(\theta|y)$ 這個 posterior pdf 的 ++median++。
:::
> $\rightarrow$ <font color = "red">Bayes estimator 會隨著 loss function 的不同而改變！</font>

最後，關於這個例子，如果我們希望可以得到 interval estimate of $\theta$（去得到 $\theta$ 應該要落在哪個區間），那我們就會去找兩個 $y$ 的 functions $u(y)$ 和 $v(y)$，使得：

\begin{equation}
\int_{u(y)}^{v(y)} k(\theta|y) \,d\theta = 1 - \alpha
\end{equation}

其中 $\alpha$ 是一個很小的值。
> 把整個想法和這個式子的意思結合在一起講：
>
> 我們在猜測 $\theta$ 的值時，所有我們猜的值也可以視為一個 random variable $Y$，如果 $Y$ 的值是 $y$ 時（我們從所有猜測的值中選出一個 $y$），我們希望可以求出兩個 function，讓我們代入這個 $y$ 以後分別會給 output 一個 $\theta$ 的範圍，最小是 $u(y)$、最大是 $v(y)$。
>
> 因為 $k(\theta|y)$ 是 posterior probability，所以對每個可能的 $\theta$ 值，它都會 output 如果 $Y=y$ 的情況下，$\theta$ 是這個值的機率為多少。
>
> 因此，我們對它積分，從 $u(y)$ 積到 $v(y)$，意思就是我們在++加總當取到 $y$ 時，$\theta$ 的值落在這個範圍的機率，且這個機率的值會是 $1- \alpha$++。

假設我們令 $\alpha = 0.05$，也就是我們要求的範圍是 $95\%$ 的機率會落在這裡面。

這個例子裡，因為 posterior pdf 是 normal 的，所以透過查表我們可以知道要涵蓋 $95\%$ 的機率，我們的範圍要取在 mean 左右約 $1.96$ 個 standard deviation。

詳細如下圖：

> 因此我們就得到了兩個 function $u(y)$ 和 $v(y)$ 來告訴我們當 $Y=y$ 是什麼值時，我們的區間要取在哪才會讓 $\theta$ 有 $95\%$ 的機率會落在裡面。

---

這個小節的最後，我們要來說明，其實除了用一個 statistic $Y$，我們也可以選擇用我們的 sample observations $X_1,X_2,...,X_n$。

那我們原本用的 $Y$ 的 pdf $g(y|\theta)$ 就可以換成 likelihood function：

\begin{equation}
L(\theta) = f(x_1|\theta)f(x_2|\theta)...f(x_n|\theta)
\end{equation}

也就是 given $\theta$ 的情況下，$X_1,X_2,...,X_n$ 的 joint pdf。
> 因為這些 sample observations 彼此獨立，所以它們的 joint pdf 為個別的 pdf 相乘。
>> 一點點說明可參考筆記「[補充：Conditional Distributions](https://hackmd.io/@pipibear/HJB-y57wR)」最後的 "Independence" 部分。

而原本的 posterior pdf 式子 $k(\theta|y) \propto g(y|\theta)h(\theta)$ 就會變成：

\begin{equation}
k(\theta|x_1,x_2,...,x_n) \propto h(\theta)f(x_1|\theta)f(x_2|\theta)...f(x_n|\theta) = h(\theta)L(\theta)
\end{equation}

這樣一來，在我們有 data $x_1,x_2,...,x_n$ 的情況下，$k(\theta|x_1,x_2,...,x_n)$ 包含了所有關於 $\theta$ 的資訊。

因此，根據 loss function，我們就能像上面的例子那樣去選擇我們的 Bayes estimate of $\theta$（例如例子中的選 mean 或 median。）

有趣的是，我們會觀察到：

:::warning
如果我們的 loss function 設為：

「在真正的 $\theta$ 附近的 small neighborhood 為 $0$，否則，在這之外是一個很大的常數。」

那麼我們的 Bayes estimate $w(x_1,x_2,...,x_n)$ 就應該設為 conditional pdf $k(\theta|x_1,x_2,...,x_n)$ 的 mode。
:::
> 關於 mode 稍微清楚一點的定義寫在下面解釋的圖中。
> 在這裡我們看個圖簡單講一下：
> 
> 在這個圖中， $h(x)$ 的最大值為 $0.3$，發生在當 $x$ 介於 $(1.5,2.5)$ 這個區間時。
>
> 我們稱 interval $(1.5,2.5)$ 為 <font color = "snake">modal class</font>（有最大 class height 的 interval） $x=2$ 為 <font color = "snake">mode</font>（modal class 對應的 class mark）

所以再把上面的話講得更清楚一點，意思是在上述條件下，為了讓 expected loss 最小化，我們應該要把 $w(x_1,x_2,...,x_n)$ 設為讓 $k(\theta|x_1,x_2,...,x_n)$ 最大化時的 $\theta$。
> 更清楚一點的舉個例子：
>
> $k(\theta|x_1,x_2,...,x_n)$ 是當我們的 outcome $X_1=x_1,...,X_n = x_n$ 時 $\theta$ 在各個值的機率分佈，假如我們在 $\theta=1$ 時有最大的 probability：
>
> $k(\theta=1|x_1,x_2,...,x_n)=0.6$
>> 意思就是在 $X_1=x_1,...,X_n = x_n$ 時，$\theta$ 最有可能是 $1$，是 $1$ 的機率高達 $0.6$。
>
> 上面的話的意思就是，如果是這樣，那我們就應該把我們估計 $\theta$ 的值 $w(x_1,x_2,...,x_n)$ 設為 $1$。

理由如下圖解釋：

先把 loss function 用數學式表示：



# 參考資料
- Hogg,Tanis,Zimmerman, Probability and Statistical Inference, 9th ed(2015),p.233, 288-293
> Section 6.8 Bayesian Estimation
- wiki:
- [Bayes estimator](https://en.wikipedia.org/wiki/Bayes_estimator)
- [Beta distribution](https://en.wikipedia.org/wiki/Beta_distribution)
- [Binomial distribution](https://en.wikipedia.org/wiki/Binomial_distribution)