線性判別分析

本文的閱讀等級:中級

在機器學習和模式識別中,分類 (classication) 可視為一種決策問題:給定一數據點,判斷它所屬的類別。本文介紹源自於統計學多變量分析的一個古典分類法,稱作線性判別分析 (linear discriminant analysis,簡稱 LDA)。就理論面來說,線性判別分析與費雪 (Ronald Fisher) 的判別分析 (一種應用於分類問題的降維方法,見“費雪的判別分析與線性判別分析”) 和邏輯斯回歸 (logistic regression,一種應用於分類問題的非線性模型) 有著密切的關係。就應用面而言,由於線性判別分析建立於嚴苛的假設上,它的分類效能並不突出,或許因為這個緣故,線性判別分析經常被當作與其他方法比較的基準。

 
假設給定一筆維數等於 d,樣本大小為 n,包含 k\ge 2 個類別的數據 \mathcal{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_n\}。數據點 \mathbf{x}_1,\ldots,\mathbf{x}_n 散布在 \mathbb{R}^d 空間,我們以 \mathcal{C}_1,\ldots,\mathcal{C}_k 標記類別或代表類別的指標集,例如,i\in\mathcal{C}_j 表示 \mathbf{x}_i 來自第 j 類,或直接說 \mathbf{x}_i 屬於 \mathcal{C}_j。給定一數據點 \mathbf{x},如何判定它應歸於何類?從機率 (概率) 學的觀點,貝氏定理 (Bayes’ theorem) 提供了解答 (見“貝氏定理──量化思考的利器”):

\displaystyle  P(\mathcal{C}_j\vert\mathbf{x})=\frac{p(\mathbf{x}\vert\mathcal{C}_j)P(\mathcal{C}_j)}{p(\mathbf{x})},~~~j=1,\ldots,k

其中

  1. P(\mathcal{C}_j) 是類別 \mathcal{C}_j 出現的機率,稱為先驗機率 (priori probability);
  2. p(\mathbf{x}\vert\mathcal{C}_j) 是條件密度函數,即給定類別 \mathcal{C}_j,數據點 \mathbf{x} 的機率密度函數,也稱為似然 (likelihood);
  3. p(\mathbf{x}) 是數據點 \mathbf{x} 的機率密度函數,稱為證據 (evidence),算式為

    p(\mathbf{x})=p(\mathbf{x}\vert\mathcal{C}_1)P(\mathcal{C}_1)+\cdots+p(\mathbf{x}\vert\mathcal{C}_k)P(\mathcal{C}_k)

  4. P(\mathcal{C}_j\vert\mathbf{x}) 是指在給定數據點 \mathbf{x} 的情況下,該點屬於 \mathcal{C}_j 的機率,稱為後驗機率 (posterior probability)。

貝氏最佳分類器判定數據點 \mathbf{x} 給具有最大後驗機率的類別。具體地說,若 P(\mathcal{C}_l\vert\mathbf{x})=\max_{1\le j\le k}P(\mathcal{C}_j\vert\mathbf{x}),即 l=\arg\max_{1\le j\le k}P(\mathcal{C}_j\vert\mathbf{x}),則 \mathbf{x} 歸於 \mathcal{C}_l。我們的目標要利用給定的樣本估計出正確的後驗機率,以期得到一個表現良好的分類器。

 
S型函數與 Softmax 函數

先考慮二類別問題 (k=2)。類別 \mathcal{C}_1 的後驗機率可以寫成

\displaystyle\begin{aligned}  P(\mathcal{C}_1\vert\mathbf{x})&=\frac{p(\mathbf{x}\vert\mathcal{C}_1)P(\mathcal{C}_1)}{p(\mathbf{x})}\\  &=\frac{p(\mathbf{x}\vert\mathcal{C}_1)P(\mathcal{C}_1)}{p(\mathbf{x}\vert\mathcal{C}_1)P(\mathcal{C}_1)+p(\mathbf{x}\vert\mathcal{C}_2)P(\mathcal{C}_2)}\\  &=\frac{1}{1+\frac{p(\mathbf{x}\vert\mathcal{C}_2)P(\mathcal{C}_2)}{p(\mathbf{x}\vert\mathcal{C}_1)P(\mathcal{C}_1)}}\\  &=\frac{1}{1+\exp(-a)}=f(a),  \end{aligned}

其中

\displaystyle  a=\ln\frac{p(\mathbf{x}\vert\mathcal{C}_1)P(\mathcal{C}_1)}{p(\mathbf{x}\vert\mathcal{C}_2)P(\mathcal{C}_2)}=\ln\frac{P(\mathcal{C}_1\vert\mathbf{x})}{P(\mathcal{C}_2\vert\mathbf{x})}

稱為對數賠率 (log odds),且

\displaystyle  f(a)=\frac{1}{1+\exp(-a)}

稱為S型函數 (sigmoid function)。下圖顯示S型函數的型態,它將整個實數軸映射至區間 [0,1],因此也稱為擠壓函數 (squashing function),具有下列性質:

  • 對稱性:f(-a)=1-f(a)
  • 逆函數:a=\ln(\frac{f}{1-f}),稱為 logit 函數。
  • 導數:\frac{df}{da}=f(1-f)
Sigmoid function

Sigmoid function

 
考慮類別數大於2的情況,即 k>2,後驗機率為

\displaystyle  P(\mathcal{C}_j\vert\mathbf{x})=\frac{p(\mathbf{x}\vert\mathcal{C}_j)P(\mathcal{C}_j)}{p(\mathbf{x})}=\frac{p(\mathbf{x}\vert\mathcal{C}_j)P(\mathcal{C}_j)}{\sum_{l=1}^kp(\mathbf{x}\vert\mathcal{C}_l)P(\mathcal{C}_l)}

多類別的後驗機率和二類別的後驗機率有不同的形式,我們改採下列定義:

\displaystyle  a_j=\ln\left(p(\mathbf{x}\vert\mathcal{C}_j)P(\mathcal{C}_j)\right),~~1\le j\le k

所以,

\displaystyle  P(\mathcal{C}_j\vert\mathbf{x})=\frac{\exp(a_j)}{\sum_{l=1}^k\exp(a_l)},~~1\le j\le k

上式中,歸一化的指數函數也稱為 Softmax 函數,它是 max 函數的一種平滑版本,因為當所有 i\neq ja_j\gg a_i 可推得 P(\mathcal{C}_j\vert\mathbf{x})\approx 1P(\mathcal{C}_i\vert\mathbf{x})\approx 0。請注意,以 Softmax 函數表達的後驗機率並不唯一,因為數組 (a_1+c,\ldots,a_k+c) 的 Softmax 函數有相同的回傳值。

 
線性判別分析

上述推演的主旨在將後驗機率轉換以S型函數 (二類別) 或 Softmax 函數 (多類別) 表示,表面上,這個舉措似乎無足輕重。但S型函數和指數函數都是單調遞增函數,我們可以用 a(\mathbf{x})a_j(\mathbf{x}) 來取代後驗機率 (因為我們僅在乎後驗機率的大小排序)。如果 a(\mathbf{x})a_j(\mathbf{x}) 具有簡單的表達式,那麼上述做法尤其顯得有意義。欲求出後驗機率,不論是二類別的S型函數,還是多類別的 Softmax 函數,我們都必須得到條件密度函數 p(\mathbf{x}\vert\mathcal{C}_j)。使用機率生成模型法 (generative model),假設 p(\mathbf{x}\vert\mathcal{C}_j) 服從常態分布,而且所有的類別都有相同的共變異數矩陣 (covariance matrix),如下 (見“共變異數矩陣與常態分布”):

\displaystyle  p(\mathbf{x}\vert\mathcal{C}_j)=\frac{1}{(2\pi)^{n/2}\vert\Sigma\vert^{1/2}}\exp\left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_j)^T\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu}_j)\right\}

其中 \boldsymbol{\mu}_j\in\mathbb{R}^d 是類別 \mathcal{C}_j 的平均數向量,\Sigmad\times d 階共變異數矩陣,\vert\Sigma\vert 表示 \Sigma 的行列式。下面分別討論 k=2k>2 的後驗機率表達式。

 
對於 k=2

\displaystyle\begin{aligned}  a(\mathbf{x})&=\ln p(\mathbf{x}\vert\mathcal{C}_1)+\ln P(\mathcal{C}_1)-\ln p(\mathbf{x}\vert\mathcal{C}_2)-\ln P(\mathcal{C}_2)\\  &=-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_1)^T\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu}_1)+\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_2)^T\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu}_2)+\ln P(\mathcal{C}_1)-\ln P(\mathcal{C}_2)\\  &=(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)^T\Sigma^{-1}\mathbf{x}-\frac{1}{2}\boldsymbol{\mu}_1^T\Sigma^{-1}\boldsymbol{\mu}_1+\frac{1}{2}\boldsymbol{\mu}_2^T\Sigma^{-1}\boldsymbol{\mu}_2+\ln P(\mathcal{C}_1)-\ln P(\mathcal{C}_2)\\  &=\mathbf{w}^T\mathbf{x}+w_0,  \end{aligned}

其中

\displaystyle\begin{aligned}  \mathbf{w}&=\Sigma^{-1}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)\\  w_0&=-\frac{1}{2}\boldsymbol{\mu}_1^T\Sigma^{-1}\boldsymbol{\mu}_1+\frac{1}{2}\boldsymbol{\mu}_2^T\Sigma^{-1}\boldsymbol{\mu}_2+\ln P(\mathcal{C}_1)-\ln P(\mathcal{C}_2).  \end{aligned}

二次項 \mathbf{x}^T\Sigma^{-1}\mathbf{x} 被消去,後驗機率變成吃進 \mathbf{x} 的線性函數 a(\mathbf{x}) 的S型函數的回傳值:

\displaystyle  P(\mathcal{C}_1\vert\mathbf{x})=\frac{1}{1+\exp(-(\mathbf{w}^T\mathbf{x}+w_0))}

兩個類別的邊界即為超平面 \mathbf{w}^T\mathbf{x}+w_0=0 (見“超平面”),因為邊界滿足 P(\mathcal{C}_1\vert\mathbf{x})=P(\mathcal{C}_2\vert\mathbf{x})=\frac{1}{2}。我們得到了一個線性分類器,稱為線性判別分析:若 \mathbf{w}^T\mathbf{x}\ge -w_0,則數據點 \mathbf{x} 歸於 \mathcal{C}_1,否則歸於 \mathcal{C}_2

 
對於 k>2,相同的推導步驟可得

\displaystyle  a_j(\mathbf{x})=\mathbf{w}_j^T\mathbf{x}+w_{j0},~~1\le j\le k

其中

\displaystyle\begin{aligned}  \mathbf{w}_j^T&=\Sigma^{-1}\boldsymbol{\mu}_j\\  w_{j0}&=-\frac{1}{2}\boldsymbol{\mu}_j^T\Sigma^{-1}\boldsymbol{\mu}_j+\ln P(\mathcal{C}_j).  \end{aligned}

後驗機率即是吃進 \mathbf{x} 的線性函數 a_j(\mathbf{x}) 的 Softmax 函數的回傳值,

\displaystyle  P(\mathcal{C}_j\vert\mathbf{x})=\frac{\exp(\mathbf{w}_j^T\mathbf{x}+w_{j0})}{\sum_{i=1}^n\exp(\mathbf{w}_i^T\mathbf{x}+w_{i0})},~~1\le j\le k

指數是單調遞增函數,故多類別線性判別分析法則如下:若 l=\arg\max_{1\le j\le k}\mathbf{w}_j^T\mathbf{x}+w_{j0},則數據點 \mathbf{x} 歸於 \mathcal{C}_l。如果類別 \mathcal{C}_i\mathcal{C}_ji\neq j,在 \mathbb{R}^d 所屬的區域相鄰,則邊界為超平面 \mathbf{w}_i^T\mathbf{x}+w_{i0}=\mathbf{w}_j^T\mathbf{x}+w_{j0}

 
最大似然估計

最後一個問題:如何求得先驗機率 P(\mathcal{C}_j),平均數向量 \boldsymbol{\mu}_j,以及共變異數矩陣 \Sigma?下面介紹一個統計學的經典方法──最大似然估計法 (maximum likelihood estimation)。考慮兩個類別的情況,假設樣本記為 \mathcal{X}=\{\mathbf{x}_i,r_i\}_{i=1}^n,其中 \mathbf{x}_i\in\mathbb{R}^dr_i=1 表示 i\in \mathcal{C}_1r_i=0 表示 i\in\mathcal{C}_2。為簡化記號,設 P(\mathcal{C}_1)=\pi,則 P(\mathcal{C}_2)=1-\pi,並以 \mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_j,\Sigma) 表示常態分布的密度函數。如果 \mathbf{x}_i 來自 \mathcal{C}_1

\displaystyle  p(\mathbf{x}_i,\mathcal{C}_1)=p(\mathbf{x}_i\vert\mathcal{C}_1)P(\mathcal{C}_1)=\pi\mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_1,\Sigma)

如果 \mathbf{x}_i 來自 \mathcal{C}_2

\displaystyle  p(\mathbf{x}_i,\mathcal{C}_2)=p(\mathbf{x}_i\vert\mathcal{C}_2)P(\mathcal{C}_2)=(1-\pi)\mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_2,\Sigma)

寫出似然函數

\displaystyle\begin{aligned}  \mathcal{L}(\pi,\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma|\mathcal{X})&=p(\mathcal{X}|\pi,\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma)\\  &=\prod_{i=1}^n\left(p(\mathbf{x}_i,\mathcal{C}_1)\right)^{r_i}\left(p(\mathbf{x}_i,\mathcal{C}_2)\right)^{1-r_i}\\  &=\prod_{i=1}^n\left(\pi\mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_1,\Sigma)\right)^{r_i}\left((1-\pi)\mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_2,\Sigma)\right)^{1-r_i}.  \end{aligned}

對數是單調遞增函數,因此最大化 \mathcal{L} 等價於最大化 \ln \mathcal{L},如下:

\displaystyle  \ln \mathcal{L}=\sum_{i=1}^n\left(r_i[\ln \pi+\ln \mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_1,\Sigma)]+(1-r_i)[\ln (1-\pi)+\ln \mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_2,\Sigma)]\right).

分別對 \pi\boldsymbol{\mu}_j\Sigma 求導函數,設為零即得最佳解滿足的條件式。首先,

\displaystyle  \frac{\partial \ln \mathcal{L}}{\partial \pi}=\sum_{i=1}^n\left(\frac{r_i}{\pi}-\frac{1-r_i}{1-\pi}\right)

\frac{\partial\ln \mathcal{L}}{\partial\pi}=0。由此解出 \pi 的最大似然估計

\displaystyle  \hat{\pi}=\frac{n_1}{n}=\frac{n_1}{n_1+n_2}

其中 n_1=\sum_{i=1}^nr_i 表示 \mathcal{C}_1 的數據點數,n_2=n-n_1 表示 \mathcal{C}_2 的數據點數。再來,使用純量的向量導數公式 (見“矩陣導數”,SV-8,SV-10),

\displaystyle\begin{aligned}  \frac{\partial \ln \mathcal{L}}{\partial\boldsymbol{\mu}_1}&=\sum_{i=1}^nr_i\frac{\partial \ln\mathcal{N}(\mathbf{x}_i\vert\boldsymbol{\mu}_1,\Sigma)}{\partial\boldsymbol{\mu}_1}\\  &=\sum_{i\in\mathcal{C}_1}\frac{\partial}{\partial\boldsymbol{\mu}_1}\left(-\frac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_1)^T\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_1)+\alpha\right) \\  &=\sum_{i\in\mathcal{C}_1}\frac{\partial}{\partial\boldsymbol{\mu}_1}\left(-\frac{1}{2}\boldsymbol{\mu}_1^T\Sigma^{-1}\boldsymbol{\mu}_1+\boldsymbol{\mu}_1^T\Sigma^{-1}\mathbf{x}_i  -\frac{1}{2}\mathbf{x}_i^T\Sigma^{-1}\mathbf{x}_i+\alpha\right) \\  &=\sum_{i\in\mathcal{C}_1}\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_1)\\  &=\Sigma^{-1}\sum_{i\in\mathcal{C}_1}(\mathbf{x}_i-\boldsymbol{\mu}_1),  \end{aligned}

其中 \alpha 整合了指數函數乘數的對數 (包含 \vert\Sigma\vert)。令 \frac{\partial \ln \mathcal{L}}{\partial\boldsymbol{\mu}_1}=\mathbf{0},解得 \boldsymbol{\mu}_1 的最大似然估計即為樣本平均數向量

\displaystyle  \mathbf{m}_1=\hat{\boldsymbol{\mu}}_1=\frac{1}{N_1}\sum_{i\in\mathcal{C}_1}\mathbf{x}_i

按照同樣方式或使用對稱性,可得 \boldsymbol{\mu}_2 的最大似然估計

\displaystyle  \mathbf{m}_2=\hat{\boldsymbol{\mu}}_2=\frac{1}{N_2}\sum_{i\in\mathcal{C}_2}\mathbf{x}_i

最後推導共變異數矩陣 \Sigma 的最大似然估計。將似然函數裡涉及 \Sigma 的部分挑出來,使用跡數循環不變性,可得

\displaystyle\begin{aligned}  \ln \mathcal{L}&=-\frac{1}{2}\sum_{i=1}^nr_i\ln\vert\Sigma\vert-\frac{1}{2}\sum_{i=1}^nr_i(\mathbf{x}_i-\boldsymbol{\mu}_1)^T\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_1)\\  &~~~-\frac{1}{2}\sum_{i=1}^n(1-r_i)\ln\vert\Sigma\vert-\frac{1}{2}\sum_{i=1}^n(1-r_i)(\mathbf{x}_i-\boldsymbol{\mu}_2)^T\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_2)+\beta\\  &=-\frac{1}{2}\sum_{i\in\mathcal{C}_1}\ln\vert\Sigma\vert-\frac{1}{2}\sum_{i\in\mathcal{C}_1}\text{trace}\left((\mathbf{x}_i-\boldsymbol{\mu}_1)^T\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_1)\right)\\  &~~~-\frac{1}{2}\sum_{i\in\mathcal{C}_2}\ln\vert\Sigma\vert-\frac{1}{2}\sum_{i\in\mathcal{C}_2}\text{trace}\left((\mathbf{x}_i-\boldsymbol{\mu}_2)^T\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_2)\right)+\beta\\  &=-\frac{n}{2}\ln\vert\Sigma\vert-\frac{1}{2}\sum_{i\in\mathcal{C}_1}\text{trace}\left(\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_1)(\mathbf{x}_i-\boldsymbol{\mu}_1)^T\right)\\  &~~~-\frac{1}{2}\sum_{i\in\mathcal{C}_2}\text{trace}\left(\Sigma^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_2)(\mathbf{x}_i-\boldsymbol{\mu}_2)^T\right)+\beta\\  &=-\frac{n}{2}\ln\vert\Sigma\vert-\frac{1}{2}\text{trace}\left(\Sigma^{-1}\sum_{i\in\mathcal{C}_1}(\mathbf{x}_i-\boldsymbol{\mu}_1)(\mathbf{x}_i-\boldsymbol{\mu}_1)^T+\Sigma^{-1}\sum_{i\in\mathcal{C}_2}(\mathbf{x}_i-\boldsymbol{\mu}_2)(\mathbf{x}_i-\boldsymbol{\mu}_2)^T\right)+\beta\\  &=-\frac{n}{2}\ln\vert\Sigma\vert-\frac{n}{2}\text{trace}(\Sigma^{-1}S)+\beta,  \end{aligned}

其中 \beta 囊括所有與 \Sigma 無關的數,且

\displaystyle\begin{aligned}  S&=\frac{n_1}{n}S_1+\frac{n_1}{n}S_2\\  S_j&=\frac{1}{n_j}\sum_{i\in\mathcal{C}_j}(\mathbf{x}_i-\boldsymbol{\mu}_j)(\mathbf{x}_i-\boldsymbol{\mu}_j)^T,~~j=1,2.  \end{aligned}

使用跡數與行列式的導數公式 (見“跡數與行列式的導數”,det-2,tr-18),

\displaystyle  \frac{\partial \ln \mathcal{L}}{\partial \Sigma}=-\frac{n}{2}\Sigma^{-1}+\frac{n}{2}\Sigma^{-1}S\Sigma^{-1}

\frac{\partial \ln \mathcal{L}}{\partial \Sigma}=0,解出共變異數矩陣 \Sigma 的最大似然估計為 \hat{\Sigma}=S,也就是樣本 (有偏) 共變異數矩陣 S_1S_2 的加權平均 (將上式的 \boldsymbol{\mu}_j 替換為 \mathbf{m}_j)。

 
以上二類別的估計公式可立即推廣至多類別,所有參數的最大似然估計整理於下:對於 j=1,\ldots,kn_j 表示類別 \mathcal{C}_j 的數據點數,n=\sum_{j=1}^kn_k

\displaystyle\begin{aligned}  \pi_j&=\hat{P}(\mathcal{C}_j)=\frac{n_j}{n}\\  \mathbf{m}_j&=\hat{\boldsymbol{\mu}}_j=\frac{1}{n_j}\sum_{i\in\mathcal{C}_j}\mathbf{x}_i\\  S_j&=\frac{1}{n_j}\sum_{i\in\mathcal{C}_j}(\mathbf{x}_i-\mathbf{m}_j)(\mathbf{x}_i-\mathbf{m}_j)^T  \end{aligned}

以及

\displaystyle  S=\hat{\Sigma}=\frac{n_1}{n}S_1+\cdots+\frac{n_k}{n}S_k

 
線性判別分析假設樣本遵守常態分布,而且所有的類別有相同的共變異數矩陣,這使得類別邊界具有清晰的幾何意義,並導出簡潔的分類法則。然而,當實際數據不符合這些假設條件時,線性判別分析的分類效能往往不及其他競爭對手。另一個困難是,常態分布的最大似然估計不具備強健性 (robustness),也就是說,離群值 (outlier) 極可能造成巨大的參數估計偏差 (見“樣本平均數與樣本中位數,孰優孰劣?”)。

相關閱讀:
This entry was posted in 機器學習 and tagged , , , , , . Bookmark the permalink.

4 則回應給 線性判別分析

  1. 張盛東 說:

    周老師,我有一個問題。

    既然樣本平均數和樣本共變異數矩陣很容易受到離群值影響,為何不用樣本中位數作為樣本平均值的估計並且用其計算樣本共變異數矩陣?因為在常態分佈的假設下,常態分佈的mean和median相等。這樣是否可能提高似然估計的robustness?

  2. 張盛東 說:

    周老師,S型函數與 Softmax 函數那一節最後的P(Cj)應該是P(Cj | x)吧?

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s