精华内容
下载资源
问答
  • 最大似然估计的目的就是:利用已知的样本结果,反推最有可能(最大概率)导致这样结果的参数值。 原理:极大似然估计建立在极大似然原理的基础上的一个统计方法,概率论在统计学中的应用。 ... 记已知的样本集为...

    最大似然估计的目的就是:利用已知的样本结果,反推最有可能(最大概率)导致这样结果的参数值。 原理:极大似然估计是建立在极大似然原理的基础上的一个统计方法,是概率论在统计学中的应用。 ... 记已知的样本集为: 似然函数(linkehood function):联合概率密度函数 称为相对于 的θ的似然函数

     

    连续概率密度函数

    如果抛硬币10次,出现8次正面;最大似然估计就是反推正面的概率是多少,能使这种结果出现的可能性最大;在参数theta的所有取值(0到1无数个取值可能性)中,这里的参数theta就是正面向上的概率值;最大似然估计结果:theta=p(正面)=0.8;

    展开全文
  • 3.极大似然估计的核心让产生所采样的样本出现的概率最大。即利用已知的样本结果信息,反推具有最大可能导致这些样本结果出现的模型的参数值。 既然事情已经发生了,为什么不让这个出现的结果的可能性最大呢?这也...

    1.极大似然估计中采样产生的样本需要满足一个重要假设,所有采样的样本都是独立同分布的。
    2.极大似然估计是在模型已定,参数未知的情况下,估计模型中的具体参数。
    3.极大似然估计的核心是让产生所采样的样本出现的概率最大。即利用已知的样本结果信息,反推具有最大可能导致这些样本结果出现的模型的参数值。
    既然事情已经发生了,为什么不让这个出现的结果的可能性最大呢?这也就是最大似然估计的核心。
    求最大似然函数估计值的一般步骤:
    (1)写出似然函数;似然函数值的大小意味着这组样本值出现的可能性的大小,是个概率值。
    (2)对似然函数取ln对数,并整理化简;对数函数是单调增函数,所以对数函数取最大值时,原函数也取得最大值。(对数函数y=logax,当a>1时单调递增,当0<a<1时单调递减。)
    (3)求导数,令导数为0,得到似然方程;
    (4)解似然方程,得到的参数即为所求;


    下面使用极大似然法对逻辑斯蒂回归中的参数进行估计:
    1.假设当前的样本为{(x1,y1=1),(x2,y2=0),(x3,y3=1),(x4,y4=0),(x5,y5=0)},样本是满足独立同分布的。
    2.模型为P(Y=1|x)=11+exp(wx),表示一个样本x在该模型下预测值为1的概率,满足“模型已定,参数未知”的原则。
    3.由于样本是满足独立同分布的,那么出现以上样本分布的总概率为下式,需要让产生这一组样本的概率最大。

    P=P(Y=1|x=x1)P(Y=0|x=x2)P(Y=1|x=x3)P(Y=0|x=x4)P(Y=0|x=x5)

    P(Y=1|x)=π(x),P(Y=0|x)=1π(x),则上式可以化为:
    P=i=15[π(xi)]yi[1π(xi)](1yi)

    以上即为逻辑斯蒂回归的似然函数,似然函数的值的大小意味着该样本值出现的可能性的大小。需要最大化这个似然函数,使得出现这组样本的概率最大。可以将上式化为对数似然函数,之后将实际模型带入,求得似然函数在最大值时参数的实际值。具体参考:《统计学习方法》第六章逻辑斯蒂回归与最大熵模型学习笔记

    参考:
    一文搞懂极大似然估计
    最大似然估计的学习,这个对数学原理讲的很清楚

    展开全文
  • 线性模型是什么? 线性模型是监督学习各类学习算法最基本一类学习算法,它基本特性在于“线性”性质。线性是我们能处理一类最简单情况,同时其用处也非常广泛。其假设我们样本空间χ\chiχ到标记空间...

    简单线性模型及局部加权线性模型的基本原理和python实现(参数估计的两个基本角度:几何直观和概率直观。函数最值问题的两大基本算法:梯度方法与迭代方法)


    线性模型是什么?

    线性模型是监督学习的各类学习算法最基本的一类学习算法,所谓线性值得是模型的输出与数据的各个属性值之间的关系是线性的。线性是我们最擅长处理的(无论是工程实现还是数学分析),同时利用线性进行拓展,我们还可以处理很多非线性的情况。
    简单线性模型假设我们的样本空间χ\chi与参数θ\theta之间的关系是线性的,如下所示:
    y^=θ0x0+θ1x1+...+θnxn \hat{y} = \theta_0x_0+\theta_1x_1+...+\theta_nx_n记之为y^=hθ(x)=θTx\hat{y} = h_{\theta}(x) = \theta^Tx,简单线性模型是后面我们将学习的广义线性模型的基础。

    求解参数θ\theta的两种角度:几何直观和概率直观

    估计上述参数θ\theta的方式,主要有两种,一种是从几何直观,一种是从概率直观。它们分别是:

    • 最小二乘方法 (Ordinary Squart Least, OSL)
    • 极大似然估计 (Maximum Liklihood Estimate, MLE)
    1. 最小二乘方法从几何的角度出发,认为拟合值与原值之间的平方误差和应该最小化,因此我们定义了一个损失函数(cost function) 如下:
      J(θ)=12i=1m[hθ(x(i))y(i)]2J(\theta) = \frac{1}{2}\sum_{i=1}^{m}[h_{\theta}(x^{(i)})-y^{(i)}]^2 参数θ\theta需要使得损失函数最小化

    2. 极大似然估计是从概率的角度出发,认为样本数据既然发生了,其背后的概率应该最大化。因此我们首先对模型进行概率解释,即
      y=θ0x0+θ1x1+...+θnxn+ϵy = \theta_0x_0+\theta_1x_1+...+\theta_nx_n + \epsilon然后我们做一个比较合理的假设,即误差的零均值同方差正态假设,因此有
      ϵN(0,σ2)\epsilon \sim{N(0,\sigma^2)}所以标签值yy的分布为yN(θTx,σ2)y\sim{N(\theta^Tx,\sigma^2)}因此整个样本数据的概率似然函数为
      L(θ)=i=1m12πσexp((y(i)hθ(x(i)))2σ2)L(\theta) = \sum_{i=1}^{m}\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-h_{\theta}(x^{(i)}))^2}{\sigma^2})那么参数θ\theta需要使得似然函数最大化。

    以上两种参数估计的方法最终都归结为某个函数形式的最值问题,下面我们讨论如何求解非线性函数的最值。

    函数最值问题的两大数值算法

    • 梯度方法
      • 批量梯度下降(Batch Gradient Descent, BGD)
      • 随机梯度下降(Stochastic Gradient Descent, SGD)
    • 迭代方法
      • 一般迭代法
      • 牛顿迭代法

    梯度方法

    梯度是微积分中概念。在映射f:RnRf:R^n \rightarrow R中,它是一个nn维向量,其方向是该点处切平面增长最快的方向,在微分中,因为我们用切平面的增量近似函数的增量,所以梯度的方向也是该点函数增长最快的方向。从而利用自变量空间中关于ff的梯度场,就可以保证达到某个极值,在凸性的保证,还可以保证其达到最值。

    迭代方法

    我们详细介绍迭代方法。首先我们知道最值点的必要条件是稳定点,即导数为0。因此我们实际上需要一个寻找函数零点的方法。一般来讲,迭代方法可以分为两类:

    • 直接法(连续函数的中值定理)
    • 间接法

    直接法有我们熟悉的二分法等,我们略去不谈。这里要重点介绍间接法。我们首先谈迭代法的一般理论,然后再谈广泛使用的牛顿迭代法及其多元形式。
    在一维情况下,由f(x)=0f(x)=0可以尝试得到:
    x=φ(x)x=\varphi(x)φ(x)&lt;1\mid \varphi^{&#x27;}(x)\mid &lt; 1在区间[a,b][a,b]成立时,则可以证明在由迭代式xn=φ(xn1)x_{n}=\varphi(x_{n-1})产生的迭代序列{xn}\{x_{n}\}收敛于xx^*,其中x0[a,b]x_0\in[a,b]。从而有
    limnxn=limnφ(xn1)=x\lim_{n\to\infty}x_{n}=\lim_{n\to\infty}\varphi(x_{n-1})=x^*
    现在我们讨论牛顿迭代法。其基本思想在于做切线不断逼近零点。假设我们面对的函数为f(x)f(x),则在初始值x0x_0处的切线为yy0=f(x0)(xx0)y-y_0 = f^{&#x27;}(x_0)(x-x_0)y=0y=0得到x1x_1,从而不断迭代逼近真实值,其迭代式为
    xn=φ(xn1)=xn1f(xn1)f(xn1)x_n = \varphi(x_{n-1}) = x_{n-1}-\frac{f(x_{n-1})}{f^{&#x27;}(x_{n-1})}当求最值问题时,则迭代式为
    xn=xn1f(xn1)f(xn1)x_n = x_{n-1}-\frac{f_{&#x27;}(x_{n-1})}{f^{&#x27;&#x27;}(x_{n-1})}牛顿迭代法在一定的条件下可以证明收敛 (参见《计算方法》,凹凸性),其收敛条件一般采用事后判断法,即
    xnxn1&lt;ϵ\mid x_{n}-x_{n-1}\mid &lt; \epsilon当自变量为多元情况的时候,其迭代式的形式仍是一致的。根据多元微积分,此时的一阶偏导数为一个向量,即梯度,记为xf(x)\bigtriangledown_{x}f(x),二阶偏导数为一个矩阵,即海森矩阵(Hessian)HH,每一个元素为Hi,j=2f(x)xixjH_{i,j} = \frac{\partial^2{f(x)}}{\partial{x_{i}}\partial{x_{j}}}因此多元情况下的牛顿迭代法的迭代式为
    xn=φ(xn1)=xn1H1xf(x)x_n = \varphi(x_{n-1}) = x_{n-1}-H^{-1}\bigtriangledown_{x}f(x)注:一般情况下牛顿迭代法会比梯度下降迭代次数更少,收敛速度更快。然而每一次迭代牛顿迭代法需要需要计算的量更多,因为其中有一个海森矩阵。但只要维度n不是很大,总体上牛顿迭代的速度更快。

    Newton’s method typically enjoys faster convergence than (batch) gra-dient descent, and requires many fewer iterations to get very close to the minimum. One iteration of Newton’s can, however, be more expensive than one iteration of gradient descent, since it requires finding and inverting an n-by-n Hessian; but so long as n is not too large, it is usually much faster overall. When Newton’s method is applied to maximize the logistic regres-sion log likelihood function ℓ(θ), the resulting method is also called Fisher scoring.——Andrew Ng

    殊途同归:简单线性模型下几何直观和概率直观的统一汇合

    前面我们讨论了θ\theta的两个求解角度,得到了两个函数,即损失函数和似然函数。接下来我们将看到,最小化损失函数和最大化似然函数实际上是等价的。然后我们会讨论在简单线性模型下,两种角度下的函数最值求解,其实不同数值算法也可以,因为我们可以直接进行代数分析给出最值点的数学解析式。

    注意,从OSL和MLE两个角度出发,得到损失函数和似然函数,其最值问题等价,这个事情并不平凡。在后面介绍广义线性模型的时候,我们会看到基于各种指数簇分布的回归特例(简单线性模型只是基于正态分布下的回归特例),这个事情并不成立,即它们在理论上可能得到不同的参数θ\theta。BTW,后面我们采用的都是MLE的角度,因为从MLE出发,利用梯度方法,各个回归特例的梯度更新公式都有统一的优美形式,此为后话。

    几何直观角度下(OSL)的BGD、SGD梯度方法

    OSL定义的损失函数(cost function) 如下:
    J(θ)=12i=1m[hθ(x(i))y(i)]2J(\theta) = \frac{1}{2}\sum_{i=1}^{m}[h_{\theta}(x^{(i)})-y^{(i)}]^2其关于参数的导数为J(θ)θ=i=1m(hθ(x(i))y(i))x(i)\frac{\partial{J(\theta)}}{\partial\theta}=\sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}在BSD中,每一次对θ\theta的更新需要遍历全部训练集的数据,其更新规则如下:θ:=θ+αi=1m(y(i)hθ(x(i)))x(i)\theta:=\theta+\alpha\sum_{i=1}^{m}(y^{(i)}-h_{\theta}(x^{(i)}))x^{(i)}而在SGD中,每一次仅从训练集中挑出一个训练样本对θ\theta进行更新,因此其更新规则为:
    θ:=θ+α(y(i)hθ(x(i)))x(i)\theta:=\theta+\alpha(y^{(i)}-h_{\theta}(x^{(i)}))x^{(i)}其中α\alpha表示的是学习率,或者形象地认为其为更新的步长。
    注:从更高的角度来看,我们要降低不仅仅是训练数据集中的整体误差,而是整个样本空间的泛化误差,因此真正的损失函数为J(θ)=12xχ[hθ(x(i))y(i)]2p(x)dxJ(\theta) = \frac{1}{2}\int_{x\in\chi}[h_{\theta}(x^{(i)})-y^{(i)}]^2p(x)dx而由于事实上我们无法知道总体的分布,因此按照蒙特卡洛模拟算法,应该抽取一定的样本进行近似,若取样本量M=mM = m,则为BGD,若取样本量M=1M = 1时,则为SGD,显然还可以取1&lt;M&lt;m1 &lt; M &lt; m,此时则为小批量梯度下降算法 (Mini-Batch Gradient Descent)。
    至此其实可以看到BGD与SGD本质上是一致的,即从总体中抽出样本来降低模拟误差。实际中SGD的应用会更多一些,特别是在训练数据量m较大时,其收敛速度更快,同时一定程度上能避免陷入局部最小值。

    概率直观下(MLE)的 BGD、SGD方法

    概率直观下定义的似然函数为
    L(θ)=i=1m12πσexp((yihθ(x))2σ2)L(\theta) = \sum_{i=1}^{m}\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y_{i}-h_{\theta}(x))^2}{\sigma^2})对数之并约简,即可得到
    (θ)=lnL(θ)=mln(2πσ)1σ2i=1m(yihθ(x))2\ell(\theta) = lnL(\theta)=-mln(\sqrt{2\pi}\sigma)-\frac{1}{\sigma^2}\sum_{i=1}^{m}(y^{i}-h_{\theta}(x))^2可以看到,求(θ)\ell(\theta)的最大值,等价于求i=1m[hθ(x(i))y(i)]2\sum_{i=1}^{m}[h_{\theta}(x^{(i)})-y^{(i)}]^2的最小值,这正是OLS的损失函数的形式。因此MLE对参数θ\theta的求解,理论上同OLS就是等价的。由此也可以领悟到OLS是一个非常自然的方法,它实际上是在一组假设下做极大似然估计(我们其实更偏爱概率直观,所以用MLE来解释OSL)。

    This is thus one set of assumptions under which least-squares re-gression can be justified as a very natural method that’s just doing maximum likelihood estimation.----Andrew Ng

    两种角度下最优参数的数学解析式

    通过代数分析,我们可以直接在理论上给出θ=argminθΘJ/L(θ)\theta = {\arg\min}_{\theta\in\Theta}J/L(\theta)时的数学表达式。为了方便进行代数分析,我们引入下列记号:

    • 样本矩阵XX
      [1x1(1)x2(1)xn(1)1x1(2)x2(2)xn(2)1x1(3)x2(3)xn(3)1x1(m)x2(m)xn(m)] \left[ \begin{matrix} 1 &amp; x_1^{(1)} &amp; x_2^{(1)} &amp; \cdots &amp; x_n^{(1)}\\ 1 &amp; x_1^{(2)} &amp; x_2^{(2)} &amp; \cdots &amp; x_n^{(2)}\\ 1 &amp; x_1^{(3)} &amp; x_2^{(3)} &amp; \cdots &amp; x_n^{(3)}\\ \vdots &amp; &amp; &amp; \ddots &amp; \vdots\\ 1 &amp; x_1^{(m)} &amp; x_2^{(m)} &amp; \cdots &amp; x_n^{(m)} \end{matrix} \right]

    • 对参数向量θ\theta的求导运算:

    1. 向量求导:ATθθ=A\frac{\partial{A^T\theta}}{\partial{\theta}}=A,其中AA是一个常数列向量。
    2. 向量求导:θTBθθ=2Bθ\frac{\partial{\theta^TB\theta}}{\partial{\theta}}=2B\theta,其中BB是一个对称矩阵。

    至此,则损失函数可以表示为矩阵代数的形式:
    J(θ)=12u^Tu^=12(yXθ)T(yXθ)=12(yTy2yTXθ+θTXTXθ)\begin{aligned} J(\theta)=&amp; \frac{1}{2}\hat{u}^T\hat{u} \\ =&amp;\frac{1}{2}(y-X\theta)^T(y-X\theta)\\ =&amp;\frac{1}{2}( y^Ty-2y^TX\theta+\theta^TX^TX\theta) \end{aligned}J(θ)=XTXθXTy=0J^{&#x27;}(\theta) = X^TX\theta-X^Ty = 0,即可解出
    θ=(XTX)1XTy\theta = (X^TX)^{-1}X^Ty

    简单线性模型的进一步拓展:局部加权线性回归(Local Weighted Liner Regression,LWLR)

    局部加权线性回归采用一种非参数的方法来预测未知的样本。在之前的方法中,一旦参数向量θ\theta学习出来后,则可以将历史的样本抛弃。而在LWLR中,历史收集到的样本需要始终保留下来,每一次对新样本进行估计时,都会动态得到结合样本数据和新样本的位置,给出一个参数θ\theta的一个新的估计结果。换言之,θ\theta是为新样本量身定做的。下面我们主要从几何直观的角度讨论LWLR。
    注:我们其实仍然可以从概率直观的角度来解释LWLR,有兴趣的读者可以自己思考(从异方差的角度)。

    下面讨论LWLR的具体做法。对于一个新进的样本x(j)x^{(j)}

    • 利用已有的数据(经验),最小化i=1mw(i)(hθ(x(i))y(i))2\sum_{i=1}^{m}w^{(i)}(h_{\theta}(x^{(i)})-y^{(i)})^2,得到θj^\hat{\theta_j}
    • 给出预测值y^(j)=θj^Tx(j)\hat{y}^{(j)} = \hat{\theta_j}^{T}x^{(j)}

    其中,w(i)w^{(i)}可以由核光滑函数来确定,即:
    w(i)=exp((x(i)x(j))22τ2)w^{(i)} = exp(-\frac{(x^{(i)}-x^{(j)})^2}{2\tau^2})其中τ\tau称为宽度(bandwidth),控制了权重随距离下降的快慢。从几何直观的角度看,越靠近x(j)x^{(j)}x(i)x^{(i)},其误差平方将会被给予更大的权重。
    总体上可以认为对y(j)y^{(j)}的估计时,不同的训练样本不在被一视同仁,越靠近x(j)x^{(j)}的样本点越被重视,从而得到的直线将会更好地符合x(j)x^{(j)}周围的样本点的情况。

    在局部加权最小二乘法中,利用梯度方法进行θ\theta的迭代求解,此时的更新公式是:

    • BGD:θ:=θ+αi=1m(y(i)hθ(x(i)))x(i)w(i)\theta:=\theta+\alpha\sum_{i=1}^{m}(y^{(i)}-h_{\theta}(x^{(i)}))x^{(i)}w^{(i)}
    • SGD:θ:=θ+α(y(i)hθ(x(i)))x(i)w(i)\theta:=\theta+\alpha(y^{(i)}-h_{\theta}(x^{(i)}))x^{(i)}w^{(i)}

    利用矩阵代数分析,同样可以得到θj^\hat{\theta_{j}}的数学解析式。首先令
    权重矩阵WW
    [w(1)000w(2)000w(m)] \left[ \begin{matrix} w^{(1)} &amp; 0 &amp; \cdots &amp; 0\\ 0 &amp; w^{(2)} &amp; \cdots &amp; 0\\ \vdots &amp; &amp; \ddots &amp; \vdots \\ 0 &amp; 0 &amp; \cdots &amp; w^{(m)} \end{matrix} \right] 则此时的损失函数可以表示为:
    J(θ)=12i=1mw(i)(hθ(x(i))y(i))2=12i=1m(w(i)hθ(x(i))w(i)y(i))2=12i=1m(hθ(w(i)x(i))w(i)y(i))2\begin{aligned} J(\theta) = &amp; \frac{1}{2}\sum_{i=1}^{m}w^{(i)}(h_{\theta}(x^{(i)})-y^{(i)})^2\\ = &amp; \frac{1}{2}\sum_{i=1}^{m}(\sqrt{w^{(i)}}h_{\theta}(x^{(i)})-\sqrt{w^{(i)}}y^{(i)})^2\\ = &amp; \frac{1}{2}\sum_{i=1}^{m}(h_{\theta}(\sqrt{w^{(i)}}x^{(i)})-\sqrt{w^{(i)}}y^{(i)})^2 \end{aligned}因此原来的样本数据等价变换后的样本数据,其属性值为w(i)x(i)\sqrt{w^{(i)}}x^{(i)},标记值为w(i)y(i)\sqrt{w^{(i)}}y^{(i)}。利用矩阵的记法,则样本矩阵记为W12XW^{\frac{1}{2}}X,标记矩阵记为W12yW^{\frac{1}{2}}\vec{y},故代入原来的解析式
    θ=(XTX)1XTy\theta = (X^TX)^{-1}X^T\vec{y}可得
    θj^=(XTWX)1XTWy\hat{\theta_{j}} = (X^TWX)^{-1}X^TW\vec{y}一般的,当样本属性空间是多维度的时候,即xx是多维向量时,根据高维欧式空间距离的推广式:xy=(xy)T(xy) \parallel x-y \parallel = (x-y)^T(x-y) 核光滑函数可以表示为:
    w(i)=exp((x(i)x(j))T(x(i)x(j))2τ2)w^{(i)} = exp(-\frac{(x^{(i)}-x^{(j)})^T(x^{(i)}-x^{(j)})}{2\tau^2})

    实战项目:简单线性模型下BGD、SGD、Mini-Batch GD、数学解析式的实现和LWLR的实现 (python)

    特别说明:在梯度下降算法中,需要特别留意学习率α\alpha(步长)的大小设置。以BGD为例,如果α\alpha设置过大,则会出现总误差不减反增的现象,最终导致溢出。通过查看中间输出可看到,其过程中误差函数的偏导数出现正负号轮流出现的规律,可以推断这是由于α\alpha过大导致的震荡现象,从而使得偏导数越来越大,总误差趋向无穷。
    SGD中步长更合适的方法应该是采取自适应方式,即不应当将其设为一个常值,原因在于在最后收敛处,较大的步长会影响其得到稳定的优解,而在一开始迭代时,较小的步长会影响到其迭代的速度,因此总体上步长的设置应该逐步减小。

    myLm.py

    import numpy as np
    from numpy import dot
    from numpy.linalg import inv
    import random
    
    def Lm(X,Y,method,alpha = 0.001,epsilon = 0.0000001,o = None,tau = None):
        '''
        X:样本矩阵 [1,x1,x2,x3,...,xn]
        Y:样本标签值向量 [y]
        method:训练参数时使用的方法,有OLS、BGD、SGD、MBGD、LWLM
        alpha: 学习速率(步长)
        epsilon:迭代终止阈值
        return:本函数通过使用相应的算法,给出简单线性回归模型的参数估计值
        '''
        if(method == 'OLS'):
            theta_0 = dot(dot(inv(dot(X.T,X)),X.T),Y)
            return theta_0
        elif(method == 'BGD'):
            #初始值设置
            theta_0 = np.ones(shape = X.shape[1])
            #计算当前损失函数的总误差
            error_0 = 0.5*dot((dot(X,theta_0)-Y).T,dot(X,theta_0)-Y)
            #启动循环
            error = epsilon+1
            #迭代次数n
            n = 0
            while(error > epsilon):       
                #计算全部样本误差函数的梯度
                G = dot(dot(X.T,X),theta_0)-dot(X.T,Y)
                #更新theta
                theta_0 = theta_0 - alpha*G
                #计算步长更新后的损失函数的总误差
                error_1 = 0.5*dot((dot(X,theta_0)-Y).T,dot(X,theta_0)-Y)
                #计算前后二者的差
                error = abs(error_1-error_0)
                #将新一轮的误差保存用于下一轮
                error_0 = error_1
                # 输出迭代过程
                # print('---------------第',n,'轮迭代--------------')
                # print('theta:',theta_0)
                # print('error:',error)
                # print('G:',G,'\n')
                n = n + 1
            print('BGD Iteration Times: ',n)
            return theta_0
            
        elif(method == 'SGD'):
            #初始值设置
            theta_0 = np.ones(shape = X.shape[1])
            #计算当前初始值下损失函数的总误差
            error_0 = 0.5*dot((dot(X,theta_0)-Y).T,dot(X,theta_0)-Y)
            #启动循环
            error = epsilon+1
            #样本量
            m = X.shape[0]
            #遍历列表
            k = list(range(m))
            #迭代次数
            n = 0
            while(error > epsilon):
                #选择一个随机数[] r = random.randint(0,m-1)
                #遍历乱序列表
                random.shuffle(k)
                #计算当前样本下的梯度
                for r in k:
                    G = (dot(X[r,:],theta_0)-Y[r])*X[r,:]
                    #更新theta的值
                    theta_0 = theta_0 - alpha*G
                #计算步长更新后theta_1下损失函数的总误差
                error_1 = 0.5*dot((dot(X,theta_0)-Y).T,dot(X,theta_0)-Y)
                #计算误差修正值大小
                error = abs(error_1-error_0)
                #将新一轮的误差保存用于下一轮
                error_0 = error_1
                # 输出迭代过程
                # if(n%100 == 0):
                    # print('---------------第',n,'轮迭代--------------')
                    # print('theta:',theta_0)
                    # print('error:',error)
                    # if(n > 8000):
                        # alpha = 0.00001
                    # print('r:',r)
                    # print('G:',G,'\n')
                n = n + 1
            print('SGD Iteration Times: ',n)
            return theta_0
        elif(method == 'MBGD'):
            #初始值设置
            theta_0 = np.ones(shape = X.shape[1])
            #计算当前损失函数的总误差
            error_0 = 0.5*dot((dot(X,theta_0)-Y).T,dot(X,theta_0)-Y)
            #启动循环
            error = epsilon+1
            #样本量
            m = X.shape[0]
            #遍历列表
            k = list(range(m))        
            #迭代次数n
            n = 0
            while(error > epsilon):       
                #随机选取一定批量的样本
                random.shuffle(k)
                minBatch = k[:m//10]
                minX = X[minBatch]
                minY = Y[minBatch]  
                #计算小批量样本误差函数的梯度
                G = dot(dot(minX.T,minX),theta_0)-dot(minX.T,minY)
                #更新theta
                theta_0 = theta_0 - alpha*G
                #计算步长更新后的损失函数的总误差
                error_1 = 0.5*dot((dot(X,theta_0)-Y).T,dot(X,theta_0)-Y)
                #计算前后二者的差
                error = abs(error_1-error_0)
                #将新一轮的误差保存用于下一轮
                error_0 = error_1
                # 输出迭代过程
                # print('---------------第',n,'轮迭代--------------')
                # print('theta:',theta_0)
                # print('error:',error)
                # print('G:',G,'\n')
                n = n + 1
            print('MBGD Iteration Times: ',n)
            return theta_0
        elif(method == 'LWLM'):
        	#局部加权线性回归
            XCopy = X[:,1:]-o
            d = []
            for i in range(XCopy.shape[0]):
                d.append(dot(XCopy[i],XCopy[i]))
            d = np.array(d)
            w = np.exp(-d**2/(2*tau**2))
            W = np.diag(w)
            theta_0 = dot(dot(dot(inv(dot(dot(X.T,W),X)),X.T),W),Y)
            return theta_0
        else:
            print('No Such Method In this Func.')
    

    myLmtest.py

    import myLm
    import numpy as np
    import pandas as pd
    
    if __name__ == '__main__':
        samData = pd.read_table('mlData/regreData/simpleRegre.txt',header = None)
        #转换为二维数组
        sample = np.array(samData)
        #将属性值与标签值分开
        sampleX = sample[:,:np.shape(sample)[1]-1]
        sampleY = sample[:,sample.shape[1]-1]
        x = np.array([0])
        #Testing
        theta = myLm.Lm(sampleX,sampleY,method = 'LWLM',o = x,tau = 1)
        print('LWLM:',theta)
        theta = myLm.Lm(sampleX,sampleY,method = 'SGD')
        print('SGD:',theta)
        theta = myLm.Lm(sampleX,sampleY,method = 'BGD')
        print('BGD:',theta)
        theta = myLm.Lm(sampleX,sampleY,method = 'MBGD')
        print('MBGD:',theta)
        theta = myLm.Lm(sampleX,sampleY,method = 'OLS')
        print('OLS:',theta)
    
    展开全文
  • 在统计计算中,最大期望(EM)算法在概率模型中寻找参数最大似然估计或者最大后验估计的算法,其中概率模型依赖于无法观测的隐性变量。最大期望算法经常用在机器学习和计算机视觉得数据聚类领域。EM算法...

    一、简介

    EM算法

    最大期望算法(Expectation-maximization algorithm,简称EM,又译期望最大化算法)在统计中被用于寻找依赖于不可观察的隐性变量的概率模型中,参数的最大似然估计。在统计计算中,最大期望(EM)算法是在概率模型中寻找参数最大似然估计或者最大后验估计的算法,其中概率模型依赖于无法观测的隐性变量。最大期望算法经常用在机器学习和计算机视觉得数据聚类领域。EM算法的标准计算框架由E步和M步交替组成,算法的收敛可以确保迭代至少逼近局部极大值。EM算法被广泛应用于处理数据的缺失值,以及很多机器学习算法,包括高斯混合模型和隐马尔科夫模型的参数估计。

    Jensen不等式(又名琴生不等式、詹森不等式)

    琴生不等式是以丹麦技术大学数学家约翰·延森(Johan Jensen)命名,它给出积分的凸函数值和凸函数的积分值间的关系。Jensen不等式是关于凸函数性质的不等式,它和凸函数的定义是息息相关的。

    47c3d7ea3433

    设f是定义域为实数的函数,如果对于所有的实数x,f(x)的二次导数大于等于0,那么f是凸函数

    如果f是凸函数,X是随机变量,那么满足Jensen不等式:

    math?formula=E%5Bf(x)%5D%20%5Cgeq%20f(E%5Bx%5D)

    Jensen不等式应用于凹函数时,不等号方向反向。

    二、算法原理

    1、实例解析

    47c3d7ea3433

    这里有两个硬币分别是A和B,我们进行抛硬币的实验,假设这两个硬币做过手脚,出现正反面的概率不同,初始假设为:

    A:0.6几率正面

    B:0.5几率正面

    那么我们分别抛10次硬币,会出现5正5反的概率:

    选择硬币A的概率:

    选择硬币B的概率:

    math?formula=1-P(A)%3D0.55

    我们跟着箭头一步一步看:

    第一步:给出了初始假设,选择两个硬币的正面概率分别是0.6和0.5

    第二步:E-step

    1、初始值选择A的概率是:0.45,此时有5次H,5次T,可以计算出

    math?formula=E(H)%3D0.45%20%5Ctimes%205%20%3D%202.2H%EF%BC%8CE(T)%3D0.45%20%5Ctimes%205%20%3D2.2T

    选择B的概率是:0.55,此时有5次H,5次T,可以计算出

    math?formula=E(H)%3D0.55%20%5Ctimes%205%20%5Capprox%202.8H%EF%BC%8CE(T)%3D0.55%20%5Ctimes%205%20%5Capprox%202.8T

    2、先来看一下投掷出9次正面1次反面的概率

    math?formula=P(A)%3DC_%7B10%7D%5E9(0.6%5E9)(0.4%5E1)%20%5Capprox%200.04%20%EF%BC%8CP(B)%3DC_%7B10%7D%5E9(0.5%5E9)(0.5%5E1)%20%5Capprox%200.0098

    选择硬币A的概率

    math?formula=%5Cfrac%7BP(A)%7D%7BP(A)%2BP(B)%7D%3D0.8%2C选择硬币B的概率

    math?formula=1-P(A)%3D0.2

    我们分别计算硬币A正反面的期望:

    math?formula=E(H)%3D%200.8%20%5Ctimes%209%20%3D7.2H%2CE(T)%3D0.8%20%5Ctimes%201%20%3D0.8T

    硬币B正反面的期望:

    math?formula=E(H)%3D%200.2%20%5Ctimes%209%20%3D1.8H%2CE(T)%3D0.2%20%5Ctimes%201%20%3D0.2T

    3、依次计算出所有的样本,并分别将A和B硬币出现H和T的期望求和,

    A硬币:

    math?formula=%5Csum_%7BH%7D%3D(2.2%2B7.2%2B5.9%2B1.4%2B4.5)H%3D21.3H

    math?formula=%5Csum_%7BT%7D%3D(2.2%2B0.2%2B1.5%2B2.1%2B1.9)T%3D8.6T

    B硬币:

    math?formula=%5Csum_%7BH%7D%3D(2.8%2B1.8%2B2.1%2B2.6%2B2.5)H%3D11.7H

    math?formula=%5Csum_%7BH%7D%3D(2.8%2B0.2%2B0.5%2B3.9%2B1.1)H%3D8.4H

    第三步:M-step

    通过E-step计算的结果分别求出两个硬币的分布,

    math?formula=%5Cwidehat%20%5Ctheta_%7BA%7D%5E%7B(1)%7D%20%3D%5Cfrac%7B21.3%7D%7B21.3%2B8.6%7D%20%5Capprox%200.71

    math?formula=%5Cwidehat%20%5Ctheta_%7BB%7D%5E%7B(1)%7D%20%3D%5Cfrac%7B11.7%7D%7B11.7%2B8.4%7D%20%5Capprox%200.58

    将上述两个结果带入第一步进行迭代,不断更新P(A)和P(B),直到结果收敛为止就得到了我们需要的结果。

    2、公式推导

    (1) 问题:样本集

    math?formula=%7B%5Clbrace%20x_%7B1%7D%2Cx_%7B2%7D%2C...%2Cx_%7Bm%7D%5Crbrace%7D,包含m个独立的样本,其中每个样本

    math?formula=i对应的类别

    math?formula=z_%7B(i)%7D是未知的,所以很难用最大似然求解

    math?formula=l(%5Ctheta)%3D%5Csum_%7Bi%3D1%7D%5E%7Bm%7D%5Clog%20p(x_%7Bi%7D%3B%5Ctheta)%3D%5Csum_%7Bi%3D1%7D%5E%7Bm%7D%20%5Clog%20%5Csum_%7BZ%7D%20p(x_%7Bi%7D%2Cz%3B%5Ctheta)

    上式中,要考虑每个样本在各个分布中的情况,本来正常求偏导就可以了,但是现在

    math?formula=%5Clog后面还有求和,这就难解了。

    (2) 右式分子分母同时乘

    math?formula=Q(z)(

    math?formula=Q(z)

    math?formula=Z的分布函数)

    math?formula=%5Clog%20%5Csum_%7BZ%7D%20p(x_%7Bi%7D%2Cz%3B%5Ctheta)%3D%5Clog%20%5Csum_%7BZ%7DQ(z)%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D为什么这么做?目的是凑Jensen不等式

    (3) 利用Jensen不等式求解

    由于

    math?formula=%5Csum_%7BZ%7DQ(z)%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D

    math?formula=%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D的期望

    假设

    math?formula=Y%3D%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D ,则

    math?formula=%5Clog%20%5Csum_%7Bz%7DQ%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ%7D%3D%5Clog%20%5Csum_%7BY%7DP(Y)Y%3D%5Clog%20E(Y)

    此时:

    math?formula=%5Clog%20E(Y)%20%5Cgeq%20E(%5Clog%20Y)%20%3D%5Csum_%7BY%7DP(Y)%5Clog%20Y%3D%5Csum_%7BZ%7DQ(z)%5Clog%20%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D

    (4) 结论

    math?formula=l(%5Ctheta)%3D%5Csum_%7Bi%3D1%7D%5E%7Bm%7D%5Clog%5Csum_%7BZ%7Dp(x_%7Bi%7D%2Cz%3B%5Ctheta)%20%5Cgeq%20%5Csum_%7Bi%3D1%7D%5E%7Bm%7D%20%5Csum_%7BZ%7DQ(z)%5Clog%20%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D

    下界比较好求,所以我们要优化这个下界来使得似然函数最大。

    47c3d7ea3433

    如何能使得等式成立呢?(取等号)

    根据Jensen中等式成立的条件是随机变量是常数,则:

    math?formula=Y%3D%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BQ(z)%7D%3DC

    由于

    math?formula=Q(z)是z的分布函数,则关于分布函数的和一定为1:

    math?formula=%5Csum_%7Bz%7DQ(Z)%3D%5Csum_%7BZ%7D%20%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BC%7D%3D1

    由上式可得,所有的分子和等于常数C(分母相同),即C就是

    math?formula=p(x_%7Bi%7D%2Cz%3B%5Ctheta)对z求和

    math?formula=Q(Z)%3D%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7BC%7D%3D%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7B%5Csum_%7Bz%7Dp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%3D%5Cfrac%7Bp(x_%7Bi%7D%2Cz%3B%5Ctheta)%7D%7Bp(x_%7Bi%7D)%7D%3Dp(z%7Cx_%7Bi%7D%3B%5Ctheta)

    math?formula=Q(z)代表第

    math?formula=i个数据是来自

    math?formula=z_%7Bi%7D的概率,在固定

    math?formula=%5Ctheta后,使下界拉升的

    math?formula=Q(z)的计算公式,解决了

    math?formula=Q(z)如何选择的问题。这一步就是E-step,建立

    math?formula=l(%5Ctheta)的下界。接下来的M-step,就是在给定

    math?formula=Q(z)后,调整

    math?formula=%5Ctheta,去极大化

    math?formula=l(%5Ctheta)的下界。

    3、总结

    EM算法流程

    1、初始化分布参数

    math?formula=%5Ctheta

    2、E-step:根据参数

    math?formula=%5Ctheta计算每个样本属于

    math?formula=z_%7Bi%7D的概率(也就是我们的Q)

    3、M-step:根据Q,求出含有

    math?formula=%5Ctheta的似然函数的下界并最大化它,得到新的参数

    math?formula=%5Ctheta

    4、不断地迭代更新下去,直到收敛。

    三、GMM(高斯混合模型)

    GMM(Gaussian mixture model),高斯混合模型,也可以简写成MOG.高斯模型就是用高斯概率密度函数(正态分布曲线)精确地量化事物,将一个事物分解为若干的基于高斯概率密度函数(正态分布曲线)形成的模型。GMM已经在数值逼近、语音识别、图像分类、图像去噪、图像重构、故障诊断、视频分析、邮件过滤、密度估计、目标识别与跟踪等领域取得了良好的效果。实际上,GMM的目的就是找到一个合适的高斯分布(也就是确定高斯分布的参数μ,Σ),使得这个高斯分布能产生这组样本的可能性尽可能大(即:拟合样本数据)。高斯混合模型也​被视为一种聚类方法,是机器学习中对“无标签数据”进行训练得到的分类结果。其分类结果由概率表示,概率大者,则认为属于这一类。

    1、算法适合的情形:

    数据内部有多个类别,每个类别的数据呈现不同的分布,这时我们就可以使用GMM算法,尤其是有隐变量,隐变量有不同的分布。

    2、公式推导

    一般化的描述为:假设混合高斯模型由K个高斯模型组成(即数据包含K个类),则GMM的概率密度函数如下:

    math?formula=p(x)%3D%5Csum_%7Bk%3D1%7D%5E%7BK%7Dp(k)p%7B(x%7Ck)%7D%3D%5Csum_%7Bk%3D1%7D%5E%7BK%7D%5Cpi_%7Bk%7DN(x%7C%5Cmu_%7Bk%7D%2C%5Csum%7Bk%7D)其中,

    math?formula=p(x%7Ck)%3DN(x%7C%5Cmu_%7Bk%7D%2C%5Csum%20k)是第k个高斯模型的概率密度函数,可以看成选定第K个模型后,该模型产生x的概率;

    math?formula=P(k)%3D%5Cpi%20k是第k个高斯模型的权重,称作选择第K个模型的先验概率,且满足,

    math?formula=%5Csum_%7Bk%3D1%7D%5E%7BK%7D%5Cpi_%7Bk%7D%3D1

    概率密度知道后,我们看一下最大化对数似然函数的意义:

    首先直观化地解释一下最大化对数似然函数要解决的是什么问题。

    假设我们采样得到一组样本yt ,而且知道变量Y服从高斯分布(本文只提及高斯分布,其他变量分布模型类似),数学形式表示为Y∼N(μ,Σ)。采样的样本如下图所示,我们的目的就是找到一个合适的高斯分布(也就是确定高斯分布的参数μ,Σ),使得这个高斯分布能产生这组样本的可能性尽可能大。那怎么找到这个合适的高斯分布呢(在如下图中的表示就是1~4哪个分布较为合适)?这时候似然函数就闪亮登场了。

    47c3d7ea3433

    似然函数数学化:设有样本集

    math?formula=Y%3Dy_%7B1%7D%2Cy_%7B2%7D%2C...%2Cy_%7BN%7D%2Cp(y_%7Bn%7D%7C%5Cmu%2C%CE%A3)是高斯分布的概率分布函数,表示变量

    math?formula=Y%3Dy_%7Bn%7D的概率。假设样本的抽样是独立的,那么我们同时抽到这N个样本的概率是抽到每个样本概率的乘积,也就是样本集Y的联合概率。此联合概率即为似然函数:

    math?formula=L(%5Cmu%2C%CE%A3)%3DL(y_%7B1%7D%2Cy_%7B2%7D%2C...%2Cy_%7BN%7D%3B%5Cmu%2C%CE%A3)%3D%5Cprod_%7Bn%3D1%7D%5E%7BN%7Dp(y_%7Bn%7D%3B%5Cmu%2C%CE%A3)

    对上式进行求导并令导数为0(即最大化似然函数,一般还会先转化为对数函数再最大化),所求出的参数就是最佳的高斯分布对应的参数。所以最大化似然函数的意义就是:通过使得样本集的联合概率最大来对参数进行估计,从而选择最佳的分布模型。

    我们尝试使用最大化对数似然函数来解决GMM模型,

    math?formula=lnL(%5Cmu%2C%CE%A3%2C%5Cpi)%3D%5Csum_%7Bn%3D1%7D%5E%7BN%7Dln%5Csum_%7Bk%3D1%7D%5E%7BK%7D%5Cpi_%7Bk%7DN(y_%7Bn%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D)然后求导,令导数为0,就可以得到模型参数(μ,Σ,π) ,然而仔细观察可以发现,对数似然函数里面,对数里面还有求和。实际上没有办法通过求导的方法来求这个对数似然函数的最大值。此时,就需要借助EM算法来解决此问题了

    接下来我们使用EM算法求解GMM模型,

    EM算法可以用于解决数据缺失的参数估计问题(隐变量的存在实际上就是数据缺失问题,缺失了各个样本来源于哪一类的记录)。

    引入隐变量

    math?formula=r_%7Bjk%7D,它的取值只能是1或者0。

    取值为1:第j个观测变量来自第k个高斯分量

    取值为0:第j个观测变量不是来自第k个高斯分量

    那么对于每一个观测数据

    math?formula=y_%7Bi%7D都会对应于一个向量变量,

    math?formula=%5CGamma_%7Bj%7D%3D%7Br_%7Bj1%7D%2C...%2Cr_%7Bjk%7D%7D,那么有:

    math?formula=%5Csum_%7Bk%3D1%7D%5E%7BK%7Dr_%7Bjk%7D%3D1

    math?formula=p(%5CGamma_%7Bj%7D)%3D%5Cprod_%7Bk%3D1%7D%5E%7BK%7D%5Calpha_%7Bk%7D%5E%7Br_%7Bjk%7D%7D

    其中,K为GMM高斯分量的个数,

    math?formula=%5Calpha_%7Bk%7D 为第k个高斯分量的权值。因为观测数据来自GMM的各个高斯分量相互独立,而

    math?formula=%5Calpha_%7Bk%7D 刚好可以看做是观测数据来自第k个高斯分量的概率,因此可以直接通过连乘得到整个隐变量

    math?formula=%5CGamma_%7Bj%7D的先验分布概率。

    ​得到完全数据的似然函数

    对于观测数据

    math?formula=y_%7Bi%7D,当已知其是哪个高斯分量生成的之后,其服从的概率分布为:

    math?formula=p(y_%7Bi%7D%7Cr_%7Bjk%7D%3D1%3B%5Ctheta)%3DN(y_%7Bj%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D) 由于观测数据从哪个高斯分量生成这个事件之间的相互独立的,因此可以写成:

    math?formula=p(y_%7Bi%7D%7C%5CGamma_%7Bj%7D%3B%5Ctheta)%3D%5Cprod_%7Bk%3D1%7D%5E%7BK%7DN(y_%7Bj%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D)%5E%7Br_%7Bjk%7D%7D

    这样我们就得到了已知

    math?formula=%5CGamma_%7Bj%7D的情况下单个测试数据的后验概率分布。结合之前得到的

    math?formula=%5CGamma_%7Bj%7D的先验分布,则我们可以写出单个完全观测数据的似然函数为:

    math?formula=p(y_%7Bj%7D%2C%5CGamma_%7Bj%7D%2C%5Ctheta)%3D%5Cprod_%7Bk%3D1%7D%5E%7BK%7D%5Calpha_%7Bk%7D%5E%7Br_%7Bjk%7D%7DN(y_%7Bj%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D)%5E%7Br_%7Bjk%7D%7D取对数,得到对数似然函数为:

    math?formula=lnp(y%2C%5CGamma_%7Bj%7D%3B%5Ctheta)%3D%5Csum_%7Bj%3D1%7D%5E%7BN%7D%5Csum_%7Bk%3D1%7D%5E%7BK%7D(r_%7Bjk%7D%5Cln%5Ctheta_%7Bk%7D%2Br_%7Bjk%7D%5Cln%20N(y_%7Bj%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D))得到各个高斯分量的参数计算公式,

    首先,我们将上式的

    math?formula=%5Cln%20N(y_%7Bj%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D)根据单高斯的向量形式的概率密度函数的表达形式展开:

    math?formula=%5Cln%20N(y_%7Bj%7D%7C%5Cmu_%7Bk%7D%2C%CE%A3_%7Bk%7D)%3D-%5Cfrac%7BD%7D%7B2%7D%5Cln%20(2%5Cpi)-%5Cfrac%7B1%7D%7B2%7D%5Cln%20%7C%CE%A3_%7Bk%7D%7C-%5Cfrac%7B1%7D%7B2%7D(y_%7Bj%7D-%5Cmu_%7Bk%7D)%5ET%CE%A3_%7Bk%7D%5E%7B-1%7D(y_%7Bj%7D-%5Cmu_%7Bk%7D)假设我们已经知道隐变量

    math?formula=r_%7Bjk%7D的取值,对上面得到的似然函数分别对

    math?formula=%5Calpha_%7Bk%7D

    math?formula=%CE%A3_%7Bk%7D求偏导并且偏导结果为零,可以得到:

    math?formula=%5Cmu_%7Bk%7D%3D%5Cfrac%7B%5Csum_%7Bj%3D1%7D%5E%7BN%7D%5Csum_%7Bk%3D1%7D%5E%7BK%7Dr_%7Bjk%7Dy_%7Bj%7D%7D%7B%5Csum_%7Bj%3D1%7D%5E%7BN%7D%5Csum_%7Bk%3D1%7D%5E%7BK%7Dr_%7Bjk%7D%7D

    math?formula=%CE%A3_%7Bk%7D%3D%5Cfrac%7B%5Csum_%7Bj%3D1%7D%5E%7BN%7D%5Csum_%7Bk%3D1%7D%5E%7BK%7Dr_%7Bjk%7D(y_%7Bj%7D-%5Cmu_%7Bk%7D)(y_%7Bj%7D-%5Cmu_%7Bk%7D)%5ET%7D%7B%5Csum_%7Bj%3D1%7D%5E%7BN%7D%5Csum_%7Bk%3D1%7D%5E%7BK%7Dr_%7Bjk%7D%7D现在参数空间中剩下一个

    math?formula=%5Calpha_%7Bk%7D还没有求。这时一个约束满足问题,因为必须满足约束

    math?formula=%CE%A3_%7Bk%7D%3D%5Csum_%7Bk%3D1%7D%5E%7BK%7D%5Calpha_%7Bk%7D%3D1.我们使用拉格朗日乘数法结合似然函数和约束条件对

    math?formula=%5Calpha_%7Bk%7D求偏导,可以得到:

    math?formula=%5Calpha_%7Bk%7D%3D%5Cfrac%7B%5Csum_%7Bj%3D1%7D%5E%7BN%7Dr_%7Bjk%7D%7D%7B-%5Clambda%7D将上式的左右两边分别对k=1,2,...,K求和,可以得到:

    math?formula=%5Clambda%3D-N

    math?formula=%5Clambda带入,最终得到:

    math?formula=%5Calpha_%7Bk%7D%3D%5Cfrac%7B%5Csum_%7Bj%3D1%7D%5E%7BN%7Dr_%7Bjk%7D%7D%7BN%7D至此,我们在隐变量已知的情况下得到了GMM的三种类型参数的求解公式。

    得到隐变量的估计公式

    根据EM算法,现在我们需要通过当前参数的取值得到隐变量的估计公式也就是说隐变量的期望的表现形式。即如何求解

    math?formula=E%7Br_%7Bjk%7D%7Cy%2C%5Ctheta%7D

    47c3d7ea3433

    四、案例分析

    案例一

    这里的数据集是每天经过美国西雅图弗莱蒙特桥上自行车的数量,时间间隔是[2012/10/02-2014/5/31]有需要的同学可以私信。

    1、先来看一下数据集

    import pandas as pd

    import matplotlib.pyplot as plt

    from sklearn.decomposition import PCA

    from sklearn.mixture import GaussianMixture

    def dataset():

    """

    读取数据并处理

    :return:

    """

    data = pd.read_csv('../../../数据集/机器学习/EM算法/FremontHourly.csv',index_col='Date',parse_dates=True)

    print(data.head(20))

    if __name__=="__main__":

    dataset()

    47c3d7ea3433

    2、将数据以图的形式展示

    data.plot()

    plt.show()

    47c3d7ea3433

    这时我们发现数据太多密集,我们以"周(week)"的形式展示

    #在时间上进行重采样

    data.resample('w').sum().plot()

    plt.show()

    47c3d7ea3433

    我们还可以以"天(day)"的形式将一年的数据呈现出来

    data.resample('D').sum().rolling(365).sum().plot()

    plt.show()

    47c3d7ea3433

    我们按照每天的时间进行求平均值,看一下一天的那个时间点是峰值

    data.groupby(data.index.time).mean().plot()

    plt.xticks(rotation=45)

    plt.show()

    47c3d7ea3433

    我们按照桥的两侧,分别计算每个时间点的停车量

    data.columns=['East','West']

    data['Total']=data['East']+data['West']

    pivoted = data.pivot_table('Total',index=data.index.time,columns=data.index.date)

    print(pivoted.iloc[:5,:5])

    pivoted.plot(legend=False,alpha=0.01)

    plt.xticks(rotation=45)

    plt.show()

    47c3d7ea3433

    47c3d7ea3433

    由于数据量比较少,上图不是很清晰,但是可以隐约看到是有两种分布的数据存在。

    接下来我们看一下数据结构,

    print(pivoted.shape)

    47c3d7ea3433

    通过上图我们得知数据是24 x 607的矩阵,这对于聚类来说是不合理的,数据量太少特征值过多,所以这里我们将矩阵进行转置,

    X = pivoted.fillna(0).T.values

    print(X.shape)

    47c3d7ea3433

    这时数据处理完后,我们发现特征值有点多,所以我们利用PCA进行降维处理,降到两维,我们看一下图像,

    #PCA降维成2维

    pca = PCA(n_components=2)

    Y = pca.fit_transform(X)

    # print(Y.shape)

    plt.scatter(Y[:,0],Y[:,1])

    plt.show()

    47c3d7ea3433

    降维结束后,我们开始使用GMM算法进行聚类,首先我们先看一下聚成两类之后的概率值,由于数据集比较明显,所以分出的概率值差异比较大,

    def gmm(data,pivoted):

    """

    利用GMM聚类算法进行聚类

    :return:

    """

    gmm = GaussianMixture(n_components=2)#两种高斯分布

    gmm.fit(data)

    labels = gmm.predict_proba(data)

    x = gmm.predict(data)

    print(labels)

    if __name__=="__main__":

    data,pivoted = dataset()

    gmm(data,pivoted)

    47c3d7ea3433

    我们输出一下分类的结果看一下,

    print(x)

    47c3d7ea3433

    接下来我们画图看一下分类的结果,

    plt.scatter(data[:,0],data[:,1],c=x,cmap='rainbow')

    plt.show()

    47c3d7ea3433

    最后我们看一下原始数据集分类后的分布图

    fig, ax = plt.subplots(1, 2, figsize=(14, 6))

    pivoted.T[labels == 0].T.plot(legend = False,alpha=0.1,ax=ax[0])

    pivoted.T[labels == 1].T.plot(legend=False, alpha=0.1, ax=ax[1])

    ax[0].set_title('Purple Cluster')

    ax[1].set_title('Red Cluster')

    plt.show()

    47c3d7ea3433

    这里是所有的代码

    import pandas as pd

    import matplotlib.pyplot as plt

    from sklearn.decomposition import PCA

    from sklearn.mixture import GaussianMixture

    def dataset():

    """

    读取数据并处理

    :return:

    """

    data = pd.read_csv('../../../数据集/机器学习/EM算法/FremontHourly.csv',index_col='Date',parse_dates=True)

    # print(data.head(30))

    # data.plot()

    # plt.show()

    # #在时间上进行重采样

    # data.resample('w').sum().plot()

    # data.resample('D').sum().rolling(365).sum().plot()

    # plt.show()

    # data.groupby(data.index.time).mean().plot()

    # plt.xticks(rotation=45)

    # plt.show()

    data.columns=['East','West']

    data['Total']=data['East']+data['West']

    pivoted = data.pivot_table('Total',index=data.index.time,columns=data.index.date)

    print(pivoted.iloc[:5,:5])

    pivoted.plot(legend=False,alpha=0.01)

    plt.xticks(rotation=45)

    plt.show()

    # print(pivoted.shape)

    X = pivoted.fillna(0).T.values

    print(X.shape)

    #PCA降维成2维

    pca = PCA(n_components=2)

    Y = pca.fit_transform(X)

    # print(Y.shape)

    plt.scatter(Y[:,0],Y[:,1])

    plt.show()

    return Y,pivoted

    def gmm(data,pivoted):

    """

    利用GMM聚类算法进行聚类

    :return:

    """

    gmm = GaussianMixture(n_components=2)#两种高斯分布

    gmm.fit(data)

    labels = gmm.predict_proba(data)

    x = gmm.predict(data)

    print(labels)

    plt.scatter(data[:,0],data[:,1],c=x,cmap='rainbow')

    plt.show()

    fig, ax = plt.subplots(1, 2, figsize=(14, 6))

    pivoted.T[labels == 0].T.plot(legend = False,alpha=0.1,ax=ax[0])

    pivoted.T[labels == 1].T.plot(legend=False, alpha=0.1, ax=ax[1])

    ax[0].set_title('Purple Cluster')

    ax[1].set_title('Red Cluster')

    plt.show()

    return None

    if __name__=="__main__":

    data,pivoted = dataset()

    gmm(data,pivoted)

    案例二

    这里我们举例对比一下GMM算法与K-means算法的聚类效果。

    首先我们先制作一些数据集,同时画图展示一下。

    from sklearn.datasets.samples_generator import make_blobs

    import matplotlib.pyplot as plt

    from sklearn.cluster import KMeans

    from sklearn.mixture import GaussianMixture

    import numpy as np

    def datasets():

    """

    使用datasets包产生一些数据

    :return:

    """

    plt.rcParams['axes.unicode_minus'] = False#解决不显示负数问题

    X,y_true = make_blobs(n_samples=800,centers=4,random_state=11)

    plt.scatter(X[:,0],X[:,1])

    plt.show()

    if __name__ == "__main__":

    datasets()

    47c3d7ea3433

    接下来我们先用K-means算法聚类一下,看看效果。

    def GmmKmean(data):

    """

    GMM算法与Kmeans算法对比

    :return:

    """

    kmeans = KMeans(n_clusters=4)

    kmeans.fit(data)

    y_kmeans = kmeans.predict(data)

    plt.scatter(data[:,0],data[:,1],c=y_kmeans,s=50,cmap='viridis')

    plt.show()

    centers = kmeans.cluster_centers_

    print(centers)

    return None

    if __name__ == "__main__":

    data = datasets()

    GmmKmean(data)

    47c3d7ea3433

    47c3d7ea3433

    我们再用GMM算法实现一下,

    gmm = GaussianMixture(n_components=4,random_state=1)

    gmm.fit(data)

    labels = gmm.predict(data)

    plt.scatter(data[:,0],data[:,1],c=labels,s=40,cmap='viridis')

    plt.show()

    47c3d7ea3433

    通过上面两个算法的对比,我们发现比较简单的数据集上,两者的分类效果差不多。

    接下来,我们再制作一些数据集对比一下聚类效果,

    plt.show()

    rng = np.random.RandomState(13)

    Y = np.dot(X,rng.randn(2,2))

    plt.scatter(Y[:,0],Y[:,1])

    plt.show()

    47c3d7ea3433

    我们的目的是希望能分出四个类,左上角两个,右下角两个

    首先我们先看一下K-means算法的效果,

    kmeansy = KMeans(n_clusters=4,random_state=1)

    kmeansy.fit(dataY)

    datay_kmeans = kmeansy.predict(dataY)

    plt.scatter(dataY[:,0],dataY[:,1],c=datay_kmeans,s=40,cmap='viridis')

    plt.show()

    47c3d7ea3433

    通过上图我们发现,K-means算法能够聚类出来,但是左上角两个类分出的效果并不是很理想,接下来我们用GMM算法再聚类一下,

    gmm.fit(dataY)

    labelsY = gmm.predict(dataY)

    plt.scatter(dataY[:,0],dataY[:,1],c=labelsY,s=40,cmap='viridis')

    plt.show()

    47c3d7ea3433

    这时我们发现GMM算法的聚类效果就非常好。

    代码在这里展示一下,

    from sklearn.datasets.samples_generator import make_blobs

    import matplotlib.pyplot as plt

    from sklearn.cluster import KMeans

    from sklearn.mixture import GaussianMixture

    import numpy as np

    def datasets():

    """

    使用datasets包产生一些数据

    :return:

    """

    plt.rcParams['axes.unicode_minus'] = False#解决不显示负数问题

    X,y_true = make_blobs(n_samples=800,centers=4,random_state=11)

    plt.scatter(X[:,0],X[:,1])

    plt.show()

    rng = np.random.RandomState(13)

    Y = np.dot(X,rng.randn(2,2))

    plt.scatter(Y[:,0],Y[:,1])

    plt.show()

    return X,Y

    def GmmKmean(data,dataY):

    """

    GMM算法与Kmeans算法对比

    :return:

    """

    kmeans = KMeans(n_clusters=4)

    kmeans.fit(data)

    y_kmeans = kmeans.predict(data)

    plt.scatter(data[:,0],data[:,1],c=y_kmeans,s=50,cmap='viridis')

    plt.show()

    centers = kmeans.cluster_centers_

    print(centers)

    gmm = GaussianMixture(n_components=4,random_state=1)

    gmm.fit(data)

    labels = gmm.predict(data)

    plt.scatter(data[:,0],data[:,1],c=labels,s=40,cmap='viridis')

    plt.show()

    kmeansy = KMeans(n_clusters=4,random_state=1)

    kmeansy.fit(dataY)

    datay_kmeans = kmeansy.predict(dataY)

    plt.scatter(dataY[:,0],dataY[:,1],c=datay_kmeans,s=40,cmap='viridis')

    plt.show()

    # gmmY = GaussianMixture(n_components=4)

    gmm.fit(dataY)

    labelsY = gmm.predict(dataY)

    plt.scatter(dataY[:,0],dataY[:,1],c=labelsY,s=40,cmap='viridis')

    plt.show()

    return None

    if __name__ == "__main__":

    data,dataY = datasets()

    GmmKmean(data,dataY)

    好了,EM算法与GMM算法的学习到这里就结束了,本节的公式推导比较多,希望大家能够手动推一下!

    展开全文
  • 什么是极大似然估计?参数估计就是通过若干次试验,已知某个参数能使这个样本出现的概率最大,...极大似然估计的原理?一个随机试验如有若干个可能的结果A,B,C,…。若在仅仅作一次试验中,结果A出现,则一般认为...
  • 卡尔曼滤波是什么 卡尔曼滤波适用于估计一个动态系统最优状态。即便是观测到系统状态参数含有噪声,观测值不准确,卡尔曼滤波也能够完成对状态真实值最优估计。网上大多数教程讲到卡尔曼数学公式推导,会...
  • Ordinary least squares是什么意思?

    千次阅读 2017-03-06 14:53:32
    应用最多的参数估计方法,也从最小二乘原理出发其他估计方法基础。 英语解释:In statistics, ordinary least squares (OLS) or linear least squares is a method for estimating the unknown ...
  • 1、什么是参数估计???? 参数估计:本质对未知参数作出估计。又分为点估计和区间估计两种类型。 设总体 X分布函数形式已知,但它一个或多个参数未知,借助于总体一个样本来估计总体未知参数...
  • 展开全部因为要提高e69da5e887aa62616964757a686964616f31333433656130参数估计的无偏性,两阶段最小二乘法用于检验有内生性变量的回归模型。工具变量法对于恰好识别的结构方程有效的。但对过度识别方程虽然能够给...
  • 1.什么是VIO估计初始化? 初始化就是把变量赋值为默认值,把控件设置为默认状态,把每准备好准备好。这么看来,初始化问题似乎并不那么困难,但不幸的是,把初始化这个定义用在整个系统层面来讲就不是那么...
  • 广义似然比原理

    千次阅读 2018-03-09 15:24:21
    广义似然比检验检验的是什么参数,分布正态分布: μ=μ0?或μ&lt;-[μ1, μ2]?概率有多大假设检验:构造统计量困难原理:根据极大似然估计构造统计量,-2lnλ服从卡方分布。似然比检验实质是在比较有约束...
  • ​ 极大似然估计建立在极大似然原理一种参数估计方法。其目的利用已知样本结果,反推最有可能导致这样结果参数值。 通俗地说,就是通过若干次试验,观察其结果,利用试验结果得到某个参数值能够使样本...
  • 极大似然估计-大白话

    2021-04-11 09:56:56
    极大似然估计-大白话预备知识——排列组合极大似然估计预备...(一种使用最广泛的参数估计方法) 例子: 参考链接已知实验结果,求小球概率。 问题: 这个函数为什么是凸函数,为什么求导=0,就可找到使这个函数最
  • 极大似然思想原理

    千次阅读 2015-08-14 19:39:00
    极大似然估计,顾名思义一种估计方法。既然一种估计方法,我们至少必须搞清楚几个问题:估计什么?需要什么前提或假设?...这一方法基于这样的思想:我们所估计的模型参数,要使得产生这个给定样本的可能性最大
  • 卡尔曼滤波原理及实现 前一段时间,做项目研究了一下卡尔曼滤波,并且在项目当中实现了一个物体跟踪功能,所以,借着新鲜...即便观测到系统状态参数含有噪声,观测值不准确,卡尔曼滤波也能够完成对状态真实...
  • 最大似然估计方法介绍

    千次阅读 2018-10-06 14:46:58
    然而面对繁琐数学公式和复杂推理过程,使我对概念非常模糊,也不懂到底是什么原理,但通过后来慢慢学习使我对最大似然估计有了不一样认识。 最大似然估计是使用概率模型、利用已知样本结果,反推一件事情...
  • 极大似然估计也可以用于估计,总结下极大似然估计的原理和方法,以及如何用于经典线性回归方程的估计。 首先,什么是极大似然估计? 虽然很多理论基于极大似然估计的,但是,要强调的,极大似然估计忽略了小...
  • 在机器学习中,尤其是回归模型,经常用到梯度下降法和最小二乘法,这里把最小二乘法的原理及代码实现总结处理。 1 最小二乘法原理 首先要清楚,最小二乘法要解决的是什么问题呢?根据前面的线性回归,我们知道线性...
  • 于是,Linux下可以使用emacs --deamon来启动Emacs作为一个守护进程,但该参数不支持windows平台(虽然说在windows平台使用Emacs一件感觉很别扭事情),估计是使用了windows平台没有特性,具体使用的什么特性我...
  • 梯度下降法原理

    2020-03-29 23:31:28
    (该博文为一网友所写,非常详细易懂,故搬运过来以后方便回忆学习) 一、为什么需要梯度下降法 ...在绝大多数情况下,损失函数很复杂(比如逻辑回归),根本无法得到参数估计表达式。因此需要一种对大...

空空如也

空空如也

1 2 3 4 5 ... 11
收藏数 210
精华内容 84
关键字:

参数估计的原理是什么