網誌搬家了 I moved to github!

我的網誌移到 github 上囉! 以後的新文章都會發佈在 https://yenhsunlin.github.io/

留言

張貼留言

這個網誌中的熱門文章

IA9a - K-means分群與EM算法: 理論

圖片
市場內分門別類的產品 原則上 K-means分群 是 EM算法 (expectation-maximization algorithm) 的一個特例，相較於EM從似然函數 (likelihood) 的方式出發並導出不同資料點對不同群 (或稱聚類、簇，cluster) 的貢獻度，K-means單純用歐氏距離來分類，雖然這加快了K-means算法的計算效率，但顯然的這麼naïve的方法會產生明顯的缺陷。但仍讓我們會先從簡單的K-means開始，並探討K-means算法的結果在某些情況下和人類肉眼辨識的結論會產生明顯的歧義，接著再推進到通用的EM算法，除了 python 的demo外，我也會證明在數學上如何從EM得到K-means的極限 (僅限高斯混合的情況)，而在這個part內，讓我們先從數學理論的角度來面。 一、K-means分群 今天如果我有一大坨資料共120個點如圖1所示，從資料的散佈趨勢來看，很明顯的分成兩群，那我們是不是可以有一個演算法來自動化將資料點如同我們臆測的一樣分成兩袋？一個最簡單的發法便是使用K-means分群，而這個大寫K表示群的數量，以圖1而言，我們猜測K=2，所以是個2-means分群問題。 圖1. 資料散佈圖 (普魯士藍) 與初始化的紅綠兩點 　所以為了達成分群，我們需要先初始化兩個點，這兩個點代表兩群的中心。事實上，初始化的位置不必要一定剛好在兩群的中心，可以是任意挑選的地方，但為了計算方便，我們會將起始化的兩點開始在靠近資料點的平均值附近，如同圖1中紅綠兩點。我們將這兩個初始化的點叫作\(p_k\)其中\(k=1,2\)表示第\(k\)群。另外將資料點的集合用\(\mathbf{D}=\{d_1,d_2\cdots d_n \cdots d_{120}\}\)表示，所以我們計算每個點\(d_n=(x_{d_n},y_{d_n})\)和\(p_k=(x_{p_k},y_{p_k})\)的歐氏距離 \[ \ell_{n,k}=\sqrt{(x_{d_n}-x_{p_k})^2+(y_{d_n}-y_{p_k})^2}. \tag{1} \] 若該點和\(p_1\)的距離小於和\(p_2\)的距離 (\(\ell_{n,1}<\ell_{n,2}\))，則我們說\(d_n\)比較靠近\(p_1\)且應該...
閱讀完整內容

IA1 - 布豐投針，圓周率\(\pi\)的估計

圖片
Mathematica®生成400根針的隨機投放示意圖 布豐投針問題 (Buffon's needle problem) 是這樣敘述的： 有一房間的地板由數片平行並列的木板所構成，板與板之間溝紋的距離皆相同，如上圖棕色隔線的間距。試問隨機扔一針，該針觸及溝紋的機率為何？ 這項問題主要由18世紀的布豐伯爵所提出的，並由義大利數學家Mario Lazzarini在1901年實地進行拋針實驗，共拋出了3,408根針(🥴)並計算出\(\pi\)的近似值約為\(355/113\)，參見 維基百科 。在我們繼續追本溯源這個歷史問題以前，讓我們先看看當代最常見估計圓周率的方法。 一、灑點估計法 這個方法概念非常的簡單，我們將產生\(N\)個座標\((x,y)\)隨機散佈在\(0\)和\(1\)之間的點，所這些點將被包在以原點為中心，邊長為2的正方形的第一象限中 (first quarter)。接著我們計算出每個點距原點的長度\(d=\sqrt{x^2+y^2}\)，若\(d<1\)則紀錄，反之則丟棄，最後我們紀錄了\(n\)個點，這些點都被包在以原點為中心，半徑為\(1\)的\(1/4\)圓內，見圖1。 圖1. 在第一象限中藍色正方形的面積為\(1\)而紅色區域的面積為圓的\(1/4\) 　現在我們相信這些處於單位圓內點的數量和全部點的數量的比值應等於那\(1/4\)圓面積和正方形面積之比，即 \[ \frac{n}{N}=\frac{A_圓}{A_方}.\tag{1} \] 而\(A_正 = 1^2\)和\(A_圓=\pi 1^2/4\)，所以我們便得到了 \[ \pi = 4\frac{n}{N},\tag{2} \] 正是\(\pi\)的估計。我們預測若\(N\)愈大，則其估計值愈接近真正的圓周率。另外\(\pi\)可以解析表示成 \[ \pi = \int_{-1}^{1}\frac{dx}{\sqrt{1-x^2}}. \tag{3} \] 上式是由現代分析學之父卡爾·魏爾斯特拉斯 (Karl T. W. Weierstraß, 1815-1897) 在1841年所定義的。 　我們接著可以透過 rand 來實現這個簡單的蒙地卡羅模擬，腳本如下： import numpy as np import matp...
閱讀完整內容

IA9b - K-means分群與EM算法: 實作

圖片
這一篇文章延續了上一篇IA9a的內容，但側重在 python 的實作，所以所有的節數、圖案以及公式編號都會依照IA9a的編號接序下來，所以重複的內容就直接引用，不在這裡重新贅述。 四、實作EM 4.1 高斯混合分佈 讓我們從1.1節的協方差矩陣開始，生成的圖形就如圖4 (腳本內加了seed random方便我們每次追蹤)，另外也看到了執行5次K-means算法後 (再多次也不會變了) 得到的分群結果如圖5，這顯然與我們的預期有極大的落差！ 　但我們將按照3.3結給出EM算法的偽代碼，假設資料群是高斯混合分佈的形式，重新執行一遍分群，我們將看到EM算法能很好的分出我們預期的兩群的結果。我們預計會有兩群，所以\(k=1,2\)，且第一步需要初始化參數組\(\boldsymbol{\theta}=\{\mu_k,\Sigma_k,\pi_k\}\)。初始化的過程是給定一個起始點，原則上隨機選擇的方式是合理的，不過為了計算上方便，可以讓初始的協方差矩陣都為單位矩陣，\(\pi_k\)則為每個群均分 (這裡就是每個都是0.5)，如下代碼 # Initializing parameters sigma1,sigma2 = np.eye(2),np.eye(2) pi1 = 0.5 pi2 = 1-pi1 N = len(cluster) 完成初始的參數選擇後，我們用偽代碼的第二步，也就是E-step計算置信度\(\gamma (z_{nk})\)。這裡我將Eq. (10)列出來方便解說 \[ \gamma (z_{nk})=\frac{\pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}{\sum_{k=1}^2 \pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}\tag{10}. \] 因為我們是假設兩群都是高斯分佈，所以對第一個群而言的分子就是 # Starting E-step: Calculating responsibilities # Numerator for k=1 k1 = pi1*multivariate_normal.pdf(cluster,p1,sigma1) 直接透過 python 腳本我們很容易理解了分子的意義便是計算每個資料點分別對應到\(\{...
閱讀完整內容