2. 程式人生 > >em算法系列1：python實現三硬幣模型

em算法系列1：python實現三硬幣模型

阿新 • • 發佈：2018-12-13

三硬幣模型

假設有3枚硬幣，分別記做A，B，C。這些硬幣正面出現的概率分別是π,p和q。進行如下擲硬幣實驗：先擲硬幣A，根據其結果選出硬幣B或C，正面選B，反面選硬幣C；然後投擲選重中的硬幣，出現正面記作1，反面記作0；獨立地重複n次（n=10)，結果為1111110000 我們只能觀察投擲硬幣的結果，而不知其過程，估計這三個引數π,p和q。

EM演算法

可以看到投擲硬幣時到底選擇了B或者C是未知的。我們設隱藏變數Z 來指示來自於哪個硬幣，Z={z1,z2,…,zn}Z=\{z_1,z_2,…,z_n\}，令θ={π,p,q}θ=\{π,p,q\}，觀察資料X

=x1,x2,…,xnX={x_1,x_2,…,x_n}。

寫出生成一個硬幣時的概率： P(x∣θ)=∑zP(x,z∣θ)=∑zP(z∣π)P(x∣z,θ)=πpx(1−p)1−x+(1−π)qx(1−q)1−x\begin{aligned}P(x|θ) &=∑zP(x,z|θ) = ∑zP(z|π)P(x|z,θ) \\ &= πp^x(1−p)^{1−x}+(1−π)q^x(1−q)^{1−x} \end{aligned}

x(1−q)1−x​ 有了一個硬幣的概率，我們就可以寫出所有觀察資料的log似然函式： L(θ∣X)=logP(X∣θ)=∑j=1nlog[πpxj(1−p)1−xj+(1−π)qxj(1−q)1−xj]\begin{aligned} L(θ|X)=logP(X|θ)=∑_{j=1}^nlog[πp^{x_j}(1−p)^{1−x_j} + (1−π)q^{x_j}(1−q)^{1−x_j}] \end{aligned} 然後求似然函式最大值 θ∗=argmaxL(θ
∣X) θ^*= arg\quad maxL(θ|X) 其中L(θ∣X)=logP(X∣θ)=log∑ZP(X,Z∣θ)L(θ|X)=logP(X|θ)=log∑_ZP(X,Z|θ)。因為loglog裡面帶著加和所以這個極大似然是求不出解析解的。 如果我們知道隱藏變數Z的話，求以下的似然會容易很多: L(θ|X,Z)=logP(X,Z|θ) 但是隱藏變數Z的值誰也不知道，所以我們轉用EM演算法求它的後驗概率期望EZ∣X,θ(L(θ∣X,Z))E_{Z|X,θ}(L(θ|X,Z))的最大值。

  • E步：：假設當前模型的引數為π,p,qπ,p,q時，隱含變數來自於硬幣B的後驗概率 μ=P(Z∣X,θ)=πpxi(1−p)1−xiπpxi(1−p)1−xi+(1−π)qxi(1−q)1−xi\begin{aligned} μ=P(Z|X,θ)= \frac{πp^{x_i}(1−p)^{1−x_i}}{πp^{x_i}(1−p)^{1−x_i}+(1−π)q^{x_i}(1−q)^{1−x_i}} \end{aligned} 那麼隱含變數來自於硬幣C的後驗概率自然為1−μ1−μ。
def cal_u(pi1,p1,q1,xi):
    return pi1*math.pow(p1,xi) * math.pow(1-p1,1-xi)/\
           float(pi1* math.pow(p1,xi) * math.pow(1-p1,1-xi)+
                 (1-pi1) * math.pow(q1,xi) * math.pow(1-q1,1-xi))
  • M步： 先寫出似然關於後驗概率的期望，它是似然期望下界函式的最大值， EZ∣X,θ(L(θ∣X,Z))=∑j=1n[μlogπpxi(1−p)1−xiμ+(1−μ)log(1−π)qxi(1−q)1−xi1−μ]E_{Z|X,θ}(L(θ|X,Z))=∑_{j=1}^n [μlog \frac{πp^{x_i}(1−p)^{1−x_i}}{μ} + (1−μ)\frac{log(1−π)q^{x_i}(1−q)^{1−x_i}}{1−μ}] 要注意這裡把μμ看做固定值。然後我們分別求偏導，獲得引數π,p,qπ,p,q的新估計值 π=1n∑j=1nμjp=∑j=1nμjxj∑j=1nμjq=∑j=1n(1−μj)xj∑j=1n(1−μj) \begin{aligned} π &=\frac{1}{n} ∑_{j=1}^nμ_j \\ p&=\frac{∑_{j=1}^nμ_jx_j}{∑_{j=1}^nμ_j }\\ q&=\frac{∑^n_{j=1}(1−μ_j)xj}{∑^n_{j=1}(1−μ_j)} \end{aligned} 
def cal_u(pi, p, q, xi):
    """
      u值計算
    :param pi: 下一次迭代開始的 pi
    :param p:  下一次迭代開始的 p
    :param q:  下一次迭代開始的 q
    :param xi: 觀察資料第i個值，從0開始
    :return: 
    """
    return pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) / \
           float(pi * math.pow(p, xi) * math.pow(1 - p, 1 - xi) +
                 (1 - pi) * math.pow(q, xi) * math.pow(1 - q, 1 - xi))

def e_step(pi,p,q,x):
    """
        e步計算
    :param pi: 下一次迭代開始的 pi
    :param p:  下一次迭代開始的 p
    :param q:  下一次迭代開始的 q
    :param x: 觀察資料
    :return:
    """
    return [cal_u(pi,p,q,xi) for xi in x]

演算法首先選取引數初始值θ(0)=π(0),p(0),q(0)θ(0)=π(0),p(0),q(0)，然後迭代到收斂為止

python實現

# -*- coding: utf-8 -*-
import math

def cal_u(pi, p, q 
 

            
  

           

              
              

            

            
相關推薦

			   
            
            
            
 

    


    
    
em算法系列1：python實現三硬幣模型

     
 
							
							
							三硬幣模型
假設有3枚硬幣，分別記做A，B，C。這些硬幣正面出現的概率分別是π,p和q。進行如下擲硬幣實驗：先擲硬幣A，根據其結果選出硬幣B或C，正面選B，反面選硬幣C；然後投擲選重中的硬幣，出現正面記作1，反面記作0；獨立地重複n次（n=10)，結果為1111 



  
 

    


    
    
算法系列- 棧：如何實現瀏覽器的前進和後退功能？

     
  
 
 
 整理自極客時間-資料結構與演算法之美。原文內容更完整具體，且有音訊。購買地址： 
  
 1.如何理解棧 
 關於“棧”，有一個非常貼切的例子，就是一摞疊在一起的盤子。我們平時放盤子的時候，都是從下往上一個一個放；取的時候，我們也是從上往下一個一個地依次取，不能從中間任意抽出。後進先出，先進後出 



  
 

    


    
    
機器學習經典算法具體解釋及Python實現--線性回歸（Linear Regression）算法

     
 ica   single   方便   最好的   而且   ==   show   des   fun   

（一）認識回歸
回歸是統計學中最有力的工具之中的一個。機器學習監督學習算法分為分類算法和回歸算法兩種，事實上就是依據類別標簽分布類型為離散型、連續性而定義的。顧名思義。分類算法用於離散型分布 



  
 

    


    
    
常見的查找算法的原理及python實現

     
 put   arch   img   字典   python實現   需要   技術   ()   one   
順序查找
二分查找
練習


一、順序查找

data=[1,3,4,5,6]
value=1
def linear_search(data,value):
    flag=False
     



  
 

    


    
    
梅森算法生成隨機數的Python實現

     
 urn   func   __init__   object   set   ret   function   obj   int   import time
class Util(object):    def __init__(self):        self.index = 624        s 



  
 

    


    
    
資料結構與算法系列1--簡介

     
  
  
  
 為什麼要學習資料結構和演算法？ 
 1.直接好處是能夠有寫出效能更優的程式碼 2.演算法，是一種解決問題的思路和方法，有機會應用到生活和事業的其他方面。 3.長期來看，大腦思考能力是個人最重要的核心競爭力，而演算法是為數不多的能夠有效訓練大腦思考能力的途徑之一。 
 什麼是資料結構？什麼是演 



  
 

    


    
    
白話經典算法系列之一 氣泡排序的三種實現

     
                  氣泡排序是非常容易理解和實現，，以從小到大排序舉例：設陣列長度為N。1．比較相鄰的前後二個數據，如果前面資料大於後面的資料，就將二個數據交換。2．這樣對陣列的第0個數據到N-1個數據進行一次遍歷後，最大的一個數據就“沉”到陣列第N-1個位置。3．N=N-1，如果N不為0就重複 



  
 

    


    
    
算法系列-佇列：佇列線上程池等有限資源池中的應用

     
  
 
 
 整理自極客時間-資料結構與演算法之美。原文內容更完整具體，且有音訊。購買地址： 
  
 1.如何理解佇列 
 佇列可以把它想象成排隊買票，先來的先買，後來的人只能站末尾，不允許插隊。先進者先出，這就是典型的“佇列”。 
 佇列跟棧非常相似，支援的操作也很有限，最基本的操作也是兩個：入隊 enq 



  
 

    


    
    
算法系列——陣列：為什麼很多程式語言中陣列都從0開始編號？

     
  
 
 
 整理自極客時間-資料結構與演算法之美。原文內容更完整具體，且有音訊。購買地址： 
  
 1.什麼是陣列 
 陣列（Array）是一種線性表資料結構。它用一組連續的記憶體空間，來儲存一組具有相同型別的資料。 
 這個定義裡有幾個關鍵詞，理解了這幾個關鍵詞，我想你就能徹底掌握陣列的概念了。 
 
 



  
 

    


    
    
[閉目洞察算法系列之二]插入排序三種實現

     
 
                
直接插入排序的思想是：將一個待排序的數字按照指定的排序規則，插入到一個前面已經排好序的的子序列中， 從而形成一個新的子序列，等待下一個數字的插入， 直至所有數字排序結束。


假設有陣列a[0 ... n - 1], 排序流程為：
a[0]位第一個有序子序列， a[1 .. 



  
 

    


    
    
C語言實現TOP K算法系列之快速排序實現

     
 
								
								            
							
							
							TOP K演算法的實現有很多種方式，其中類似於快排的實現是非常棒的，堆的實現也是非常好的，其中就是關於快排的實現會得到一個TOP K的集合，而這個集合不一定保證裡面的資料都是有序的。

下面就獻上TOP 



  
 

    


    
    
練習程式：算法系列14：狼、羊、菜和農夫過河問題

     
 
							
							
							描述、思路





補充知識

參考2：函式物件：關於仿函式、函式物件、ptr_fun 
參考3：bind2nd使用

實現1



#include <iostream>
#include <deque>
#include <s 



  
 

    


    
    
算法系列之一 ：Google方程式

     
 
                算法系列之一 ： Google方程式    有一個字元組成的等式：WWWDOT - GOOGLE = DOTCOM，每個字元代表一個0－9之間的數字，WWWDOT、GOOGLE和DOTCOM都是合法的數字，不能以0開頭。請找出一組字元和數字的對應關係，使它們互相替換，並且替換 



  
 

    


    
    
算法系列1 初識演算法 演算法複雜性模型 演算法複雜度的計算

     
 ## 算法系列1 初識演算法

## 什麼是演算法？
![在這裡插入圖片描述](https://img-blog.csdnimg.cn/20200919091641177.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10, 



  
 

    


    
    
機器學習：Python實現聚類算法(三)之總結

     
 .fig   ask   class   ted   ssi   缺點   處理   blob   ron         考慮到學習知識的順序及效率問題，所以後續的幾種聚類方法不再詳細講解原理，也不再寫python實現的源代碼，只介紹下算法的基本思路，使大家對每種算法有個直觀的印象，從而可以更好的理解函數中 



  
 

    


    
    
python：遞歸算法求和1-100

     
 true   In   pen   數字   IE   return   app   遞歸算法   clas   
#遞歸算法求和1-100def qiehe():    def he(lis):        if lis == []:            return 0        return l 



  
 

    


    
    
算法系列-連結串列：如何實現LRU快取淘汰演算法

     
 
                整理自極客時間-資料結構與演算法之美。原文內容更完整具體，且有音訊。購買地址：



開篇語

今天我們來聊聊“連結串列（Linked list）經典的應用場景，那就是 LRU 快取淘汰演算法。

快取是一種提高資料讀取效能的技術，在硬體設計、軟體開發中都有著非常廣泛的應用， 



  
 

    


    
    
算法系列（三）插入排序的兩種改進：規避邊界檢測和取消交換（Java實現）

     
 
                
前言：演算法第四版習題2.1.24插入排序的哨兵和習題2.1.25不需要交換的插入排序
規避邊界檢測：
在插入排序的實現中先找到最小的元素並將其置於陣列的第一個位置，可以省掉內迴圈的判斷條件 j>0 。能夠省略判斷條件的元素稱為哨兵。

public class Ex 



  
 

    


    
    
運籌系列1：線性規劃單純形法python程式碼

     
 
							
							
							1. 模型

常見的線性規劃模型如下： 
max z=cxz=cx 
s.t. Ax=bAx=b



2. 求解步驟

假設B是基變數集合，通過矩陣的線性變換，基變數可由非基變量表示： 
x′i=ci+Σj∉Bmi,jx′j,(i∈B)xi′=ci+Σj∉Bm 



  
 

    


    
    
算法系列-動態規劃(1)：初識動態規劃

     
 昨天，羅拉去面試回來，垂頭喪氣。顯然是面試不順利，我趕忙過去安慰。

經過詢問才知道，羅拉麵試掛在了動態規劃。

說到動態規劃，八哥可就來精神了，於是就結合勞拉的面試題簡單的和她介紹了動態規劃。

事情是這樣的，勞拉的面試官給了她一道題，題目如下：  

```
有一個數列，規律如下：1、1、2、3、5、8、