Question

Estimation（Parameter Estimation）

“估计”就是选定一种分布（例如常态分布、二项式分布、Gamma 分布），找到适当的参数（例如平均数、变异数、形状参数），尽量符合手边的一堆样本。

“估计”也是迴归，函数改成机率密度函数，数据改成样本。

Plot3D[PDF[MultinormalDistribution[{0, 0}, {{1, 0}, {0, 1}}], {x, y}], {x, -2, 2}, {y, -2, 2}, Boxed -> False, Axes -> False, Mesh -> Automatic, MeshFunctions -> {#3&}]

Estimator

以函数进行迴归，所谓最符合就是误差最小，例如 Least Squares。以分布进行估计，所谓最符合就是机率最大，例如 Maximum Likelihood、Maximum a Posterior。

一、已知样本，欲求分布的参数。

二、首先採用直接法，以样本直接推导参数。
　　然而十分困难，无法得到数学公式。

三、于是採用试误法，穷举各种参数，一一验证。
　甲、针对一种参数：
　　口、各个样本一一代入分布，分别求得出现机率。求得总乘积。
　　　　（此为 ML。若为 MAP，则额外乘上该分布参数的出现机率。）
　　口、记录总乘积最大者。
　乙、取 log，连乘化作连加。
　　　避免机率越乘越小、低于浮点数精确度而变成零。
　　　（取 log 最大值位置不变。）
　丙、穷举法太慢，需改用其他最佳化演算法。

范例：店面一日人潮统计

我准备好笔记本、手錶、提神饮料，从半夜 12 点开始，坐在麵店门口 24 小时，痴痴地等著客人上门，登记每位客人的到访时间，做为样本。我认为到访时间呈 Gamma 分布，我想估计平均数和形状参数是多少。

【待补图片】

平均数可能是 0 点、1 点、2 点、……、24 点，看起来每一种都有可能。这个时候就应该使用 Maximum Likelihood。

【待补图片】

不过根据我对真实世界的了解，我知道大家通常晚上六点下班，然后大家相约吃饭小酌一下。所以平均数是 19 点、20 点、21 点、22 点的机率非常高！我也知道三更半夜，不太有人吃麵，所以平均数是 2 点、3 点、4 点的机率非常低！我认为平均数呈常态分布！这个时候就应该使用 Maximum a Posterior。

【待补图片】

延伸阅读：Estimation

以统计学惯用的代数符号，重新说明“估计”。

一、已知一堆样本 X = {x1, ..., xN}。
　　已知特定分布的机率密度函数 f(x, μ, σ², λ)。
　　不知特定分布的参数Θ = {μ, σ², λ, ...}，
　　像是机率密度函数的平均数μ、变异数σ²、形状参数λ、...

二、统计学家习惯把已知与未知写成条件机率。
　　p( μ,σ,λ,... | x1,...,xN,f ) 或者 p( Θ | X,f )

三、所谓最符合，就是机率越大越好。
　　max p( μ,σ²,λ,...  | x1,...,xN,f ) 或者 max p( Θ | X,f )

四、找到此时平均数θ、变异数σ、形状参数λ是多少。
　　argmax p( μ,σ²,λ,... | x1,...,xN,f ) 或者 argmax p( Θ | X,f )
　　μ,σ²,λ,...                                  Θ

五、虽然我们知道 p 函数一定存在，但是我们不知道 p 函数长什麽样，无从计算。

Maximum Likelihood：找到其中一种分布参数，在此参数下，各个样本的机率，乘积最大。

一、ML 是找到其中一种分布参数，在此参数下，这堆样本的机率最大。
　　argmax f(X|Θ)
　　  Θ

二、推定样本之间互相独立、不互相影响，就可以套用乘法定理。
　　argmax f(X|Θ)
　　  Θ
　= argmax [ f(x1|Θ) * ... * f(xN|Θ) ]
　　  Θ

三、取 log 将连乘化作连加。取 log 后最大值位置仍相同。
　　argmax log f(X|Θ)
　　  Θ
　= argmax log [ f(x1|Θ) * ... * f(xN|Θ) ]
　　  Θ
　= argmax     [ log f(x1|Θ) + ... + log f(xN|Θ) ]
　　  Θ

四、求得函数 log f(X, Θ) 的最大值。

Maximum a Posterior：找到其中一种分布参数暨样本，出现机率的乘积最大。

一、MAP 是计算各个样本暨各个分布参数的出现机率，令总乘积最大。
　　argmax p(X,Θ)
　　  Θ

二、套用贝式定理
　　argmax p(X,Θ) = argmax { p(X|Θ) * p(Θ) }
　　  Θ               Θ

三、推定样本之间互相独立、不互相影响，就可以套用乘法定理。
　　argmax p(X,Θ) = argmax { p(X|Θ) * p(Θ) }
　　  Θ               Θ
　= argmax { f(x1|Θ) * ... * f(xN|Θ) * p(Θ) }
　　  Θ

四、后面都跟 Maximum Likelihood 的步骤完全一样，多了一项 p(Θ) 而已。

　　argmax log f(X,Θ) = argmax log { p(X|Θ) * p(Θ) }
　　  Θ                   Θ
　= argmax log { f(x1|Θ) * ... * f(xN|Θ) * p(Θ) }
　　  Θ
　= argmax     { log f(x1|Θ) + ... + log f(xN|Θ) + log p(Θ) }
　　  Θ

MAP 是 ML 的通例。ML 假设各种分布参数的出现机率均等，呈 uniform distribution。MAP 更加仔细考虑分布参数的出现机率，不见得要均等。

Model Selection / Model Validation

https://en.wikipedia.org/wiki/Model_selection

http://en.wikipedia.org/wiki/Regression_model_validation

我要怎麽知道一开始选择的分布是对的？我要如何判断到访时间比较像 Gamma 分布，或者比较像 Poisson 分布呢？

这属于统计学的范畴，就此打住。我没有研究。似乎大家都是自由心证。

AIC = -2*ln(likelihood) + 2*k,
BIC = -2*ln(likelihood) + ln(N)*k
k = model degrees of freedom
N = number of observations

Bias / Variance

假设一开始选定的分布是对的！如果不对，此段落没啥好谈。

Bias 是指“穷举各种样本组成，分别估计参数，所有结果的平均数”与“真实参数”的差值。

Bias 是衡量估计参数对不对的指标。不好的 Estimator，可以证明估计参数铁定失准。

Variance 就是变异数，此处我们是算“穷举各种样本组成，分别估计参数，所有结果的变异数”。

Variance 此处用来衡量估计参数的浮动范围。我们希望对于奇葩的样本组成，估计结果仍然差不多，浮动范围越小越好。

Bias 和 Variace 是两件事情。即便正确，还是可以有浮动范围。

仔细推导 Bias 和 Variance 的关係式。平方误差的平均数，由 Bias 和 Variance 组成。完美的估计，令平方误差达到极小值、为定值，而此时 Bias 和 Variance 此消彼长，鱼与熊掌不可兼得。

mean square error = (bias)² + variance

儘管我们不可能知道真实参数是多少，不过却可以得到鱼与熊掌不可兼得的结论：无论採用哪种 Estimator，Bias 和 Variance 无法同时令人满意。

演算法（样本平均数、样本变异数）

经典的分布，估计平均数、变异数，採用 Maximum Likelihood，一次微分等于零、二次微分大于零，推导公式解，公式解多半是所有样本的“母体平均数”、“母体变异数”。母体二字常省略。

μ = (x₁ + ... + xN) / N
σ² = [ (x₁ - μ)² + ... + (xN - μ)² ] / N

母体平均数没有问题，其 Bias 等于零，证明省略。

母体变异数则有问题，其 Bias 不是零，Maximum Likelihood 不可靠。分母设定成 N-1，其 Bias 才是零。证明省略。

这称做“样本平均数”、“样本变异数”。

x = (x₁ + ... + xN) / N
s² = [ (x₁ - x)² + ... + (xN - x)² ] / (N-1)

另外补充一下。根据大数定律，当样本无限多、样本互相独立（样本随机取得），则母体平均数趋近分布平均数。大可不必透过 Maximum Likelihood。

但是我查不到母体变异数趋近分布变异数的任何资料。

演算法（Expectation-Maximization Algorithm）

经典的分布，诸如二项式分布、常态分布，估计时採用 ML 或 MAP，可以推导确切的函数式子，甚至微分式子。运气好，推导公式解；运气差，套用最佳化演算法。

专为 ML 和 MAP 设计的最佳化演算法，找到机率最大值。

http://www.seanborman.com/publications/EM_algorithm.pdf

http://www.cs.cmu.edu/~awm/10701/assignments/EM.pdf

一、凹（凸）函数定义可以写成内插。
　　内插之后函数值必然上升（下降）。
　　注：凹函数的外观是随时向上凸，定义不太直觉。
二、机率函数的期望值就是一种内插！
　　如果机率函数是凹（凸）函数，
　　想求极值，那就好办，不断求期望值即可！
三、改变 ML 函数、移动 log 位置，变成一个凹函数。
　　证明此凹函数小于等于原式，是 ML 函数的下界。
四、凹函数求期望值、往上爬，函数值严格上升。
　　ML 函数的函数值必然同时跟著上升。
五、根据现在位置，
　　不断求一个新的凹函数，不断求期望值、往上爬。
　　最后就会得到区域极值，类似 Hill Climbing 演算法。

Normal Regression

Normal Distribution

机率论课程有教，就不多介绍了。

            e-(k-μ)²/2σ²
P(X = k) = ―――――――――――
              √2σ²

若选定“ Normal Distibution ”进行估计，参数μ为样本平均数，参数σ^2 为样本变异数。

Linear Regression with Normal Error

http://www.robots.ox.ac.uk/~fwood/teaching/W4315_Fall2011/Lectures/lecture_3/lecture_3.pdf

线性迴归，迴归函数添上误差项，误差项推定为常态分布。此即讯号学的“ Additive White Gaussian Noise ”。

现在要尽量符合数据。由于迴归函数含有常态分布，只好从迴归改为估计，採用 Maximum Likelihood。

由于是经典分布，式子漂亮，得直接运用“一次微分等于零”，推导公式解。

可以发现：添加误差项、採用 Maximum Likelihood 实施估计，未添加误差项、採用 Least Squares 实施迴归，两者结果恰好一致！

也就是说：普通的线性迴归，已经内建“误差项呈常态分布”的效果！我们直接实施普通的线性迴归即可，大可不必设定误差项，自寻烦恼。

网路上有人把这件事情解读成：Maximum Likelihood 是 Least Squares 的通例。虽然不那麽正确，但是也不能说完全错误。

Normal Regression

换个角度解释方才的事情。

一、每笔数据，最后一个维度，由不同的 1D Normal Distribution 产生。
　　每笔数据，前面几个维度，用于各个 1D Normal Distribution 的平均数 μ。
二、推定各个 1D Normal Distribution 的平均数 μ 呈线性成长。
　　μ = ax + b
三、推定各个 1D Normal Distribution 的变异数 σ² 皆相同。

迴归函数是线性函数，代表各个常态分布的平均数。

为何我们迫令变异数皆相同呢？N+1 个关係式（N 笔数据、平均数线性成长）、N+1 个未知数（N 种平均数、一种变异数），得到唯一解。

Poisson Regression

Log-linear Model

http://en.wikipedia.org/wiki/Generalized_linear_model

线性函数进行迴归，容易推导公式解。于是统计学家以线性函数为基础，产生各种函数。

此处要使用的是 log-linear model，取 log 之后是线性函数。

y = eax + b , log(y) = ax + b

Poisson Distribution

机率论课程有教，就不多介绍了。

P(X = k) = λk e-k / k!

若选定“ Poisson distibution ”进行估计，参数λ为样本平均数。

Poisson Regression

一、每笔数据，最后一个维度，由不同的 1D Poisson Distribution 产生。
　　每笔数据，前面几个维度，用于各个 1D Poisson Distribution 的平均数 λ。
二、推定各个 1D Poisson Distribution 的平均数 λ 呈指数成长。
　　λ = eax + b , log(λ) = ax + b

数据的维度可以推广到三维以上。

2D: (x, f(x))           ==> λ = eax + b
3D: (x, y, f(x,y))      ==> λ = eax + by + c
4D: (x, y, z, f(x,y,z)) ==> λ = eax + by + cz + d

数据可以分组。一组数据，共用一个 Poisson Distribution。

http://biostat.tmu.edu.tw/page/tmudata/ch9.pptx

Poisson Regression 推定平均数呈指数成长，符合真实世界的常见情况。当然也可以更换成别种成长方式。

演算法

【待补文字】

Logistic Regression

Logistic Function

“ Logistic Function ”初始成长缓慢、中途成长迅速、饱和成长缓慢，是真实世界的常见情况。

“ Logit Function ”是其反函数。

此函数源于物流学，理应翻译为“物流函数”，却翻译为“逻辑函数”。莫名奇妙。

Bernoulli Distribution

掷一次硬币，正面机率 p，反面机率 1-p。

P(X = 1) = p
P(X = 0) = 1-p

若选定“ Bernoulli Distribution ”进行估计，参数 p 为样本平均数。注意到：数据只能是 0 或 1。

Bernoulli Regression（Logistic Regression）

http://www.mc.vanderbilt.edu/gcrc/workshop_files/2004-11-12.pdf

一、每笔数据，最后一个维度，由不同的 1D Bernoulli Distribution 产生。
　　每笔数据，前面几个维度，用于各个 1D Bernoulli Distribution 的出现机率 p。
二、推定各个 1D Bernoulli Distribution 的出现机率 p 呈物流函数成长：
　　p = 1/(1+e-(ax + b))
　　另一种解释，推定其胜算（odds）呈指数成长：
　　p/(1-p) = eax + b
　　另一种解释，推定其 Logit Function 呈线性成长：
　　log(p/(1-p)) = ax + b , log(p) - log(1-p) = ax + b

演算法

採用 Maximum Likelihood，经过偏微分，可以找到某两个式子成立时，有最大值。

http://www.real-statistics.com/logistic-regression/basic-concepts-logistic-regression/
http://www.real-statistics.com/logistic-regression/finding-logistic-regression-coefficients-using-newtons-method/logistic-regression-using-newtons-method-detailed/

此二式难以推导公式解，故採用最佳化演算法 Newton's Method。

UVa 10727

Gaussian Process Regression

Gaussian Process Regression（Kriging）

多项式内插，引入机率的概念，让内插函数拥有浮动范围。

http://en.wikipedia.org/wiki/Kriging