最大似然估计(Maximal Likelihood Estimation)

(一)基本思想

给定一个概率分布 D,假设其概率密度函数是 f_{D},它与一个未知参数 \theta 相关。我们可以从这个分布中抽取 n 样本 x_{1},x_{2},...,x_{n},我们就可以得到这个概率是

P(x_{1},...,x_{n}) = f_{D}(x_{1},...,x_{n}|\theta).

但是,在这里我们并不知道参数 \theta 的值。如何估计参数 \theta 的取值就成为了关键之处。一个简单的想法就是从这个分布中随机抽取样本 x_{1},...,x_{n},然后利用这些数据来估算 \theta 的值。

最大似然估计 (maximal likelihood estimator) 算法会计算参数 \theta 的最可能的值,也就是说参数的选择会使得这个采样的概率最大化

用数学的语言来说,首先我们需要定义似然函数

L(\theta) = f_{D}(x_{1},...,x_{n}|\theta),

并且在 \theta 的所有取值上,使得这个函数的取值最大化。换言之,也就是函数 L(\theta) 的一阶导数等于零。这个使得 L(\theta) 最大化的参数 \hat{\theta} 称为 \theta最大似然估计

Remark. 最大似然函数不一定是唯一的,甚至不一定是存在的。

(二)基本算法

求解最大似然函数估计值的一般步骤

(1)定义似然函数;

(2)对似然函数求导数,或者说对似然函数的对数求导数,目的都是为了更加方便地计算一阶导数;

(3)令一阶导数等于零,得到关于参数 \theta 的似然方程;

(4)求解似然方程,得到的参数就是最大似然估计。在求解的过程中,如果不能够直接求解方程的话,可以采取牛顿法来近似求解。

(三)例子

(i)Bernoulli 分布(Bernoulli Distribution)

假设我们有 n 个随机样本 x_{1},...,x_{n}. 如果第 i 个学生没有自行车,那么 x_{i}=0; 否则 x_{i}=1. 并且假设 x_{i} 是满足未知参数 p 的 Bernoulli 分布的。我们此时的目标是计算最大似然估计 p,也就是全体学生中拥有自行车的比例。

如果 \{x_{i}:1\leq i\leq n\} 是相互独立的 Bernoulli 随机变量,那么对每一个 x_{i} 而言,它的概率函数则是:

f(x_{i};p)=p^{x_{i}}(1-p)^{1-x_{i}}, \text{ for } x_{i}=0 \text{ or } 1 \text{ and } 0<p<1.

因此,似然函数 L(p) 可以定义为:

L(p)=\prod_{i=1}^{n}f(x_{i};p)=p^{\sum_{i=1}^{n}x_{i}}(1-p)^{n-\sum_{i=1}^{n}x_{i}}.

为了计算参数 p 的值,可以对 ln(L(p)) 求导:

\ln(L(p))=(\sum_{i=1}^{n}x_{i})\ln(p) + (n-\sum_{i=1}^{n}x_{i})\ln(1-p)

\frac{\partial\ln(L(p))}{\partial p} = \frac{\sum_{i=1}^{n}x_{i}}{p}-\frac{n-\sum_{i=1}^{n}x_{i}}{1-p}

\frac{\partial\ln(L(p))}{\partial p} = 0,可以得到

p=\sum_{i=1}^{n}x_{i}/n.

也就是说最大似然估计是

\hat{p}=\sum_{i=1}^{n}x_{i}/n.

(ii)Gaussian Distribution

假设 x_{1},...,x_{n} 满足正态分布,并且该正态分布的参数 \mu\sigma^{2} 都是未知的。目标是寻找均值 \mu 和方差 \sigma^{2} 的最大似然估计。

如果 \{x_{1},...,x_{n}\} 是满足正态分布的,那么对于每一个变量 x_{i} 的概率密度函数就是:

f(x_{i};\mu,\sigma^{2}) = \frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x_{i}-\theta)^{2}}{2\sigma^{2}}).

似然函数就是:

L(\mu,\sigma) = \prod_{i=1}^{n} f(x_{i};\mu,\sigma^{2}) = (2\pi)^{-n/2}\sigma^{-n}\exp(-\frac{\sum_{i=1}^{n}(x_{i}-\mu)^{2}}{2\sigma^{2}})

\ln(L(\mu,\sigma)) = -\frac{n}{2}\ln(2\pi) - n\ln(\sigma) - \frac{\sum_{i=1}^{n}(x_{i}-\mu)^{2}}{2\sigma^{2}}

\frac{\partial \ln(L(\mu,\sigma))}{\partial \mu} = \frac{\sum_{i=1}^{n}(x_{i}-\mu)}{\sigma^{2}}

\frac{\partial \ln(L(\mu,\sigma))}{\partial \sigma} = -\frac{n}{\sigma} + \frac{\sum_{i=1}^{n}(x_{i}-\mu)^{2}}{\sigma^{3}}

\frac{\partial \ln(L(\mu,\sigma))}{\partial \mu} =0 和  \frac{\partial \ln(L(\mu,\sigma))}{\partial \sigma}=0,可以求解方程组得到:

\hat{\mu}= \sum_{i=1}^{n}x_{i}/n,

\hat{\sigma}^{2} = \sum_{i=1}^{n}(x_{i}-\mu)^{2}/n.

(iii) Weibull 分布(Weibull Distribution)

首先,我们回顾一下 Weibull 分布的定义。Weibull 分布(Weibull Distribution)是连续型的概率分布,其概率密度函数是:

f(x;\lambda,k) = \frac{k}{\lambda}(\frac{x}{\lambda})^{k-1}\exp^{-(x/\lambda)^{k}} \text{ for } x\geq 0, f(x;\lambda,k)=0 \text{ for } x<0.

其中,x 是随机变量,\lambda>0 是 scale parameter,k>0 是 shape parameter。特别地,当 k=1 时,Weibull 分布就是指数分布;当 k=2 时,Weibull 分布就是 Rayleigh 分布。

Weibull 分布的累积分布函数

F(x;k,\lambda) = 1- \exp^{-(x/\lambda)^{k}} \text{ for } x\geq 0,

F(x;k,\lambda) = 0 \text{ for } x<0.

Weibull 分布的分位函数(quantile function, inverse cumulative distribution)是

Q(p;k,\lambda) = \lambda(-\ln(1-p))^{1/k} \text{ for } 0\leq p <1.

其次,我们来计算最大似然估计。

假设 \{x_{1},...,x_{n}\} 满足 Weibull 分布,其未知参数是 k,\lambda. 那么对于每一个 x_{i} 而言,概率密度函数是:

p(x_{i};k,\lambda) = \frac{k}{\lambda}(\frac{x_{i}}{\lambda})^{k-1}\exp(-(\frac{x_{i}}{\lambda})^{k}).

定义似然函数为:

L(k,\lambda) = \prod_{i=1}^{n}p(x_{i};k,\lambda)

取对数之后得到:

\ln(L(k,\lambda)) = n\ln(k) - nk\ln(\lambda) + (k-1)\sum_{i=1}^{n}\ln(x_{i}) - \sum_{i=1}^{n}x_{i}^{k}/\lambda^{k}.

计算一阶偏导数得到:

\frac{\partial \ln(L(k,\lambda))}{\partial \lambda} = - \frac{nk}{\lambda} + \frac{k\sum_{i=1}^{n}x_{i}^{k}}{\lambda^{k+1}},

\frac{\partial \ln(L(k,\lambda))}{\partial k} = \frac{n}{k} - n\ln(\lambda) + \sum_{i=1}^{n}\ln(x_{i}) -\sum_{i=1}^{n}(\frac{x_{i}}{\lambda})^{k}\ln(\frac{x_{i}}{\lambda}).

可以计算得出:

\lambda^{k}=\frac{\sum_{i=1}^{n}x_{i}^{k}}{n},

\frac{1}{k} = \frac{\sum_{i=1}^{n}x_{i}^{k}\ln(x_{i})}{\sum_{i=1}^{n}x_{i}^{k}} -\frac{\sum_{i=1}^{n}\ln(x_{i})}{n}.

其中第一个式子可以计算出 \lambda 的最大似然估计。第二个式子是关于 k 的隐函数,不能够直接求解,需要使用 Newton’s method 来计算。

f(k) = \frac{\sum_{i=1}^{n}x_{i}^{k}\ln(x_{i})}{\sum_{i=1}^{n}x_{i}^{k}} - \frac{\sum_{i=1}^{n}\ln(x_{i})}{n} - \frac{1}{k}.

求导得到:

f^{'}(k)= \frac{\sum_{i=1}^{n}x_{i}^{k}(\ln(x_{i}))^{2}}{\sum_{i=1}^{n}x_{i}^{k}}-(\frac{\sum_{i=1}^{n}x_{i}^{k}\ln(x_{i})}{\sum_{i=1}^{n}x_{i}^{k}})^{2} + \frac{1}{k^{2}}.

根据 Cauchy’s Inequality, 可以得到:

(\sum_{i=1}^{n}x_{i}^{k}(\ln(x_{i}))^{2})\cdot(\sum_{i=1}^{n}x_{i}^{k})\geq (\sum_{i=1}^{n}x_{i}^{k}\ln(x_{i}))^{2}.

所以,f^{'}(k)>0 \text{ for all } k>0,换言之,f(k) 是关于 k 的递增函数,并且

\lim_{k\rightarrow 0^{+}}f(k) = -\infty,

\lim_{k\rightarrow +\infty}f(k) > 0 \text{ if } \forall x_{i}>1.

那么对于递增函数 f(k) 而言,就必定有一个零点。因此使用 Newton’s Iteration 的时候,初始点可以从靠近零的整数开始,比如 k_{0}=0.0001。如果从一个比较大的数开始的时候,可能使用 Newton 法的时候,会与负轴相交。但是如果从一个较小的数开始,就必定只与正数轴相交。其中 Newton 法的公式是:

k_{0}= 0.0001,

k_{n+1} = k_{n}- \frac{f(k_{n})}{f^{'}(k_{n})} \text{ for all } n\geq 0.

n 的次数足够大的时候,k_{n} 就可以被当作最大似然估计。

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s