2020年3月15日 星期日

Kernel Density Estimation 的入門篇

Kernel Density Estimation 的入門篇

對於像我這種統計基礎很虛浮的人來說,第一次遇到 Kernel Density Estimation 這樣的詞彙難免一陣混亂,這次就從機率角度出發,一路旅行到 Kernel Density Estimation,順便重整一下機率的基礎

Definition of Probability

從古典機率來探討一下機率的定義

先用n(A) 代表某個集合 A 的數量,後定義 一個事件 E 和樣本空間 S

機率基本定義就出來了

\[ P(E) = \frac{n(E)}{n(S)} \]

Cumulative Distribution Function

不管是骰子也好硬幣也好,只要當成事件來看,都可以用基本的定義公式算出個機率來,骰子定義6個事件,硬幣定義2個事件就可以。 但是當我們想要探討連續型的機率時,這個定義就很不夠用了。此時就要使用後來的發明,也就是累積機率分佈。

如果有一個隨機變數可以是0 到1的任何一個有理數,那事件不管怎麼定義都定義不完的。此時就可以採用一個累積的概念,像是

  • 數字在 0 ~ 0.3 之間的機率
  • 數字在 0.3 ~ 0.6 之間的機率
  • 數字在 0.6 ~ 1 之間的機率

當然這個定義仍然是怪怪的,會直接重複到數字為 0.3 和 0.6 的機率

為了更好解釋,將這個隨機變數定義為\(X\),同時引入區間的概念

  • 數字在 \(0 \leq X \leq 0.3\) 之間的機率
  • 數字在 \(0.3 < X \leq 0.6\)之間的機率
  • 數字在 \(0.6 < X \leq 1\) 之間的機率

有了區間後,做一些數學符號上的簡化

  • 數字在 \([0 \quad ; 0.3]\) 之間的機率
  • 數字在 \((0.3 \quad ; 0.6]\)之間的機率
  • 數字在 \((0.6 \quad ; 1]\) 之間的機率

準備好了就可以定義 CDF 的公式 \[ F_X(x) = P(X \leq x) \]

這個\(F_X(x)\) 讓一些 原本需要定義一堆事件的機率可以用運算來表示出來

數字在 \([0 \quad ; 0.3]\) 之間的機率 =\(F_X(0.3)\)

數字在 \((0.3 \quad ; 0.6]\)之間的機率 =\(F_X(0.6) - F_X(0.3)\)

數字在 \((0.6 \quad ; 1]\) 之間的機率 =\(F_X(1) - F_X(0.6)\)

probability density function

即使CDF 讓機率可以被計算出來,但是想要定義某一個區間仍然要做一番加減,後面就出現了可以讓機率更單純的方式,就是機率密度函數PDF

關於PDF,用畫圖會更容易解釋

下面這張圖,每個柱體的高度並不是機率,那是什麼呢?

對於機率密度函數來說,每一塊 機率密度的關鍵就在於,總面積為1,因此在每個組距為0.2 的情況下計算柱體的面積時

\[ \begin{aligned} 0.2 \times 0.75 &\quad= 0.15 \\ 0.2 \times 0.75 &\quad= 0.15 \\ 0.2 \times 1 &\quad= 0.2 \\ 0.2 \times 1.45 &\quad= 0.29 \\ 0.2 \times 1.05 &\quad= 0.21 \\ \end{aligned} \]

\[0.15 + 0.15 + 0.2 + 0.29+ 0.21 = 1\]

即使將組距改成0.1,相當於每個柱體又分成兩個小的柱體

第1個大柱體分割的小柱體

  • \(0.1 \times 0.9 = 0.9\)
  • \(0.1 \times 0.6 = 0.6\)
  • \(0.9 + 0.6 = 0.15\)

第2個大柱體分割的小柱體

  • \(0.1 \times 0.7 = 0.7\)
  • \(0.1 \times 0.8 = 0.8\)
  • \(0.7 + 0.8 = 0.15\)

這就是機率密度函數的意義

Kernel Density Estimation

PDF 其實是幫助我們更了解連續型隨機變數的情況下任意值出現的機率是多少

上面的資料就用了常態分佈來擬合一個最接近的機率分佈 \[ pr(x | \mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-(x-\mu)^2/2\sigma^2} \]

常態分佈是統計的一個利器,每一個剛接觸常態分佈的人都會問到一個問題:

如果資料不是常態分佈該怎麼辦呢?

因此KDE 採用的是不假設機率分佈的方法,即使不知道資料的分佈,一樣可以擬出一個曲線 \[ P(x;h) = \frac{1}{nh} \sum_{i=1}^n K\Big( \frac{x-x_i}{h} \Big) \]

  • h 代表 bandwidth ,也就是帶寬
  • K 代表 Kernel
  • n 代表資料的總筆數
  • x 代表隨機變數對應到數線上的值

不管是怎樣的分佈,KDE在設定 kernel 和 bandwidth 下都可以擬出一條分佈

下面用常態分布先產生一堆亂數,然後再分別用常態分佈和KDE 畫出曲線

set.seed(1234567)
samp <- rnorm(n = 100, mean = 50, sd = 1)
x <- seq(0,100,by=1)

hist(samp,probability = TRUE,breaks=10)
curve(dnorm(x,mean=50,sd=1),from=0,to=100,n=1000,col=2,add = T)

hist(samp,probability = TRUE,breaks=10)
lines(density(samp,bw=.5),col='red')

可以看到曲線非常的接近

最後來看一下bandwidth 對KDE 的影響

hist(samp,probability = TRUE,breaks=10)
lines(density(samp,bw=.1),col='darkgreen',lwd=2)
lines(density(samp,bw=.5),col='red',lwd=2)
lines(density(samp,bw=.9),col='blue',lwd=2)
legend("topleft", 
  legend = c("bw=0.1", "bw=0.5","bw=0.9"), 
  col = c('darkgreen','red','blue'),lty=1)

沒有留言:

張貼留言

使用LibreOffice製作代碼瀑布

Matrix_Digital_Rain_in_Excel_1 Matrix Digital Rain in LibreOffice 之前在facebook 看到一個沒有見過的效果,一時興起就研究了一下 如何使用LibreOffie...