2017年12月24日 星期日

[凸分析] 任意遞增函數為 quasiconcave 與 quasiconvex

=============
Theorem:
令 $f : \mathbb{R} \to \mathbb{R}$ 為遞增函數,則 $f$ 同時為 quasiconcave 與 quasiconvex。
=============

Proof:
考慮 $x,y \in \mathbb{R}$ 且 $\lambda \in (0,1)$,在不失一般性情況我們假設 $x>y$,則
\[
x > \lambda x + (1-\lambda)y > y
\]因為 $f$ 遞增,我們有
\[
f(x) > f( \lambda x + (1-\lambda)y) > f(y) \;\;\;\; (**)
\]由上述不等式,我們可寫
 $$
f(x) = \max\{f(x),f(y)\}
$$故由第一部分的不等式,我們有
\[
 \max\{f(x),f(y)\} > f( \lambda x + (1-\lambda)y)
\]此表明 $f$ 為 quasiconvex。

同理,由 $(**)$ 我們亦可寫下
\[
f(y) := \min\{f(x),f(y)\}
\]故由第二部分的不等式,我們有
\[
 f( \lambda x + (1-\lambda)y) >\min\{f(x),f(y)\}
\]此表明 $f$ 為 quasiconcave 至此證明完畢。$\square$

2017年12月19日 星期二

[凸分析] 凸性 與 齊次性 的關聯 (2):一些常見的結果

此文接續前篇對於凸函數 與 齊次函數的討論,主要是給出一些常見的結果。閱讀此文之前,建議讀者先回憶 凸性 與 齊次函數的定義,關於齊次函數以及一些相關例子,讀者可參閱 [凸分析] 凸性 與 齊次性 的關聯 (1):一些常見例子

以下我們首先給出 一組齊次函數的 連乘  或 連加 仍保持為 齊次函數 的條件:

================================
Theorem: Homogeneous Functions Algebra
1. 令 $f_1,...,f_m$ 為一組 定義在 convex cone $C \subset \mathbb{R}^n$ 上的 homogeneous functions,且對於 $i=1,...,m$ 而言, $f_i$ 具有 homogeneous of degree $\alpha_i$。則
\[
z(x) :=  \prod_{i=1}^m f_i(x) = f_1(x) \cdot f_2(x) \cdots f_m(x)
\]為 homogenous of degree $(\alpha_1 + ... +\alpha_m)$

2. 令 $f_1,...,f_m$ 為一組 定義在 convex cone $C \subset \mathbb{R}^n$ 上的 homogeneous functions,且對於 $i=1,...,m$ 而言, $f_i$ 具有相同 homogeneous of degree $\alpha$。則
\[z(x): = {\left( {\sum\limits_{i = 1}^m {{f_i}(x)} } \right)^\beta }\]為 homogenous of degree $(\alpha \beta)$
================================

Proof 1: 令 $x \in C$ 觀察
\[
z(tx) = \prod_{i=1}^m f_i(tx)
\]由於 $f_i$ 為 homogeneous of degree $\alpha_i$,我們有 $$f_i(tx) = t^{\alpha_i} f(x)$$ 且 $\forall t>0, i=1,2,...,m$,因此

\begin{align*}
  z(tx) &= \prod\limits_{i = 1}^m {{f_i}} (tx) \hfill \\
   &= \prod\limits_{i = 1}^m {{t^{{\alpha _i}}}{f_i}(x)}  \hfill \\
   &= {t^{{\alpha _1}}}{f_1}(x){t^{{\alpha _2}}}{f_2}(x)...{t^{{\alpha _m}}}{f_m}(x) \hfill \\
   &= {t^{\sum\limits_{i = 1}^m {{\alpha _i}} }}\left( {{f_1}(x){f_2}(x)...{f_m}(x)} \right) \hfill \\
   &= {t^{\sum\limits_{i = 1}^m {{\alpha _i}} }}z\left( x \right)
\end{align*}
上式對任意 $t>0$ 成立,此表明 $z(x)$ 為 homogeneous of degree $(\alpha_1 + ... +\alpha_m)$,至此證畢。$\square$


Proof 2: 令 $x \in C$ 觀察
\[z(tx): = {\left( {\sum\limits_{i = 1}^m {{f_i}(tx)} } \right)^\beta }\]由於 $f_i$ 為 homogeneous of degree $\alpha$,我們有 $$f_i(tx) = t^{\alpha} f(x)$$ 且 $\forall t>0, i=1,2,...,m$,因此

\begin{align*}
  z(tx)&: = {\left( {\sum\limits_{i = 1}^m {{f_i}(tx)} } \right)^\beta } \hfill \\
   &= {\left( {\sum\limits_{i = 1}^m {{t^\alpha }{f_i}(x)} } \right)^\beta } \hfill \\
   &= {\left( {{t^\alpha }\sum\limits_{i = 1}^m {{f_i}(x)} } \right)^\beta } \hfill \\
   &= \left( {{t^{\alpha \beta }}} \right){\left( {\sum\limits_{i = 1}^m {{f_i}(x)} } \right)^\beta } \hfill \\
  & = \left( {{t^{\alpha \beta }}} \right)z\left( x \right) \hfill \\
\end{align*}
上式對任意 $t>0$ 成立,此表明 $z(x)$ 為 homogeneous of degree $(\alpha \beta)$,至此證畢。$\square$



下面這個結果表明 一次齊次函數 如果具有 次可加性(subadditivity) 則 此函數必定為 convex。

================================
Theorem: Linear Homogeneity With Subadditivity Produces Convexity
令 $f: C \subset \mathbb{R}^n \to \mathbb{R}$ 為 linearly homogeneous function; i.e., $(f(tx) = t f(x), \; \forall t>0)$ 則 $f$ 為 convex 若且唯若 $f$ 滿足 subadditivity,亦即,對任意 $x,y \in C$,$$
f(x+y) \leq f(x) + f(y)
$$================================


Proof: $(\Rightarrow)$ 給定 $f$ 為 convex 且 linearly homogeneous 要證明對任意 $x,y \in C$,$$
f(x+y) \leq f(x) + f(y)
$$故給定 $x,y \in C$ 並且觀察
\begin{align*}
  f\left( {x + y} \right) &= f\left( {2\left( {\frac{{x + y}}{2}} \right)} \right) \hfill \\
   &= 2f\left( {\frac{{x + y}}{2}} \right)\;\;\;\; (*) \hfill \\
\end{align*} 上述最後一條等式成立因為我們使用了 $f$ 為 linearly homogeneous。接著由於 $f$ 為 convex 且 $1/2 \in (0,1)$ 故我們有
\[f\left( {\frac{{x + y}}{2}} \right) \leqslant \frac{1}{2}f\left( x \right) + \frac{1}{2}f\left( y \right)\]故
\[2f\left( {\frac{{x + y}}{2}} \right) \leqslant 2\left( {\frac{1}{2}f\left( x \right) + \frac{1}{2}f\left( y \right)} \right) = f\left( x \right) + f\left( y \right) \;\;\; (**)\]由 $(*)$ 與 $(**)$ 我們有
\[f\left( {x + y} \right) \leqslant f\left( x \right) + f\left( y \right)\]

$(\Leftarrow)$ 接著我們令 $x,y \in C$ 滿足 subadditivity: $f(x+y) \leq f(x) + f(y)$,我們要證明 $f$ 是 convex。現在給定 $\lambda \in (0,1)$ ,注意到由於 $x,y \in C$ 且 $C$ 為 convex cone,故對任意 $\lambda \in (0,1)$,我們有 $\lambda x \in C$ 與 $(1-\lambda y) \in C$。現在利用 已知的 subadditivity,我們 觀察
\[
f(\lambda x + (1- \lambda) y)  \leq f(\lambda x) + f( (1-\lambda)y) \;\;\; (\star)
\]利用 $f$ 為 linearly homogeneous,$f(\lambda x)=\lambda f( x)$ 且 $f( (1-\lambda)y) = (1-\lambda)f( y)$亦即,
\[f(\lambda x) + f((1 - \lambda )y) = \lambda f\left( x \right) + \left( {1 - \lambda } \right)f\left( y \right) \;\;\;\; (\star \star)\]由 $(\star)$ 與 $(\star \star)$ 式,我們可得
\[f(\lambda x + (1 - \lambda )y) \leqslant \lambda f\left( x \right) + \left( {1 - \lambda } \right)f\left( y \right)\]此表明 $f$ 為 convex in $C$。$\square$

[凸分析] 凸性 與 齊次性 的關聯 (1):定義 與 一些常見例子


Definition:  Homogenous Function of Degree Alpha
令 $C \subset \mathbb{R}^n$ 為 convex cone。我們說函數 $f: C \to \mathbb{R}$ 為  $\alpha$ 次齊次函數 (homogeneous of degree $\alpha \in \mathbb{R}$) 若下列條件成立:
對任意 $x \in C$,
\[
f(t x) = t^\alpha f(x),\;\;\; \forall t >0
\]

Comments:
1.上述 齊次函數 定義可以推廣到不是在 convex cone上而是任意向量空間,但一般做 convex cone的假設是為了 其他在凸分析上的 用途。在比較深入的凸分析教材中,可能會探討所謂 廣義凸性(generalized convexity)比如 quasi-convexity, quasi-concavity, semi-strictly quasi-convexity 等等,則此時函數定義域 需要是凸集。

2. 注意到 degree of homogeneity $\alpha \in \mathbb{R}$,意指 此 $\alpha$ 為任意實數,正數,負數,零 或者其他都可以。

3. 我們說 $f$ 為 homogenous of degree $0$ 若對任意 $x \in C$,
\[
f(tx) = f(x), \forall t>0
\]若 $f$ 為 homogenous of degree $1$ 一般稱之為 linearly homogeneous ,亦即 對任意 $x \in C$
\[
f(tx) = t f(x), \forall t>0
\] 以下我們看幾個例子:

Example 1
考慮 需求函數 (Demand Function)
\[
D(p,R) := \frac{R}{p}
\]其中 $p > 0$ 為 price of a good 且 $R>0$ 為 income,試證此函數為 homogeneous of degree 0

Proof: 令 $C :=\{(p,r): p>0, R>0\}$,則此集合為一個 convex cone,現在觀察對任意 $(p,r) \in C$,
\[
D(tp, tR) = \frac{tR}{tp} = \frac{R}{p} = t^0 D(p,r)
\]上式對任意 $t>0$ 成立,故由定義可知 $D(p,r)$ 為 homogenous of degree zero,至此證畢。$\square$


Example 2:
考慮 生產函數 (Production Function)
\[
f(L,K) := L^{1/3} K^{2/3}
\]其中 $L>0$ 為投入的勞力 (labour) 且 $K>0$ 為 投入資本 capital,試證生產函數 $f$ 為 homogenous of degree $1$

Proof: 令 $C:= \{(L,K): L>0, K>0\}$ ,則可知 $C$ 為一個 convex cone。現在取任意 $(L,K) \in C$,我們有
\[f(tL,tK): = {\left( {tL} \right)^{1/3}}{\left( {tK} \right)^{2/3}} = t{\left( L \right)^{1/3}}{\left( K \right)^{2/3}} = t f(L,K)
\]上式對任意 $t>0$ 成立,故由定義可知 $D(p,r)$ 為 homogenous of degree zero,至此證畢。$\square$

Comments:
上述的生產函數 在經濟學中被稱為 Cobb-Douglas Function,在一般 計量經濟 中通常記作
\[
f(K,L) := A K^\alpha L^{1- \alpha}
\]此式子命名來自 兩位美國學者 C. W. Cobb 與 P. H. Douglas 於 1927 提出,此函數用以估計生產量,但事實上此式早在 1900之前就由 瑞士經濟學家 Knut Wicksell 已率先提出。此為軼事與本文無關只是單純提及。


接著我們看個上述生產函數的推廣,在(個體)經濟學中常見的 Cobb-Douglas Function:

===============
Theorem: Generalized Cobb-Douglas Function is Homogenous of Degree Alpha
令 $A>0, x_i>0, \alpha_i>0, i =1,2,...,n$,定義 Cobb-Douglas 函數
\[
f(x): = Ax_1^{{\alpha _1}}x_2^{{\alpha _2}} \cdots x_n^{{\alpha _n}}
\]則 $f$ 為 Homogeneous of degree $\alpha := \sum_{i=1}^n \alpha_i$
===============

Proof: 令 $C:=\{x=(x_1,x_2,...,x_n): x_i >0\}$ ,則不難發現 $C$ 為  convex cone,現在取任意 $x = (x_1,...,x_n) \in C$,我們觀察
\begin{align*}
  f(tx) &: = A\left( {t{x_1}} \right)_{}^{{\alpha _1}}\left( {t{x_2}} \right)_{}^{{\alpha _2}} \cdots \left( {t{x_n}} \right)_{}^{{\alpha _n}} \hfill \\
   &= A{t^{\sum\limits_{i = 1}^n {{\alpha _i}} }}\left( {{x_1}} \right)_{}^{{\alpha _1}}\left( {{x_2}} \right)_{}^{{\alpha _2}} \cdots \left( {{x_n}} \right)_{}^{{\alpha _n}} \hfill \\
   &= {t^{\sum\limits_{i = 1}^n {{\alpha _i}} }}\underbrace {A\left( {{x_1}} \right)_{}^{{\alpha _1}}\left( {{x_2}} \right)_{}^{{\alpha _2}} \cdots \left( {{x_n}} \right)_{}^{{\alpha _n}}}_{ = f\left( x \right)} \\
&= {t^\alpha }f\left( x \right)
\end{align*} 其中 $\alpha:= \sum_{i=1}^n \alpha_i$,上述等式對任意 $t>0$ 成立,故 $f$ 為 Homogeneous of degree $\alpha := \sum_{i=1}^n \alpha_i$,至此證畢。$\square$


Comments:
上述 Cobb-Douglas function $f$ 亦俱備 log-linear 性質,亦即對 $f$ 取 $\log (.)$ 之後為線性函數:觀察
\begin{align*}
 \log \left( {f\left( x \right)} \right) &= \log \left( {Ax_1^{{\alpha _1}}x_2^{{\alpha _2}} \cdots x_n^{{\alpha _n}}} \right) \hfill \\
   &= \log \left( A \right) + {\alpha _1}\log \left( {x_1^{}} \right) + {\alpha _2}\log \left( {x_2^{}} \right) + ... + {\alpha _n}\log \left( {x_n^{}} \right) \hfill \\
\end{align*} 由上述結果不難看出 $\log(f)$ 為 linear functions of $\log(x_1), \log(x_2),...,\log(x_n)$




[訊號處理] 離散時間 的 Parserval's Theorem

Discrete-Time Parserval's Theorem:
令 $x[n]$ 為 取實數值的離散時間訊號; i.e., 對任意整數 $n$ 而言,
 $x[n] \in \mathbb{R}$ 且令 $X(\omega) := DTFT(x[n])$,則下列等式成立  \[\sum\limits_{n =  - \infty }^\infty  {{{\left| {x\left[ n \right]} \right|}^2} = \frac{1}{{2\pi }}\int_0^{2\pi } {{{\left| {X\left( \omega  \right)} \right|}^2}d\omega } } \]其中 $DTFT(x[n])$ 表示 對 $x[n] $ 取 離散時間傅立葉轉換 (Discrete Time Fourier Transform, DTFT),亦即
\[
X(\omega) := \sum_{n = -\infty}^\infty x[n] e^{-j \omega n}
\]

Comments:
1. 在訊號處理或者訊號分析的領域,上述定理有可被賦予的物理意義:一般我們把 Parserval's theorem 等式左方 稱為 訊號 $x[n]$ 的 總能量 (total energy),現在觀察右式表示對 $|X(\omega)|$ 做 $2 \pi$ 週期積分,因為 DTFT 為 $2\pi$ 週期,故右式 事實上是對所有的頻率積分等於 總能量,那麼被積分項 $|X(\omega)|$ 很自然的被稱作 頻譜密度 spectral (power) density
2. 考慮 $x(t)$ 為訊號,則 $\sum\limits_{n =  - \infty }^\infty  {{{\left| {x\left[ n \right]} \right|}^2}} $ 稱作此訊號的 total energy
3. 數學上的觀點,上述提及的 總能量  = 自己與自己在 函數空間做 內積 。
4. 另一種數學上的觀點:若 $x[n]$ 想成向量空間的一個向量,則總能量 $$\sum\limits_{n =  - \infty }^\infty  {{{\left| {x\left[ n \right]} \right|}^2}} $$ 可想成 $x[n]$ "長度" 的平方 ($l_2$-norm 平方)。
5. Fourier 轉換 是保持長度(norm preserving)的一種 線性轉換。
6. 上述 定理中我們雖僅考慮任意 實數值訊號 $x[n] \in \mathbb{R}, \;\;\; \forall \; n$ ,但事實上此定理對 $x[n] \in \mathbb{C}$ 亦成立。
7. 當然,有 離散時間傅立葉轉換,亦有 反轉換,稱作 Inverse Discrete Time Fourier Transform, IDTFT。


以下我們給出上述定理的證明

Proof: 給定 $x[n]$ 與其對應的 DTFT $X(\omega)$,我們要證明 Parserval's theorem 等式成立。故首先定義一個 輔助函數
\[
y\left[ m \right]: = \sum\limits_{n =  - \infty }^\infty  {x\left[ n \right]x\left[ {n - m} \right]}  \;\;\; (\star)
\] 上述 $y$ 一般稱作是 $x$ 的 自相關函數 (auto-correlation function) 。定義此函數的好處是我們可以利用 DTFT的各種性質來求證 Parserval's theorem。現在觀察 \[
 y\left[ 0 \right] = \sum\limits_{n = - \infty }^\infty {x\left[ 0 \right]x\left[ {n - 0} \right]} = \sum\limits_{n = - \infty }^\infty {{x^2}\left[ n \right]}
 \]此為要證明的 Parserval's 等式左方。

接著我們觀察 auto-correlation function $y[n]$ 為時域訊號,可透過 IDTFT 來表示 (why?),亦即
\[y[n] = \frac{1}{{2\pi }}\int_0^{2\pi } Y (\omega ){e^{j\omega n}}d\omega \]注意到當 $n=0$,我們有
\[y[0] = {\left. {\left( {\frac{1}{{2\pi }}\int_0^{2\pi } Y (\omega ){e^{j\omega n}}d\omega } \right)} \right|_{n = 0}} = \frac{1}{{2\pi }}\int_0^{2\pi } Y (\omega )d\omega \]至此不難發現若 $Y(\omega) = |X(\omega)|^2$ 則 Parseval's theorem得證。要達成此目標,我們觀察上述  $(\star)$式 事實上等價為
\[ y[n] = x[n] * x[-n]
\]其中 $*$ 表示 convolution。 (讀者應回憶 convolution 定義並自行驗證。) 回憶 DTFT 的 convolution 性質:亦即 時域 convolution 等價為 頻域做 multiplication。 故我們對 $(\star)$式 兩邊同取 DTFT 可得
\[y\left[ m \right] = x\left[ m \right]*x\left[ { - m} \right] \Rightarrow Y\left( \omega  \right) = X\left( \omega  \right){X^*}\left( \omega  \right) = |X(\omega)|^2\]
注意到在此我們使用了另一個 FACT: 若 $X(\omega) := DTFT(x[n])$,則 $DTFT(x[-n]) = X^*(\omega)$。

現在比較 $y[0]$ 可得\[\left\{ \begin{gathered} y\left[ 0 \right] = \sum\limits_{n = - \infty }^\infty {{x^2}\left[ n \right]} \hfill \\ y[0] = \frac{1}{{2\pi }}\int_0^{2\pi } {{{\left| {X\left( \omega \right)} \right|}^2}} d\omega \hfill \\ \end{gathered} \right. \Rightarrow \frac{1}{{2\pi }}\int_0^{2\pi } {{{\left| {X\left( \omega \right)} \right|}^2}} d\omega = \sum\limits_{n = - \infty }^\infty {{x^2}\left[ n \right]} \]至此證明完畢。$\square$

2017年12月6日 星期三

[機率論] 關於含有 Factorial 函數的求導的注意事項 - Erlang 分佈為例

在某些情況,我們可能會希望對含有 Factorial (比如說 $k!$, $k \in \mathbb{N}$) 的函數 取導數,但在求導 的過程中有些細微的部分需要多加留意。以下我們用一個例子來體現。

令 $m \in \mathbb{N}$,考慮隨機變數 $X$ 配備 Erlang 分佈
$$
F_X(x) := 1 - \sum_{k=0}^{m-1} \frac{(\lambda x)^k}{k!} e^{-\lambda x}, x>0
$$ 試證 其 機率密度函數 (Probability Density Function, pdf) $f_X$ 滿足
\[{f_X}\left( x \right) = \frac{{{{(\lambda x)}^{\left( {m - 1} \right)}}}}{{\left( {m - 1} \right)!}}\lambda {e^{ - \lambda x}}\]

(FALSE) Proof:
首先注意到分佈函數可導,故我們可利用 分佈函數的導數 為 密度函數 的性質 ($F'(x) = f(x)$),來求得 $f_X$。現在對 $F_X$ 求導
\[\frac{d}{{dx}}{F_X}(x) =  - \sum\limits_{k = 0}^{m - 1} {\left( {\underbrace {\frac{{k{{(\lambda x)}^{k - 1}}\lambda }}{{k!}}{e^{ - \lambda x}}}_{**} + \frac{{{{(\lambda x)}^k}}}{{k!}}\left( { - \lambda } \right){e^{ - \lambda x}}} \right)} \]注意到summation的第一項 $(**)$,讀者可能會很自然地認為 $**$ 可寫成
\[\frac{{k{{(\lambda x)}^{k - 1}}\lambda }}{{k!}}{e^{ - \lambda x}} = \frac{{k{{(\lambda x)}^{k - 1}}\lambda }}{{k\left( {k - 1} \right)!}}{e^{ - \lambda x}}\]
然後試圖對 分子分母的 $k$對消。但注意到此項 是在 summation內部,若對 分子 與 分母 進行 對消將產生問題,因為當 $k=0$ 時候會出現 難以處理未定義的 $-1!$ 。到此我們無法繼續進行,該怎麼避免這種問題呢?我們必須將可能出問題的 $k=0$ 項次分開討論。

Proof:
首先改寫

\[{F_X}(x): = 1 - \sum\limits_{k = 0}^{m - 1} {\frac{{{{(\lambda x)}^k}}}{{k!}}} {e^{ - \lambda x}} = 1 - (\lambda x){e^{ - \lambda x}} - \sum\limits_{k = 1}^{m - 1} {\frac{{{{(\lambda x)}^k}}}{{k!}}} {e^{ - \lambda x}}\]再取導數
\[\small
 \begin{align*}
  \frac{d}{{dx}}{F_X}(x) &= 0 - \left[ {\left( { - \lambda } \right){e^{ - \lambda x}} + \sum\limits_{k = 1}^{m - 1} {\left( {\frac{{k{{(\lambda x)}^{k - 1}}\lambda }}{{k!}}{e^{ - \lambda x}} + \frac{{{{(\lambda x)}^k}}}{{k!}}\left( { - \lambda } \right){e^{ - \lambda x}}} \right)} } \right] \hfill \\
   & =   - \left[ {\left( { - \lambda } \right){e^{ - \lambda x}} + \left( {\left( {\frac{{{{(\lambda x)}^0}}}{{1!}} - \frac{{{{(\lambda x)}^1}}}{{1!}}} \right) + \left( {\frac{{2{{(\lambda x)}^1}}}{{2!}} - \frac{{{{(\lambda x)}^2}}}{{2!}}} \right) + ... + \left( {\frac{{\left( {m - 1} \right){{(\lambda x)}^{\left( {m - 1} \right) - 1}}}}{{\left( {m - 1} \right)!}} - \frac{{{{(\lambda x)}^{\left( {m - 1} \right)}}}}{{\left( {m - 1} \right)!}}} \right)} \right)\lambda {e^{ - \lambda x}}} \right]\hfill \\
  & = \lambda {e^{ - \lambda x}} - \left( {\frac{{{{(\lambda x)}^0}}}{{1!}}\underbrace { - \frac{{{{(\lambda x)}^1}}}{{1!}} + \frac{{2{{(\lambda x)}^1}}}{{2!}}}_{ = 0}\underbrace { - \frac{{{{(\lambda x)}^2}}}{{2!}} + \frac{{3{{(\lambda x)}^2}}}{{3!}}}_{ = 0}\underbrace { - \frac{{{{(\lambda x)}^3}}}{{3!}} + }_{ = 0}...\underbrace { + \frac{{\left( {m - 1} \right){{(\lambda x)}^{\left( {m - 1} \right) - 1}}}}{{\left( {m - 1} \right)!}}}_{ = 0} - \frac{{{{(\lambda x)}^{\left( {m - 1} \right)}}}}{{\left( {m - 1} \right)!}}} \right)\lambda {e^{ - \lambda x}} \hfill \\
  & = \lambda {e^{ - \lambda x}} - \left( {1 - \frac{{{{(\lambda x)}^{\left( {m - 1} \right)}}}}{{\left( {m - 1} \right)!}}} \right)\lambda {e^{ - \lambda x}} \hfill \\
   &= \frac{{{{(\lambda x)}^{\left( {m - 1} \right)}}}}{{\left( {m - 1} \right)!}}\lambda {e^{ - \lambda x}} \hfill \\
\end{align*} \]上述第三行等式為 telescoping sum,中間各項等於 $0$。至此證畢。$\square$

2017年12月2日 星期六

[機率論] 一類條件機率取極限的問題

Theorem: 令 $X,Y$ 為 jointly continuous,
\[
\lim_{h \to 0} P((X,Y) \in A| x < X \leq x +h) = \int_{-\infty}^\infty 1_A(x,y) f_{Y|X}(y|x)dy
\]其中 $1_A(\cdot)$ 為 indicator function

Proof: 首先觀察
\begin{align*}
  P((X,Y) \in A|x < X \leqslant x + h) &= \frac{{P(\left\{ {(X,Y) \in A} \right\},\left\{ {x < X \leqslant x + h} \right\})}}{{P(x < X \leqslant x + h)}} \hfill \\
   &= \frac{{P(\left\{ {\left( {X,Y} \right) \in A} \right\} \cap \left\{ {X \in \left( {x,x + h} \right]} \right\})}}{{P(x < X \leqslant x + h)}} \hfill \\
   &= \frac{{P(\left\{ {\left( {X,Y} \right) \in A} \right\} \cap \left\{ {\left( {X,Y} \right) \in \left( {x,x + h} \right] \times \mathbb{R}} \right\})}}{{P(x < X \leqslant x + h)}} \hfill \\
   &= \frac{{P(\left( {X,Y} \right) \in A \cap \left\{ {\left( {x,x + h} \right] \times \mathbb{R}} \right\})}}{{P(x < X \leqslant x + h)}} \hfill \\
   &= \frac{{\iint\limits_{A \cap \left\{ {\left( {x,x + h} \right] \times \mathbb{R}} \right\}} {{f_{XY}}\left( {t,s} \right)dtds}}}{{\int_x^{x + h} {{f_X}\left( t \right)dt} }}
\end{align*} 由於 $1_{A \cap B} = 1_A \cdot 1_B$ 我們有
\begin{align*}
  P((X,Y) \in A|x < X \leqslant x + h) &= \frac{{\iint\limits_{} {{1_A}\left( {t,s} \right){1_{\left( {x,x + h} \right] \times \mathbb{R}}}\left( {t,s} \right){f_{XY}}\left( {t,s} \right)dtds}}}{{\int_x^{x + h} {{f_X}\left( t \right)dt} }} \hfill \\
   &= \frac{{\int_{ - \infty }^\infty  {\int_x^{x + h} {{1_A}\left( {t,s} \right){f_{XY}}\left( {t,s} \right)dtds} } }}{{\int_x^{x + h} {{f_X}\left( t \right)dt} }} \hfill \\
   &= \frac{{\int_x^{x + h} {\int_{ - \infty }^\infty  {{1_A}\left( {t,s} \right){f_{XY}}\left( {t,s} \right)dsdt} } }}{{\int_x^{x + h} {{f_X}\left( t \right)dt} }} \hfill \\
\end{align*} 上下同乘 $1/h$ 我們得到
\[P((X,Y) \in A|x < X \leqslant x + h) = \frac{{\frac{1}{h}\int_x^{x + h} {\int_{ - \infty }^\infty  {{1_A}\left( {t,s} \right){f_{XY}}\left( {t,s} \right)dsdt} } }}{{\frac{1}{h}\int_x^{x + h} {{f_X}\left( t \right)dt} }}\]讓 $h \to 0$ 我們有
\begin{align*}
  P((X,Y) \in A|x < X \leqslant x + h) &= \frac{{\frac{1}{h}\int_x^{x + h} {\int_{ - \infty }^\infty  {{1_A}\left( {t,s} \right){f_{XY}}\left( {t,s} \right)dsdt} } }}{{\frac{1}{h}\int_x^{x + h} {{f_X}\left( t \right)dt} }} \hfill \\
   &= \frac{{\int_{ - \infty }^\infty  {{1_A}\left( {x,s} \right){f_{XY}}\left( {x,s} \right)ds} }}{{{f_X}\left( x \right)}} \hfill \\
   &= \int_{ - \infty }^\infty  {{1_A}\left( {x,s} \right)\frac{{{f_{XY}}\left( {x,s} \right)}}{{{f_X}\left( x \right)}}ds}  \hfill \\
   &= \int_{ - \infty }^\infty  {{1_A}\left( {x,s} \right){f_{Y|X}}\left( {s|x} \right)ds} 
\end{align*} 上述等式即為所求,至此證畢。 $\square$

[機率論] 連續隨機變數條件機率的定義

若 $X$ 是連續隨機變數,則其累積機率分配(cumulative distribution function, cdf)
\[
F_X(x) := P(X \leq x) = \int_{-\infty}^x f_X(t) dt
\]為 (對 $x$ ) 連續函數,故單點機率測度 $P(X=x) =0$。但若我們考慮條件機率的情況事情會變得稍微有點棘手,因為假設我們引入第二個連續隨機變數 $Y$且假設 $X,Y$ 為 jointly continuous,現在我們想計算 $P(Y \in C| X=x)$,由條件機率定義可知
\[
P(Y \in C| X=x) = \frac{P(X=x,Y=c)}{P(X=x)}
\]但此時我們發現因為 $P(X=x) =0$,分母是$ 0$。對於這種情況我們該怎麼對 連續隨機機變數定義其條件機率?或者更簡單的說,該怎麼計算(或者定義)$P(Y \in C| X=x)$?


要計算 $P(Y \in C| X=x)$,我們首先考慮
\[
\lim_{h \to 0} P(Y \in C| x<X\leq x + h)
\]對任意 $h>0$而言,上述條件機率可寫成
\[P(Y \in C|x < X \leqslant x + h) = \frac{{P(Y \in C,x < X \leqslant x + h)}}{{P(x < X \leqslant x + h)}}\]注意到分子部分等價為
\[P(Y \in C,x < X \leqslant x + h) = P\left( {\left( {X,Y} \right) \in \left( {x,x + h} \right] \times C} \right)\]若 $X,Y$ 為 jointly continuous,則
\begin{align*}
  P(Y \in C|x < X \leqslant x + h) &= \frac{{P(Y \in C,x < X \leqslant x + h)}}{{P(x < X \leqslant x + h)}} \hfill \\
   &= \frac{{\int_x^{x + h} {\int_C^{} {{f_{XY}}\left( {t,s} \right)dsdt} } }}{{\int_x^{x + h} {{f_X}\left( t \right)dt} }} \hfill \\
\end{align*} 對上式分子分母同除 $1/h$ 並且讓 $h \to 0$ ,利用下文中的 FACT可得
\begin{align*}
  \mathop {\lim }\limits_{h \to 0} P(Y \in C|x < X \leqslant x + h) &= \mathop {\lim }\limits_{h \to 0} \frac{{\frac{1}{h}\int_x^{x + h} {\int_C^{} {{f_{XY}}\left( {t,s} \right)dsdt} } }}{{\frac{1}{h}\int_x^{x + h} {{f_X}\left( t \right)dt} }} \hfill \\
   &= \frac{{\int_C^{} {{f_{XY}}\left( {x,s} \right)ds} }}{{{f_X}\left( x \right)dt}} \hfill \\
\end{align*} 由上述 極限,我們可定義 在給定 $X$ 條件之下 ,$Y$的條件機率密度函數,記作 $f_{Y|X}$ 如下:
================

Definition: Conditional Probability and Conditional Density: 對任意 $x$ 滿足 $f_X(x) >0$,給定 $X$ 條件之下 ,$Y$的條件機率密度函數定義為
\[
f_{Y|X}(y|x) := \frac{f_{XY}(x,y)}{f_X(x)}
\]由 $f_{Y|X} $,我們可定義 給定條件 $X=x$ 之下,事件 $Y \in C$ 的條件機率為
\[
P(Y \in C|X=x):= \int_C f_{Y|X}(y|x)dy
\]================

Comments: 1. 同理,我們可定義 Conditional CDF 記作 $F_{Y|X}$ 滿足
\[
F_{Y|X}(y|x) := P(Y \leq y| X = x) = \int_{-\infty}^y f_{Y|X}(t|x) dt
\]2. 讀者應不難驗證 $\int_{-\infty}^{\infty} f_{Y|X}(y|x)dy = 1$。

================
FACT: 令 $X$ 為隨機變數配備 機率密度函數(probability density function, pdf) $f_X$,則
\[\lim_{h \to 0}\frac{1}{h}\int_x^{x + h} {{f_X}\left( t \right)dt}  = F_X'(x) =f_X(x)\]其中 $F_X$ 為 $X$ 累積分配函數(cdf)。
================

Proof: 首先觀察
\[\frac{1}{h}\int_x^{x + h} {{f_X}\left( t \right)dt}  = \frac{1}{h}\left( {{F_X}\left( {x + h} \right) - {F_X}\left( x \right)} \right)\]現在讓 $h\to 0$我們有
\[\mathop {\lim }\limits_{h \to 0} \frac{1}{h}\int_x^{x + h} {{f_X}\left( t \right)dt}  = \mathop {\lim }\limits_{h \to 0} \frac{1}{h}\left( {{F_X}\left( {x + h} \right) - {F_X}\left( x \right)} \right) = F'(x)\]若 density 存在,則利用 pdf 是 cdf的微分的事實,
\[
F_X'(x) = f_X(x)
\] 至此證明完畢。$\square$

2017年11月22日 星期三

[數理統計] 對平均值的信賴區間

考慮一組 i.i.d. random samples $X_1,X_2,...,$ 配備 期望平均 $m$ 與 變異數 (variance) $\sigma^2$,假設變異已知,但期望平均 $m$ 為未知。我們想對 $m$ 進行估計。一般的做法是採用 sample mean 估計量
\[
M_n := \frac{1}{n}\sum_{i = 1}^n X_i
\] 則我們知道 sample mean $M_n \to m$ 當 $n \to \infty$ in probability,此性質稱作 (weak) consistency estimator,但實際上,大多情況之下我們僅僅只能量測有限 $n$ (比如只能做有限次實驗),則我們想問該如何描述紹述 $M_n$ 有多接近 $m$ 呢?此想法為建構 信賴區間的動機:對某些 $\delta >0$ 而言,我們定義
\[
P(m \in [M_n -\delta, M_n + \delta] ) = 1-\alpha
\]其中 $[M_n -\delta, M_n + \delta] $ 稱作信賴區間(Confidence interval) 且 $1-\alpha$ 稱作 信心水準 (Confidence level),因此 Confidence interval 是一個隨機集合 而信心水準是 此隨機集合包含未知參數 $m$ 的機率。一般而言,在實務上多使用 $1-\alpha \in [0.9, 0.99]$。

Comments:
1. 上述 $M_n$ 不但為為 真實平均 (或者期望值) $m$  的 consistent estimator 且 亦為不偏 (unbiased)估計量。


問題:給定 $\alpha$,我們想問該如何選取參數 $\delta$ 使得 \[
P(m \in [M_n -\delta, M_n + \delta] ) = 1-\alpha
\]成立?

要回答此問題,我們首先觀察
\[m \in [{M_n} - \delta ,{M_n} + \delta ] \Leftrightarrow {M_n} - \delta  \leqslant m \leqslant {M_n} + \delta \]亦即 $- \delta  \leqslant m - {M_n} \leqslant \delta $ 此等價為
\[
|M_n - m| \leq \delta
\]故此
\[P(m \in [{M_n} - \delta ,{M_n} + \delta ]) = P\left( {\left| {{M_n} - m} \right| \leqslant \delta } \right)\]現在我們取
\[
\delta := \frac{\sigma y}{\sqrt{n}}
\] 上述 $\delta$ 的取法給了我們極大的方便,因為
\begin{align*}
  P\left( {\left| {{M_n} - m} \right| \leqslant \delta } \right) &= P\left( {\left| {{M_n} - m} \right| \leqslant \frac{{\sigma y}}{{\sqrt n }}} \right) \hfill \\
  & = P\left( {\left| {\frac{{{M_n} - m}}{{\sigma /\sqrt n }}} \right| \leqslant y} \right) \hfill \\
  & = P\left( { - y \leqslant \frac{{{M_n} - m}}{{\sigma /\sqrt n }} \leqslant y} \right) \hfill \\
   & = {F_{\frac{{{M_n} - m}}{{\sigma /\sqrt n }}}}\left( y \right) - {F_{\frac{{{M_n} - m}}{{\sigma /\sqrt n }}}}\left( { - y} \right) \hfill \\
\end{align*} 其中 $F_X(\cdot)$ 表示隨機變數 $X$的累積機率密度函數 (cdf)。由 中央極限定理(Central Limit Theorem, CLT) 我們可知 \[{F_{\frac{{{M_n} - m}}{{\sigma /\sqrt n }}}}\left( y \right) \to \Phi \left( y \right)\]其中 $\Phi(y)$ 為 標準常態分配的累積機率密度函數 (standard normal cdf) 滿足
\[
\Phi(y) := \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^y e^{-t^2/2} dt
\]故當 $n$ 足夠大的時候,
\[\begin{gathered}
  P\left( {\left| {{M_n} - m} \right| \leqslant \delta } \right) = {F_{\frac{{{M_n} - m}}{{\sigma /\sqrt n }}}}\left( y \right) - {F_{\frac{{{M_n} - m}}{{\sigma /\sqrt n }}}}\left( { - y} \right) \hfill \\
   \approx \Phi \left( y \right) - \Phi \left( { - y} \right) \hfill \\
\end{gathered} \]由於 $\Phi(y)$ 為 even density function,由附註的 FACT ,我們有 $\Phi(-y) =1 - \Phi(y)$ 故可推得
\[P\left( {\left| {{M_n} - m} \right| \leqslant \delta } \right) \approx \Phi \left( y \right) - \Phi \left( { - y} \right) = 2\Phi \left( y \right) - 1\]現在觀察上述結果,我們得到以下結論:若我們希望
\[P\left( {\left| {{M_n} - m} \right| \leqslant \delta } \right) = 1 - \alpha \]則等價求解 $2\Phi \left( y \right) - 1 = 1- \alpha$ 亦即
\[\Phi \left( y \right) = 1 - \frac{\alpha }{2}\]注意到上式與 $n$ 無關 且 也與 $X_i$ 的 pdf 無關!此式的解我們將其記作 $y_{\alpha/2}$ 在 MATLAB 我們可求解
\[
y_{\alpha/2} = \text{norminv(1-alpha/2)}
\]以下我們將上述討論記作以下結果。

====================
Theorem:
固定信心水準 $1-\alpha$ ,則其平均值 對應的 $(1-\alpha) \%$ 信賴區間 為
\[m \in \left[ {{M_n} - \frac{1}{{\sqrt n }}\sigma {y_{\alpha /2}},{M_n}  + \frac{1}{{\sqrt n }}\sigma {y_{\alpha /2}}} \right]\]
====================


接著我們看個例子:
====================
Example: 令 $X_1,X_2...$為 i.i.d. random samples 滿足 $\sigma = 2$。若我們已知 $M_{100} = 1 $ 試求出其真實平均 落在信心水準分別為 95% 與 99% 信賴區間:
====================

Proof: 對於 $95\%$ 信心水準,其 $\alpha = 1 - 0.95 = 0.05$。接著我們求解 $y_{\alpha/2}$,利用 MATLAB :
$$
y_{\alpha/2} = \text{norminv(1 - 0.05/2)} = 1.96
$$故其對應的 95% 信賴區間為
\begin{align*}
  m &\in \left[ {{M_n} - \frac{1}{{\sqrt n }}\sigma {y_{\alpha /2}},{M_n} + \frac{1}{{\sqrt n }}\sigma {y_{\alpha /2}}} \right] \\
&= \left[ {1 - \frac{1}{{\sqrt {100} }}2\left( {1.96} \right),1 + \frac{1}{{\sqrt {100} }}2\left( {1.96} \right)} \right] \hfill \\
  & = \left[ {0.608,1.392} \right] \hfill \\
\end{align*}

對於 99% 信心區間可用同樣方法,不難求得
\[
m \in [0.4848,  1.5152]
\]細節留給讀者自行練習。

由上述例子可看出當 信心水準越大,對應的信賴區間越大。另外注意到因為
\[m \in \left[ {{M_n} - \frac{1}{{\sqrt n }}\sigma {y_{\alpha /2}},{M_n}  + \frac{1}{{\sqrt n }}\sigma {y_{\alpha /2}}} \right]\]不難看出如果想要縮減信賴區間,另一種方法則是 增加 量測值,也就是把 $n$ 提高。


Comments:
注意到上述討論我們假設 random samples 已知變異數。若變異數未知,則我們使用變異數的不偏估計量 $S_n^2$ 取代,亦即使用
\[
S_n^2 = \frac{1}{N-1} \sum_{i=1}^n (X_i -m)^2
\]則前述的信賴區間變成
\[m \in \left[ {{M_n} - \frac{1}{{\sqrt n }}{S_n}{y_{\alpha /2}},{M_n} + \frac{1}{{\sqrt n }}{S_n}{y_{\alpha /2}}} \right]
\]證明我們暫且略過。



附註
=================
FACT:
令 $X$ 為配備 even density function $f$ 的隨機變數。若 $F$ 為任意 even density function $f$ 的 cdf函數,則
$$
F(-x) = 1 - F(x)
$$=================
Proof:
令 $f$ even density function ,則由 cdf 定義可知
\[F\left( { - x} \right): = P(X \leq -x) = \int_{ - \infty }^{ - x} {f\left( t \right)dt} \]現在引入變數變換,令 $t := -y$則
\[\int_{ - \infty }^{ - x} {f\left( t \right)dt}  =  - \int_\infty ^x {f\left( { - y} \right)dy} \]由於 $f$ 為 even function 故 $f(-y) = f(y) \;\; \forall y $,所以
\begin{align*}
  \int_{ - \infty }^{ - x} {f\left( t \right)dt}  &=  - \int_\infty ^x {f\left( { - y} \right)dy}  \hfill \\
  & =  - \int_\infty ^x {f\left( y \right)dy}  \hfill \\
  & = \int_x^\infty  {f\left( y \right)dy}  = 1 - F\left( x \right) \hfill \\
\end{align*} 換言之,我們得到
\[
F(-x) = 1 - F(x)
\]至此得證。$\square$

上述 FACT 亦可用反證法求證,在此不贅述。

2017年11月6日 星期一

[機率論] 一些常用的機率上界估計-以 exponential 隨機變數為例

令 $X$ 為 exponential 隨機變數配備參數 $\lambda =1$,亦即 $X \sim exp(1)$。試求 $P(X\geq a)$ 並分別利用 Markov inequality, Chebyshev inequalty 與 Chernoff bound 估計此機率的上界。


Proof: 首先回憶 $X \sim exp(\lambda)$ 具有 機率密度
\[f_X\left( x \right): = \left\{ \begin{gathered}
  \lambda {e^{ - \lambda x}}\begin{array}{*{20}{c}}
  {}&{x \geqslant 0}
\end{array} \hfill \\
  0\begin{array}{*{20}{c}}
  {}&{x < 0}
\end{array} \hfill \\
\end{gathered}  \right.\]故 $\lambda =1$ 我們有
\[f_X\left( x \right): = \left\{ \begin{gathered}
  {e^{ - x}}\begin{array}{*{20}{c}}
  {}&{x \geqslant 0}
\end{array} \hfill \\
  0\begin{array}{*{20}{c}}
  {}&{x < 0}
\end{array} \hfill \\
\end{gathered}  \right.\] 現在我們計算
\[
P\left( {X \geqslant a} \right) = \int_a^\infty  {{e^{ - x}}dx}  =  - \left( {0 - {e^{ - a}}} \right) = {e^{ - a}}
\]

接著我們用 Markov inequality 估計:令 $a>0$則
\begin{align*}
  P\left( {X \geqslant a} \right) &\leqslant \frac{{E\left[ X \right]}}{a} \hfill \\
   &= \frac{1}{a}\int_0^\infty  {x{e^{ - x}}dx}  \hfill \\
   &= \frac{1}{a}
\end{align*}
再來是利用 Chebyshev inequalty  估計:
\begin{align*}
  P\left( {X \geqslant a} \right)& \leqslant \frac{{E\left[ {{X^2}} \right]}}{{{a^2}}} \hfill \\
   &= \frac{1}{{{a^2}}}\int_0^\infty  {{x^2}{e^{ - x}}dx}  \hfill \\
   &= \frac{2}{a^2}
\end{align*}

最後用 Chernoff bounds
\[P\left( {X \geqslant a} \right) \leqslant \mathop {\min }\limits_{s \geqslant 0} {e^{ - sa}}{M_X}\left( s \right)\]其中$M_X(s)$ 為 moment generating function 滿足
\[{M_X}\left( s \right): = E\left[ {{e^{sX}}} \right] = \int_0^\infty  {{e^{sx}}{e^{ - x}}dx}  = \frac{1}{{1 - s}}\]故欲求 Chrnoff bound 我們需要求解下列最佳化問題
\[\begin{gathered}
  \min {e^{ - sa}}{M_X}\left( s \right) \hfill \\
  {\text{subject to }}s \geqslant 0 \hfill \\
\end{gathered} \]此問題等價於
\[\begin{gathered}
  \min \frac{{{e^{ - sa}}}}{{1 - s}} \hfill \\
  {\text{subject to }}s \geqslant 0 \hfill \\
\end{gathered} \]現在我們考慮當 $s \in (0,1)$ 時,我們可求對 $\frac{{{e^{ - sa}}}}{{1 - s}}$ 求一階導數並令其為零,來得到臨界點,亦即
\begin{align*}
 & \frac{d}{{ds}} \left( \frac{1}{{1 - s}}{e^{ - sa}} \right)= 0 \hfill \\
  &\Rightarrow \frac{{{e^{ - as}}}}{{{{(1 - s)}^2}}} - \frac{{a{e^{ - as}}}}{{1 - s}} = 0 \hfill \\
  & \Rightarrow s = 1 - \frac{1}{a} \hfill \\
\end{align*} 注意到上式中我們要求 $a > 1$否則 $s<0$ 違反最佳化問題的 拘束條件。 再者我們驗證二階導數
\[{\left. {\frac{{{d^2}}}{{d{s^2}}}\left( {\frac{1}{{1 - s}}{e^{ - sa}}} \right)} \right|_{s = 1 - 1/a}} = {a^3}{e^{ - \left( {1 - \frac{1}{a}} \right)a}} > 0\]
故可知 $s=1-1/a$  為 local minimum。我們將其記作
\[
s^* := 1 - \frac{1}{a}
\]現在將其帶回 Chernoff bound 並計算
\begin{align*}
  P\left( {X \geqslant a} \right) &\leqslant \mathop {\min }\limits_{s \geqslant 0} {e^{ - sa}}{M_X}\left( s \right) \hfill \\
   &= {e^{ - {s^*}a}}{M_X}\left( {{s^*}} \right) \hfill \\
   &= {e^{1 - a}}a
\end{align*} 下圖顯示了當 $a \in [1,6]$ 時候,不同上界估計的圖形


上圖顯示在 $a in [1,6]$ 之間,各類機率上界估計互有輸贏。

現在我們考慮 $a \in [6,20]$ 並繪製 log-plot 可得下圖

上圖可發現當 $a$ 夠大,則 Chernoff bound 展現壓倒性的優勢。一般而言,當 $a$ 夠大的時候$P(X \geq a)$ 的上界估計採用 Chernoff bound 大多優於 Markov inequality 或者 Chebyshev inequality 是非常不錯的選擇。


附註:關於 Chernoff bound 可以非常容易地從 Markov inequality 推得:令 $s \geq 0$ 現在觀察以下事件
\[
\{X \geq a\} = \{s X \geq  sa\} = \{e^{sX} \geq e^{sa}\}
\]故若對上述事件取機率可得
\[
P(X \geq a) = P(e^{sX} \geq e^{sa})
\]由於$s \geq 0$我們有 $e^{sX} \geq 0$ 且 $e^{sa} >0$ 對任意 $a$,故利用 Markov inequality可推得
\[
 P(e^{sX} \geq e^{sa}) \leq \frac{E[e^{sX}]}{e^{sa}}
\]亦即
\[
P(X \geq a)  \leq \frac{E[e^{sX}]}{e^{sa}} , \;\; \forall s \geq 0
\]由於上述不等式對任意 $s\geq 0$成立,且不等式左方與 $s$ 無關,故可知
\[
P(X \geq a)  \leq \min_{s \geq 0 }\frac{E[e^{sX}]}{e^{sa}}
\]也就是 Chernoff bound。

[機率論] 關於配備 Pareto Density 隨機變數 的一個簡單例子

令 $X$ 為隨機變數配備標準 Pareto density,亦即其機率密度函數 $f_X$ 滿足
\[f_X\left( x \right): = \left\{ \begin{gathered}
  \frac{2}{{{x^3}}}\begin{array}{*{20}{c}}
  {}&{x \geqslant 1}
\end{array} \hfill \\
  0\begin{array}{*{20}{c}}
  {}&{o.w.}
\end{array} \hfill \\
\end{gathered}  \right.\]
(a) 對任意 $a \geq 1$ 試求 $P(X \geq a)$
(b) 延續 (a),利用 Markov inequality 求其上界。

Proof (a)
\begin{align*}
  P(X \geqslant a) &= 1 - P\left( {X < a} \right) \hfill \\
   &= 1 - \int_1^a {\frac{2}{{{x^3}}}dx}  \hfill \\
   &= 1 - 2\left( {\left. {\frac{{{x^{ - 2}}}}{{ - 2}}} \right|_1^a} \right) \hfill \\
   &= \frac{1}{{{a^2}}} \hfill \\
\end{align*}

Comments: 讀者可注意到上述結果亦可 透過直接計算
\[P\left( {X \geqslant a} \right) = \int_a^\infty  {\frac{2}{{{x^3}}}dx}  = 2\left( {\frac{1}{{ - 2}}\left. {{x^{ - 2}}} \right|_a^\infty } \right) = {a^{ - 2}}\]


Proof (b):
注意到 $X$ 為取值非負隨機變數,故對任意 $a \geq 1$,利用 Markov inequality, 我們有
\begin{align*}
  P(X \geqslant a) &\leqslant \frac{{E\left[ X \right]}}{a} \hfill \\
   &= \frac{1}{a}\int_{ - \infty }^\infty  {x{f_X}\left( x \right)dx}  \hfill \\
   &= \frac{1}{a}\int_1^\infty  {x\frac{2}{{{x^3}}}dx}  \hfill \\
   &= \frac{2}{a}\int_1^\infty  {\frac{1}{{{x^2}}}dx}  \hfill \\
   &=  - \frac{2}{a}\left( {\left. {{x^{ - 1}}} \right|_1^\infty } \right) =  - \frac{2}{a}\left( {0 - 1} \right) = \frac{2}{a} \hfill \\
\end{align*}

Comments:
1. 注意到上述 (b) 部分透過 Markov inequality 所得到的上界只有在 $a \geq 2$ 才有效力,因為當 $a \in [1,2]$ 之間時,我們得到 $2/a >1$。但由於機率測度不能超過 $1$,上述 $2/a$ 上界在 $a \in [1,2]$ 之間對我們的 $P(X \geq a)$ 的估計並無任何幫助。


2. 注意到 $a \geq 1$ 故不難得證
\[\frac{1}{{{a^2}}} \leqslant \frac{2}{a} \] 此說明了在此分佈之下,利用 Markov inequality 所得到的機率上界過鬆。下圖顯示了 $a \in [1,5]$ 的機率 與 Markov inequality所得的上界



2017年10月30日 星期一

[數學分析] 一類 分式與極小值 的不等式

Theorem
對任意 $i=1,...,m$,若 $a_i \geq 0$ 且 $b_i \geq 0$ 則下列不等式成立
\[
\frac{\sum_{i=1}^m a_i }{\sum_{i=1}^m b_i } \geq \min_i \frac{a_i}{b_i}
\]

Proof:
令 $i^*$ 為 某 index $i$ 使得 $\min_i \frac{a_i}{b_i}$ 成立,亦即 $i^*$ 滿足
\[\frac{{{a_{{i^*}}}}}{{{b_{{i^*}}}}} = \mathop {\min }\limits_i \frac{{{a_i}}}{{{b_i}}}\]我們要證明定理中的不等式成立。以下以各個擊破的方法來求證:

CASE 1:首先注意到若 $a_{i^*} = 0$ 則我們欲證明的不等式自動成立。

CASE 2: 故 假設 $a_{i^*} >0$,注意到若 $b_{i^*} =0$ 則我們得到兩邊不等式為無窮,故不等式仍然成立,故我們不妨假設  $a_{i^*} >0$ 且 $b_{i^*} > 0$ (*),現在觀察
\[
\frac{{\sum\limits_{i = 1}^m {{a_i}} }}{{\sum\limits_{i = 1}^m {{b_i}} }} = \frac{{{a_{{i^*}}} + \sum\limits_{i = 1}^{m - 1} {{a_i}} }}{{{b_{{i^*}}} + \sum\limits_{i = 1}^{m - 1} {{b_i}} }} = \frac{{{a_{{i^*}}}\left( {1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{a_i}}}{{{a_{{i^*}}}}}} } \right)}}{{{b_{{i^*}}}\left( {1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{b_i}}}{{{b_{{i^*}}}}}} } \right)}}\;\;\;\; (**)
\]注意到對任意 $i$ 而言,我們有
\[\frac{{{a_{{i^*}}}}}{{{b_{{i^*}}}}} \leqslant \frac{{{a_i}}}{{{b_i}}}\]又因為 $(*)$ 我們可推得對任意 $i$ 而言,下式成立
\[\frac{{{b_i}}}{{{b_{{i^*}}}}} \leqslant \frac{{{a_i}}}{{{a_{{i^*}}}}}\]故
\[\frac{{{b_i}}}{{{b_{{i^*}}}}} \leqslant \frac{{{a_i}}}{{{a_{{i^*}}}}} \Rightarrow 1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{b_i}}}{{{b_{{i^*}}}}}}  \leqslant 1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{a_i}}}{{{a_{{i^*}}}}}} \]現在將此不等式代入 $(**)$ 我們得到
\[\frac{{{a_{{i^*}}}\left( {1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{a_i}}}{{{a_{{i^*}}}}}} } \right)}}{{{b_{{i^*}}}\left( {1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{b_i}}}{{{b_{{i^*}}}}}} } \right)}} \leqslant \frac{{{a_{{i^*}}}\left( {1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{a_i}}}{{{a_{{i^*}}}}}} } \right)}}{{{b_{{i^*}}}\left( {1 + \sum\limits_{i = 1}^{m - 1} {\frac{{{a_i}}}{{{a_{{i^*}}}}}} } \right)}} = \frac{{{a_{{i^*}}}}}{{{b_{{i^*}}}}}\]亦即
\[\frac{{\sum\limits_{i = 1}^m {{a_i}} }}{{\sum\limits_{i = 1}^m {{b_i}} }} \leqslant \frac{{{a_{{i^*}}}}}{{{b_{{i^*}}}}} = \mathop {\min }\limits_i \frac{{{a_i}}}{{{b_i}}}\]至此得證。$\square$

2017年9月27日 星期三

[訊號與系統] Frequency Modulation 淺嚐:Chirp Signal

回憶一般 弦波訊號
\[
z(t) = A \cos(2 \pi f_0 t + \phi) = Re\{A e^{j (2 \pi f_0 t + \phi)}\}
\]其中 $A$ 為振幅 ,$f_0$ 為頻率,$\phi$ 為相位。由上式,我們可以定義 $z(t)$ 的 "角度" 記作
\[
\psi(t) := 2 \pi f_0 t + \phi
\]由於上式為 linear in $t$ ,我們可觀察
\[
\frac{d}{dt} \psi(t) = 2 \pi f_0 := \omega_i(t)
\] 故 $z(t)$ 的角度 隨時間的瞬時變化率 為 $2 \pi f_0$  其單位為 (rad/s) ,若我們將其除掉 $2\pi$ 即可得到 (瞬時)頻率 $f_0$ (單位為 Hz)。

上述想法可以被進一步推廣如下:

Frequency Modulation (FM) Signal
現在我們將上述 $z(t)$ 做進一步簡單的推廣:假設
\[
x(t) := A \cos(\psi (t)) = Re\{e^{j \psi(t)}\}
\]則我們可仿造前述的方法來定義 瞬時頻率 (instantaneous frequency),亦即我們先對 $\psi(t)$ 對 $t$ 微分,可得瞬時角頻率 (instantaneous angular frequency)
\[
\omega_i(t) :=  \frac{d}{dt} \psi(t)
\]若對上式兩邊同除以 $2 \pi$ ,可得瞬時頻率 (instantaneous frequency), 記作 $f_i(t)$, 如下
\[
f_i(t) :=\frac{1}{2\pi}\omega_i(t) =  \frac{1}{2 \pi} \frac{d}{dt} \psi(t)
\]單位為 Hz。

以下我們看個 FM 調頻中的一類特殊例子,假設我們想要創造一組 弦波訊號 使其 頻率可以包含一段我們感興趣頻段,比如說我們想創造一組聲音其頻率 從 300 Hz 並且一路往上到 800 Hz,一種常見的做法是採用 chirp signal 來達成,其特性如下:給定初始頻率 $f_{int}$ 與 終點頻率 $f_{end}$, chirp signal 保證訊號頻率在 $f_{int}$ 到 $f_{end}$ 之間以連續且線性方式改變(比如遞增或者遞減),故此法又稱線性調頻。

FM Signal Example: Chirp, or Linear Swept Frequency or Linear FM
以下我們考慮一類特殊的 FM 訊號,亦即我們取
\[
\psi(t) := 2\pi \mu t^2 + 2\pi f_0 t+\phi
\]則對應的瞬時頻率為
\[
f_i(t) := \frac{1}{2 \pi} \frac{d}{dt} \psi(t)  = 2 \mu t + f_0
\]亦即我們發現瞬時頻率隨時間遞增,且在 $t=0$時後我們有 起始頻率 $f_0$。讀者可使用 下列 matlab code 來聽聽看 chirp 訊號 (起始的 瞬時頻率為 300Hz 一路往上到 800Hz):

% Generate an play a chirp signal

fsamp = 11025; % sampling frequency

dur = 2;
mu = 125;
Amp = 6;
f0 = 300;

dt = 1/fsamp;
tt = 0 : dt : dur;

psi = 2*pi*(100 + f0*tt + mu*tt.*tt);
xx = real( Amp*exp(j*psi) );

soundsc( xx, fsamp );



[訊號與系統] Amplitude Modulation 淺嚐:Beat Signal

回憶一般 弦波訊號
\[
x(t) := A \cos(2 \pi f_c t + \phi)
\] 其中 $A$ 為振幅,$f_c$ 為 (載波) 頻率,$\phi$ 為相位。現在我們考慮將上述弦波做進一步簡單的推廣如下:假設 振幅 $A$ 不再是常數,而是隨時間變化的函數 比如 $A:=a(t)$ ,則我們得到\[
x(t) = a(t) \cos(2\pi f_c t +\phi)
\] 上述稱為 振幅調變 (Amplitude Modulation) 的一般形式,其中 $a(t)$ 為 依時間變化的函數 且一般而言,假設 $a(t)$ 的最高頻率 $f_a << f_c$。

Comments:
1. 震幅 $a(t)$ 隨時間變化,可看成 振幅 被 "調變"(modulation)
2. 一般實際應用上, $a(t)$ 為多半為實際帶有信息的訊號 (比如 聲音,歌聲,影像...) 且其最高頻率會遠遠低於 載子頻率 $f_c$ ,使得我們在做 AM 處理之後,$a(t)$ 訊號可被方便傳送。
3. 當然,對於 $z(t)$ 的推廣不僅僅限於頻率,我們也可以對其頻率推廣,比如說將固定 $f_c$ 改成 $f_c:=\psi(t)$ 使其成為與時間有關的函數,此法會得到所謂的 頻率調變(Frequency Modulation, FM) 我們會另外再開一篇文章描述之,在此不做贅述。

以下我們看個經典的AM例子:

AM Example: Beat Signal or Sinusoidal AM
以下我們看個特例:假設 $a(t) := A \cos(2 \pi f_a t)$ 則我們得到 AM 訊號如下
\[
x(t) = A \cos(2 \pi f_a t)  \cos(2\pi f_c t +\phi)
\]上述訊號可以透過 Inverse Euler formula 將其改寫為
\begin{align*}
  x(t) &= A\cos (2\pi {f_a}t)\cos (2\pi {f_c}t) \hfill \\
   &= A \left( \frac{{{e^{j2\pi {f_a}t}} + {e^{ - j2\pi {f_a}t}}}}{2} \right) \left( \frac{{{e^{j\left( {2\pi {f_c}t} \right)}} + {e^{ - j\left( {2\pi {f_c}t} \right)}}}}{2} \right) \hfill \\
   &= \frac{A}{4}\left( {{e^{j2\pi {f_a}t}} + {e^{ - j2\pi {f_a}t}}} \right)\left( {{e^{j\left( {2\pi {f_c}t} \right)}} + {e^{ - j\left( {2\pi {f_c}t} \right)}}} \right) \hfill \\
   &= \frac{A}{4}\left( {{e^{j2\pi \left( {{f_a} + {f_c}} \right)t}} + {e^{ - j2\pi \left( {{f_a} + {f_c}} \right)t}} + {e^{j2\pi \left( {{f_a} - {f_c}} \right)t}} + {e^{ - j2\pi \left( {{f_a} - {f_c}} \right)t}}} \right) \hfill \\
   &= \frac{A}{2}\cos \left( {2\pi \left( {{f_a} + {f_c}} \right)t} \right) + \frac{A}{2}\cos \left( {2\pi \left( {{f_a} - {f_c}} \right)t} \right) \hfill \\
\end{align*}
Comments:
1. 上述 $x(t)$ 可表為 兩弦波相加,一般又稱之為 beat signal,生活上的實際應用為比如說同時按下兩兩相鄰的鋼琴琴鍵。
2. 我們有兩種觀點看上述的 Beat Signal,首先是 $x(t)$ 可以視為是 震幅隨時間變化的弦波,故若使用 MATLAB 的 soundsc(.) 函數播放,則聲音聽起來會是漸強在接漸弱,第二種觀點則是上述 $x(t)$ 為兩個具有不同頻率 ($f_c+f_a$ 與 $f_c-f_a$)的 弦波相加,那麼聽起來便會是兩種弦波分別以不同頻率產生的聲音疊加而成。
3. 那麼該如何分辨何時只聽得到一組漸強漸弱得弦波 或者 聽到 不同頻率的弦波? 以下有一個一般性的判斷法則:令 $T$ 為 $x(t)$ 的最終持續時間,且定義 "頻寬" $B:= 2 f_a$ ,若
\[
T\cdot B <<1
\]則一般而言我們沒有辦法到底是一組弦波或者兩個不同頻率得弦波。此議題等價物理中的 Heisenberg's Uncertainty Principle。



2017年8月12日 星期六

[凸分析] 定義在凸集上的凸函數其相對極小即為全域極小

這次要介紹凸分析 或者 凸優化 中可以說是最重要的結果:

=====================
Theorem:
令 $f$ 為 convex on convex set $\Omega$。則
1. $f$ 的任意 相對極小點 $x^*$ (local minimum of $f$ on $\Omega$) 必為 全域極小點 (global minimum of $f$ on $\Omega$)
2. 上述 凸函數的極小點所成之集合為凸集,亦即
\[
S:= \{x \in \Omega: f(x) = f(x^*)\}
\]為凸集。
=====================

Proof (1):
利用反證法:令  $y \in  \Omega$ 為 $f$在 $\Omega$ 上的相對極小點,亦即存在 $\delta>0$ 使得鄰域 $$N_\delta(y):=\{z: ||z-y|| < \delta\} \subset \Omega$$ 且
\[
f(y) \leq f(x), \;\;\; \forall x \in N_\delta(y)
\] 但 $y$ 不為 global minimum :亦即 存在 $x^* \in \Omega$ 使得 \[
f(x^*) < f(y)
\]我們要證明此假設矛盾。

首先注意到 $x^* \notin N_\delta(y)$ 因為若不然,則 $f(y) \leq f(x^*)$ 此與假設不符。

現在,我們利用 $x^*$ 與 $y$ 來定義一個 新的點
\[
z(\alpha):= (1-\alpha)y +   \alpha x^*, \;\;\; \alpha := \frac{\delta}{2||y-x^*||} \in (0,1)
\]注意到 $z(\alpha) \in \Omega$ (因為 $\Omega$ 為凸集) 且我們觀察
\begin{align*}
  |z\left( \alpha  \right) - y|| &= \left\| {  (1 - \alpha ) y + \alpha x^* - y} \right\| \hfill \\
   &= \left\| {\alpha \left( {y - {x^*}} \right)} \right\| \hfill \\
   &\leqslant \frac{\delta }{{2\left\| {y - {x^*}} \right\|}}\left\| {y - {x^*}} \right\| < \delta  \hfill \\
\end{align*} 此結果表明
\[
z(\alpha) \in N_\delta(y)
\]
現在利用 $f$ 為凸函數性質,我們可知
\begin{align*}
  f(z(\alpha )) &= f\left( {\left( {1 - \alpha } \right) y + \alpha  x^*} \right) \hfill \\
   &\leqslant \alpha f\left( y \right) + \left( {1 - \alpha } \right)\underbrace {f\left( {{x^*}} \right)}_{f\left( y \right)} < f\left( y \right) \hfill \\
\end{align*} 上述不等式表明我們找到一個新的點 $z(\alpha) \neq y$  且 $z(\alpha) \in N_\delta(y)$ 使得
\[
f(z(\alpha)) < f(y)
\]此違反了 $y$ 在 $\Omega$ 上為相對極小的假設,得到矛盾。 故 $y$ 必為 global minimum。

Proof:(2)
接著我們證明上述 全域極小點 所成的集合為凸集,亦即
\[
S:= \{x \in \Omega: f(x) = f(x^*)\}
\]
首先觀察凸函數 $f$ 的 $c$-level set
\[
S_c:= \{x \in \Omega: f(x) \leq c\}
\]由下方的 FACT 可知 凸函數的 level set 為凸集。但因為 $x^*$ 為全域極小,故若取 $c:=f(x^*)$ 則
\[
S_c|_{c=f(x^*)}  =  \{x \in \Omega: f(x) \leq f(x^*)\} = \{x \in \Omega: f(x) = f(x^*)\} =S
\]由於為等式左方 $S_c$ 為凸集,故 $S$ 亦為凸集。 $\square$


=====================
FACT: 凸函數的 Sublevel Set 為凸集
令 $f: \Omega \subset \mathbb{R}^n \to \mathbb{R}$ 為凸函數,則對任意 $c \in \mathbb{R}$ ,其 $c$-level set
\[
S_c:=\{x: f(x) \leq c\}
\]為凸集
=====================

Proof: 利用凸函數定義,取 $x,y \in S_c$ 且 $\alpha \in [0,1]$ ,觀察
\[
f(\alpha x + (1-\alpha) y) \leq \alpha f(x) + (1-\alpha) f(y) \leq c 
\]故 $\alpha x + (1-\alpha) y\in S_c$,此表明 $S_c$ 為凸集。$\square$

2017年8月11日 星期五

[凸分析] 一階可導凸函數利用單點近似必定低估


Theorem: 
令 $f \in C^1$ 且 $f$ 為 convex on convex set $\Omega \subset \mathbb{R}^n$ 若且唯若 對任意 $x,y \in \Omega$ 而言,
\[
f(y) \geq f(x) + \nabla f(x) \cdot (y-x)
\]其中 $\nabla f(x) \cdot (y-x) := \nabla f(x)^T (y-x)$

給出證明之前我們先給一些直觀上的看法:

Comments:
1. 上述定理算是相當直覺,簡而言之就是說 affine (in $y$) function:
$ f(x) + \nabla f(x) (y-x)$ 可以做為 凸函數 $f$ 在 $x$ 點附近的 1 階 Taylor 近似,如下圖所示:



2. 注意到上述定理闡述的不等式對於所有 $x,y \in \Omega$ 都成立,也就是說透過 對$x$ 一階 Taylor 近似必定低估,一般 $f(x) + \nabla f(x) (y-x)$ 又稱作 global underestimaotr  of $f$。
3. 上述結果指出利用局部資訊 (一階導數) 可以得到 全域資訊 (global understametor )。
4. 若 $\nabla f(x) = 0$ 則對任意 $y \in \Omega$,我們有
\[
f(y) \geq f(x)
\]亦即 $x$ 為 全域及小點 (global minimizer) of $f$


以下我們給出證明

Proof: 先證明 $(\Rightarrow)$
令 $f \in C^1$ 且 $f$ 為 convex on convex set $\Omega \subset \mathbb{R}^n$,給定任意 $x,y \in \Omega$ ,我們要證
\[
f(y) \geq f(x) + \nabla f(x) (y-x)
\]
由於  $f$ 為 convex,令 $\alpha \in (0,1)$ 且定義
$$
z(\alpha) := \alpha x + (1-\alpha) y
$$則 $z(\alpha) \in \Omega$ 且由 $f$的凸性,我們有
\begin{align*}
  f\left( {z(\alpha )} \right) &= f\left( {\alpha x + \left( {1 - \alpha } \right)y} \right) \hfill \\
   &\leqslant \alpha f\left( x \right) + \left( {1 - \alpha } \right)f\left( y \right) \hfill \\
\end{align*} 由於 $\alpha \neq 0$ 我們可整理上式得到
\[\frac{{f\left( {\alpha x + \left( {1 - \alpha } \right)y} \right) - f\left( y \right)}}{\alpha } \leqslant f\left( x \right) - f\left( y \right)\]或者
\[\frac{{f\left( {y - \alpha \left( {y - x} \right)} \right) - f\left( y \right)}}{\alpha } \leqslant f\left( x \right) - f\left( y \right)\]取 $\alpha \to 0$,由於 $f\in C^1$ 我們不難看出上述不等式左方 為沿著 $y-x$ 的方向導數,故我們有
\[
\nabla f(y) \cdot (y-x) \leq f(x) -f(y)
\]或者
\[
f(x) \geq f(y) + \nabla f(y) \cdot (y-x)
\]上述結果對 任意 $x,y \in \Omega$ 成立,故我們將 $x,y$ 角色對換即得到定理要求的陳述。

接著我們證明$(\Leftarrow)$:
假設  對任意 $x,y \in \Omega$ 而言,
\[
f(y) \geq f(x) + \nabla f(x) (y-x) \;\;\;\;\; (**)
\]我們要證明 $f$ 為 convex。故令 $x_1, x_2 \in \Omega$ 與 $\alpha \in [0,1]$ ,並且我們額外定義
\[
\bar{x} := \alpha x_1 + (1- \alpha) x_2
'\]
則由假設可知 $x_1, x_2, \bar{x}$ 必定滿足 $(**)$,我們可寫下
\[\begin{gathered}
  f({x_1}) \geqslant f(\bar x) + \nabla f(\bar x)({x_1} - \bar x) \hfill \\
  f({x_2}) \geqslant f(\bar x) + \nabla f(\bar x)({x_2} - \bar x) \hfill \\
\end{gathered} \]現在對上述 第一條不等式 兩邊同乘上 $\alpha$ ,對 第二條不等式 兩邊乘上 $1- \alpha$ ,亦即
\begin{align*}
 & \alpha f({x_1}) \geqslant \alpha f(\bar x) + \alpha \nabla f(\bar x)({x_1} - \bar x) \hfill \\
  &\left( {1 - \alpha } \right)f({x_2}) \geqslant \left( {1 - \alpha } \right)f(\bar x) + \left( {1 - \alpha } \right)\nabla f(\bar x)({x_2} - \bar x) \hfill \\
\end{align*}
現在觀察
\begin{align*}
  \alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) &\geqslant \alpha f(\bar x) + \alpha \nabla f(\bar x)({x_1} - \bar x) \hfill \\
   & \hspace{10mm}+ \left( {1 - \alpha } \right)f(\bar x) + \left( {1 - \alpha } \right)\nabla f(\bar x)({x_2} - \bar x)
\end{align*}
將上式稍微做一下整理可得
\begin{align*}
  &\alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) \geqslant f(\bar x) + \nabla f(\bar x)\left( {\alpha ({x_1} - \bar x) + \left( {1 - \alpha } \right)({x_2} - \bar x)} \right) \hfill \\
   &\Rightarrow \alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) \geqslant f(\bar x) + \nabla f(\bar x)\underbrace {\left( {\alpha {x_1} + \left( {1 - \alpha } \right){x_2} - \bar x} \right)}_{ = 0} \hfill \\
   &\Rightarrow \alpha f({x_1}) + \left( {1 - \alpha } \right)f({x_2}) \geqslant f(\bar x) \hfill \\
\end{align*} 上述不等式表明 $f$ 為凸函數。$\square$


Comments:
1. 若 $f$ 為 $C^1$ strict convex 函數 on $\Omega$,則對任意 $x,y \in \Omega$ 而言,
\[
f(y) >f(x) + \nabla f(x) (y-x)
\]
2. 若 $f$ 為 concave 則利用 $-f$ 為 convex 特性可知 對於 concave 函數而言,定理的不等式變成: 對任意 $x,y \in \Omega$ 而言,
\[
f(y) \leq f(x) + \nabla f(x) (y-x)
\]

2017年8月3日 星期四

[最佳化理論] 有限維空間 無拘束最佳化問題 的二階充分必要條件

回顧先前我們討論過的一階必要條件,以下我們介紹有限維空間 無拘束最佳化問題的 二階必要與充分條件,首先是二階必要條件:

====================
Theorem: Second-Order Necessary Condition, SONC:
令 $S \subset \mathbb{R}^n$ 且令 $f \in C^2(S)$。若 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 則 對 在點 ${\bf x}^*$的任意可行方向 ${\bf d} \in \mathbb{R}^n$,我們有
1. $\nabla f({\bf x}^*) \cdot {\bf d} \geq 0$
2. 若 $\nabla f({\bf x}^*) \cdot {\bf d} =0$ 則 ${\bf d}^T \nabla^2 f({\bf x}^*) {\bf d} \geq 0$
====================

Proof: 由於 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 對條件 1 自動成立。我們僅需證明條件2。令  ${\bf d} \in \mathbb{R}^n$ 為 在點 ${\bf x}^*$的任意可行方向,故存在 $\bar \alpha >0$ 使得
\[
{\bf x}(\alpha) := {\bf x}^* + \alpha {\bf d} \in S, \;\;\; \alpha \in [0, \bar \alpha]
\] 利用 Taylor Theorem 對 ${\bf x}^*$ 展開到二階項,我們可得
\begin{align*}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right)
&= f\left( {{{\mathbf{x}}^*}} \right) + \nabla f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) \\
& \hspace{10mm}+ \frac{1}{2}{\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right)^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
 &  = f\left( {{{\mathbf{x}}^*}} \right) + \alpha \nabla f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + \frac{1}{2}{\alpha ^2}{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
\end{align*} 若 $\nabla f({\bf x}^*) \cdot {\bf d} =0$,則上式變成
\[\begin{gathered}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right) = f\left( {{{\mathbf{x}}^*}} \right) + \alpha \underbrace {\nabla f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}}_{ = 0} + \frac{1}{2}{\alpha ^2}{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
   = f\left( {{{\mathbf{x}}^*}} \right) + \frac{1}{2}{\alpha ^2}{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}} + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
\end{gathered} \]注意到若 ${{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}<0$ 對在 ${\bf x}^*$任意足夠小的鄰域成立,則我們得到
\[f\left( {{\mathbf{x}}\left( \alpha  \right)} \right) \leqslant f\left( {{{\mathbf{x}}^*}} \right) \]此與 ${\bf x}^*$是 local minimum 矛盾,故
\[{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}\geq 0\;\;\;\;\; \square\]

Comments:
讀者大概不難發現不論是一階 或者 二階 條件都是依賴 Taylor 定理展開式,依此來進行估計。


====================
Theorem: Second-Order Sufficient Condition, SOSC: 
令 $S \subset \mathbb{R}^n$ 且令 $f \in C^2(S)$。若 ${\bf x}^* \in S$ 為interior point 且$f$ 在該點有定義。若
1. $\nabla f({\bf x}^*) = {\bf 0}$
2. 對任意 ${\bf d} \in \mathbb{R}^n$,${\bf d}^T \nabla^2 f({\bf x}^*) {\bf d} >0$ (亦即 $\nabla^2 f({\bf x}^*)$ 為 positive definite),則 ${\bf x}^*$ 為 strict local minimum of $f$。
====================

Proof: 要證明 ${\bf x}^*$ 為 strict local minimum of $f$,我們僅需證明
\[
f({\bf x}) > f({\bf x^*}),\;\;\; {\bf x} \in \mathcal{N}({\bf x}^*)
\]其中 $\mathcal{N}({\bf x}^*)$為以 ${\bf x}^*$為中心所成之鄰域。故令  ${\bf d} \in \mathbb{R}^n$ 為 在點 ${\bf x}^*$的任意可行方向,故存在 $\bar \alpha >0$ 使得
\[
{\bf x}(\alpha) := {\bf x}^* + \alpha {\bf d} \in \mathcal{N}({\bf x}^*), \;\;\; \alpha \in [0, \bar \alpha]
\] 利用 $\nabla f({\bf x}^*) = {\bf 0}$ 與 Taylor Theorem 對 ${\bf x}^*$ 展開到二階項,我們可得
\begin{align*}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right)
&= f\left( {{{\mathbf{x}}^*}} \right) + \nabla f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) \\
& \hspace{10mm}+ \frac{1}{2}{\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right)^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right)\left( {{\mathbf{x}}\left( \alpha  \right) - {{\mathbf{x}}^*}} \right) + o\left( {{{\left\| \alpha  \right\|}^2}} \right) \hfill \\
 &  =f\left( {{{\mathbf{x}}^*}} \right) + \frac{1}{2}{\alpha ^2}\underbrace {{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}}_{ > 0} + o\left( {{{\left\| \alpha  \right\|}^2}} \right)\\
\end{align*}
在 ${\bf x}^*$ 附近鄰域而言,我們可推知
\[\begin{gathered}
  f\left( {{\mathbf{x}}\left( \alpha  \right)} \right) - f\left( {{{\mathbf{x}}^*}} \right) = \frac{1}{2}{\alpha ^2}\underbrace {{{\mathbf{d}}^T}\nabla^2 f\left( {{{\mathbf{x}}^*}} \right){\mathbf{d}}}_{ > 0} > 0 \hfill \\
\end{gathered} \]則 ${\bf x}^*$ 為 strict local minimum of $f$ 。$\square$

2017年8月2日 星期三

[最佳化理論] 有限維空間 無拘束最佳化問題 的一階必要條件

令 $f: \mathbb{R}^n \to \mathbb{R}$ 且 $S \subset \mathbb{R}^n$ 為 feasible set,現在我們考慮以下的 有限維度 (無拘束)最佳化問題
\[\begin{gathered}
  \min f\left( {\mathbf{x}} \right) \hfill \\
  s.t.{\mathbf{x}} \in S \subseteq {R^n} \hfill \\
\end{gathered} \]
我們想知道上述最佳化問題是否有解? 若有解則是哪一種解 e.g., 局部最佳解(local optimum)或者 全域最佳解(global optimum)?),以及上述的解是否能夠透過某種方法來將其描述。

對於上述有限維度最佳解的存在性問題一般可由 Weierstrass Extremum Theorem處理,在此不做贅述。以下我們討論 最佳解 存在的必要條件:更近一步地說是 最佳化理論中的求取 "局部最佳解" 的一階必要條件,在給出結果之前我們首先定義 可行方向 (Feasible Direction)如下:

======================
Definition: Feasible Direction
給定 ${\bf x} \in S \subset \mathbb{R}^n$,我們說向量 ${\bf d} \in \mathbb{R}^n$ 為 在 ${\bf x}$ 處的可行方向 (feasible direction) 若下列條件成立:存在常數 $ \bar \theta >0$ 使得對任意 $\theta \in [0, \bar\theta]$而言,
\[
{\bf x} + \theta {\bf d} \in S
\]======================

有了以上的可行方向的想法,其局部最佳解的一階必要條件有如下陳述:

======================
Theorem: (First Order Necessary Condition, FONC): 令 $S \subset \mathbb{R}^n$ 且 $f \in C^1$ on $S$ ( $f$為一階可導且導數連續)。若 ${\bf x}^*$ 為 局部極小解 (local minimum point) of $f$ over $S$ 則 對任意在 ${\bf x}^*$點上的可行方向 ${\bf d} \in \mathbb{R}^n$,我們有
\[
\nabla f({\bf x}^*) \cdot {\bf d} \geq 0
\]======================

Proof: 用反證法:假設存在 對 ${\bf x}^*$點上的可行方向 ${\bf d} \in \mathbb{R}^n$ 使得
\[
\nabla f({\bf x}^*) \cdot {\bf d} <0
\]我們要證明矛盾。由 ${\bf d}$為可行方向之定義可知: 存在 $\bar \theta>0$ 使得 對任意 $\theta \in [0, \bar\theta]$而言,
\[
{\bf x}(\theta) = {\bf x}^* + \theta {\bf d} \in S
\]現在利用 Taylor 定理對 ${\bf x}^*$ 展開 且利用已知假設 $\nabla f({\bf x}^*) \cdot {\bf d} <0$ 可得
\begin{align*}
  f({\bf x}(\theta )) &= f({{\bf x}^*}) + \nabla f({{\bf x}^*})({\bf x}(\theta ) - {{\bf x}^*}) + o(||{\theta}||) \hfill \\
   &= f({{\bf x}^*}) + \nabla f({ {\bf x}^*})\theta {\bf d} + o(||{\theta}||) \hfill \\
  &= f({{\mathbf{x}}^*}) + \theta \underbrace {\nabla f({{\mathbf{x}}^*}) \cdot {\mathbf{d}}}_{ < 0} + o(||\theta ||)\\
   &< f({{\bf x}^*}) \hfill \\
\end{align*} 上述最後一條不等式當 $\theta$ 足夠小的時候成立。此與我們假設 $ f({ {\bf x}^*})$最小 矛盾。$\square$


Comments:
1. 上述一階必要條件說明了若我們已經處在局部最佳解的位置 ${\bf x}^*$ 則 沿著任何其他可行方向移動都會增加 目標函數值,亦即 \[
\nabla f({\bf x}^*) \cdot {\bf d} \geq  0
\]
2. 關於上述使用的 Taylor Theorem 與 little-oh 符號可以參考下方補充說明或者BLOG其餘相關文章。
3. 上述 一階必要條件 使用到了 Taylor Theorem 的一階項,一般而言可視為對 ${\bf x}^*$ 的 一階近似估計。


若為無拘束情況,則我們有以下的衍生結果:

=======================
Corollary: 令 $S \subset \mathbb{R}^n$ 且 $f \in C^1$ on $S$。若 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 且若 ${\bf x}^*$ 為 $S$ 的 interior point 則
\[
\nabla f({\bf x}^*) = {\bf 0}
\]=======================

Proof Sketch: 設 ${\bf x}^*$ 為 local minimum point of $f$ over $S$ 則由 FONC可知 對任意 ${\bf d} \in \mathbb{R}^n$ 為在 ${\bf x}^*$點上的可行方向,我們有
\[
\nabla f({\bf x}^*) \cdot {\bf d} \geq 0
\]但因為 若 ${\bf x}^*$ 為 $S$ 的 interior point 故在該點${\bf x}^*$ 之可行方向 ${\bf d}$ 可在足夠小的區域選為任意方向來移動, 若 ${\bf d} \neq {\bf 0}$ 則 為了使 $\nabla f({\bf x}^*) \cdot {\bf d} = 0$ 成立,我們必定要求
\[
\nabla f({\bf x}^*)  = {\bf 0} \;\;\;\; \square
\]

Comments:
1. 上述 FONC 無拘束情況的 FONC 將原本最佳問題轉成求解有 $n$ 未知數的 $n$ 個系統方程問題。

2. 上述結果只保證 局部最佳解 的必要條件,如果想要得到充分條件,通常需要 $f$ 的二階導數的資訊也就是 需要 Hessian matrix,一般稱作 Second-Order Sufficient Condition我們在以後的文章會在提及。

3. 局部最佳解落在集合 $S$ 上,若此集合為 convex 集合且 $f$ 為凸函數,則局部最佳解 為 全域最佳解,相關結果可翻閱本 BLOG關於最佳化與凸分析的文章。

4. 上述最佳化問題可以與 變分不等式問題 (Variational Inequality Problem)等價,以下給出此結果:

=======================
Corollary 2: Relationship Between Optimization and Variational Inequality
令 $S \subset \mathbb{R}^n$ 為 convex 且 $f : \mathbb{R}^n \to \mathbb{R}$ 且 $f \in C^1(S)$。若 ${\bf x}^*$ 為 local minimum point of $f$ over $S$,則 ${\bf x}^*$ 為下列變分不等式問題的解:求 ${\bf x} \in S$ 使得 對任意 ${\bf x}' \in S$ 而言,
\[
\langle {\bf x}' - {\bf x}, \nabla f({\bf x}) \rangle \geq 0
\]=======================
Proof:
取 ${\bf d} := {\bf x}' - {\bf x}$ 即可。



補充定理
==================
Taylor Theorem with Little-oh Notation:
假設 $X \subset \mathbb{R}^n$ 為 open,且 ${\bf x} \in X$ 與 $f: X \to \mathbb{R}$ 為 $C^1$。則
\[
f({\bf x}+{\bf h}) = f({\bf x}) + Df({\bf x}) {\bf h} + o(||{\bf h}||),\;\;\; ||h|| \to 0
\]==================










2017年7月7日 星期五

[控制理論] 非線性模型預測控制 (0) - 引論

此文我們針對 非線性模型預測控制 (Nonlinear Model Predictive Control, NMPC) 做簡單的介紹: 是一套針對受控廠為 非線性系統 所發展的 回授控制理論,其控制律透過不斷求解最佳化問題而得。以下我們將其簡稱 NMPC。


基本想法:
假設你是一個圍棋高手正在與旗鼓相當的對手對弈,那麼你可能在心中盤算好幾步可能的走法 並同時 試圖來 "預測" 對手的路數,但輪到自己下子的時候,我們只能選取剛才在心中盤算出的所有走法中 " 最佳" 的那一步,並且只動那一步旗,動了之後都必須重新盤算上述過程。這個精神大概與 模型預測控制 本質幾乎相同。就是逐步最佳化,慢慢朝向目標前進。以下我們會逐步將此概念規範化,讀者可以讀完之後再回頭瞧瞧這個基本想法,也許會發現異曲同工之妙。


對 $n=0,1,2,...$,令 $x(n)$ 為當前系統狀態 且 $x^{ref}(n)$ 為 給定的參考目標軌跡。

模型預測控制 的 主要目的:
與一般控制問題相仿,模型預測控制的主要目的不外乎以下兩大類問題:
1. 鎮定(stabilizing)問題。使預測的狀態軌跡 盡可能趨近 零點。
2. 追蹤(tracking)問題。使 預測的狀態軌跡 盡可能跟隨給定的 參考目標軌跡 。

Comments:
1. 若稍有控制背景的讀者不難發現鎮定問題其實是追蹤問題的特例,更進一步地說,所謂追蹤控制問題是指:決定一組控制輸入 $u(n)$ 使得 $x(n)$ 能盡可能緊跟給定的參考狀態軌跡 $x^{ref}(n)$。換而言之,若當前狀態 $x(n)$ 與 $x^{ref}(n)$ 所差甚遠,則我們將盡可能控制此系統軌跡趨近 $x^{ref}(n)$:若當前狀態 $x(n) = x^{ref}(n)$ 則我們將盡可能控制此系統"維 持"在該狀態。

2. 關於鎮定問題提及的對零點追蹤 或者 追蹤問題的參考目標軌跡,都必須滿足隱藏的前提: 零點 或者 參考目標軌跡 必須是系統的 平衡點 (equilibrium point)。否則該系統無法進行鎮定 or 追蹤。


以下我們將介紹一類簡化的 NMPC 問題,並藉此展示此控制方法的一些基本精神。

一類簡化的非線性模型預測控制問題:
令 $x(n) \in X := \mathbb{R}^d$ 且 $u(n) \in U := \mathbb{R}^m$ 考慮 參考軌跡狀態(reference state trajectory),記作 $x^{ref}(\cdot)$,寫作
\[
x^{ref}(n) :=x_* :=0,\;\;\;\;\; n\geq 0
\]上述 對零點狀態追蹤問題 退化為 鎮定問題。

以控制理論的基本想法,我們很自然想要對 當前狀態 $x(n)$ 與 目標 $x^{ref} = 0$ 之誤差 做回授控制,為了達成此目標,我們期望控制力 $u$ 必須具備以下形式:
$$
u(n) := \mu (x(n))
$$其中 $\mu $ 為 將狀態 $x \in X$ 映到控制輸出值之集合 $U$ 之映射。

系統模型動態(System Model Dynamics):
模型預測控制的想法在於利用 系統的模型來進行未來狀態軌跡之預測,並且給出適當的目標函數後對其做最佳化求解控制力。更精確的說,系統動態 一般我們寫成下列方程
\[
x^+ = f(x,u) \;\;\;\;\; (*)
\]其中 $f: X \times U \to X$ 為已知且非線性 滿足 將狀態 $x$ 與控制輸出值 $u$ 映到下一個時刻的狀態 $x^+$ 。

預測模型動態(Predictive Model Dynamics):
現在給定當前狀態 $x(n)$,對任意給定控制序列
\[
u(0),u(1),...,u(N-1)
\]滿足控制區間 $N \geq 2$,我們可以用此 控制序列 以及 給定的當前狀態 作為初始狀態來迭代 $(*)$,如此可得 預測狀態軌跡 $x_u$ 滿足下式:
\[
x_u(0) = x(n); \;\;\;\; x_u(k+1) = f(x_u(k), u(k)),\;\;\;\; k=0,...,N-1 \;\;\;\; (**)
\]由於 $x(n)$ 被作為初始狀態來進行預測,故我們 $x_u(k)$ 實際上為狀態 $x(n+k)$ (對時刻 $t_{n+k}$ ) 的 預測狀態。因此我們得到在時刻 $t_n, t_{n+1},...,t_{n+N}$ 對系統狀態的預測 (與狀態輸入序列 $u(0),...,u(N-1)$有關)。

每時刻的最佳控制力求解:
現在我們使用最佳控制的想法來決定 $u(0),...,u(N-1)$ 使得預測狀態 $x_u$ 能盡可能靠近 $x^{ref} = 0$。為了達成此目標,我們可選定 一距離函數
\[
l(x_u(k),u(k)) : = ||x_u(k)||^2 + \lambda ||u(k)||^2
\]其中 $||\cdot||$ 為 Euclidean norm 且 $\lambda \geq 0$為控制力的權重。則現在我們可以給出最佳控制問題如下
\[\min J\left( {x\left( n \right),u\left( . \right)} \right): = \sum\limits_{k = 0}^{N - 1} {l\left( {{x_u}\left( k \right),u\left( k \right)} \right)} \]對 $u(0),...,u(N-1)$ 為admissible 且 $x_u$ 滿足 (**)。

讓我們假設上述最佳控制問題具有解,且此解由一組最佳控制數列 記作
\[
u^*(0),...,u^*(N-1)
\]亦即我們有
\[\mathop {\min }\limits_{u\left( 0 \right),...,u\left( {N - 1} \right)} J\left( {x\left( n \right),u\left( . \right)} \right) = \sum\limits_{k = 0}^{N - 1} {l\left( {{x_{{u^*}}}\left( k \right),{u^*}\left( k \right)} \right)} \]現在為了讓控制力 $u$ 具有我們需要的回授形式 $\mu(x(n))$ 我們令
\[
\mu(x(n) ) := u^*(0)
\]亦即我們只使用最佳控制序列的第一元素。接著我們重複上述過程:在時刻 $t_{n+1}$ 可量得 $x(n+1)$ 並依此再度執行最佳化求得 $\mu(x(n+1))$。同理,對於時刻 $t_{n+2}$ 可量得 $x_{n+2}$ 並依此執行最佳化得到 $\mu(x(n+2))$...


Comments:
1. 因為受控系統為非線性,一般而言為我們所得到的 回授控制律 $\mu(\cdot)$ 必須透過迭代最佳化演算 而得,沒有解析解,故計算複雜度會是模型預測控制的一大挑戰。
2. 就預測模型觀點而言,我們可看出 預測狀態軌跡 $x_u(k), \;\; k=0,1,2,...,N$ 提供了對原本系統在已知時刻 $t_n$ 對於 $t_n,...,t_{n+N}$ 的預測,以及在時刻 $t_{n+1}$ 對於 $t_{n+1},...,t_{n+N+1}$ 的預測,以及 在時刻 $t_{n+2}$ 對於 $t_{n+2},...,t_{n+N+2}$ 的預測。因此我們可看出預測區間是 "移動" 的,故模型預測控制一般又稱為 移動區間控制 (Moving Horizon Control, MHC) 或稱 (Receding Horizon Control, RHC)
3. 若上述 $f$ 為線性,則上述討論稱之為  線性模型預測控制(Linear Model Predictive Control, LMPC)或者簡稱 模型預測控制 (Model Predictive Control, MPC)


2017年6月17日 星期六

[最佳化] 無窮維 Weierstrass 極值定理

===================
Infinite Dimensional Weierstrass Extreme Theorem: 令 $X$ 為 normed vector space 且 $K \subset X$ 為 compact set 。若 $f$ 為 upper semicontinuous on $K$ 則 $f$ 在 $K$上可達到最大值,亦即 存在 $x \in K$ 使得
\[
f(x) = \sup_{z \in K} f(z)
\]=================

Proof: 要證明 $f$ 在 $K$上可達到最大值,亦即 存在 $x \in K$ 使得
\[
f(x) = \sup_{z \in K} f(z)
\]上式等價為
\[
f(x) \geq \sup_{z \in K} f(z) \;\;\; \text{and } \; f(x) \leq \sup_{z \in K} f(z)
\] 注意到 $f(x) \leq \sup_{z\in K} f(z)$ 為顯然,故以下僅需證明\[
f(x) \geq \sup_{z \in K} f(z)
\]

令 $M:= \sup_{z \in K} f(z)$ 則由 supremum 性質可知存在數列 $\{x_n\} \subset K$ 使得
\[
f(x_n) \to M \;\;\;\; (*)
\] 由於 $K$ 為 compact 故必定存在 $\{x_n\} $ 的子數列 記作 $\{x_{n_k}\} \subset K$  使得其收斂在  $x \in K$ ,另外由於 式子 $(*)$ ,我們可推知子數列亦滿足
\[
f( x_{n_k} ) \to M \;\;\;\; \text{ as $k \to \infty$}
 \]由於 $f$ 為 upper semicontinuous on $K$,可知
\[
\limsup_{z \to x} f(z) \leq f(x)
\]由 $\limsup$ 的 數列與函數數列性質 ,我們有
\[
\limsup_{k \to \infty} f(x_{n_k}) \leq \limsup_{z \to x} f(z)  \;\;\;\; (\star)
\]又因為 $f(x_{n_k}) \to M$ (as $k \to \infty$) 故
\[
\limsup_{k \to \infty} f(x_{n_k}) = \lim_{k \to \infty} f(x_{n_k}) = M \;\;\;\; (**)
\]由 $(\star)$ 與 $(**)$,我們有
\[
M=\lim_{k \to \infty} f(x_{n_k}) \leq f(x)
\]即為所求 $\square$

2017年5月12日 星期五

[泛函分析] Bessel's Inequality

下列 Bessel's inequality 在 泛函分析 與 Fourier 分析中扮演重要角色,此不等式表示給定任意在 Hilbert Space 中的點 $x$ (e.g., $L^2$ 空間上 封閉子集的函數), 且給定一組 Hilbert Space 上的正交基底函數 (e.g., complex exponential function, $\{e_i\}$) 則 任意點 $x$ 透過 $\{e_i\}$ 作為基底展開的係數平方和 有上界 且此上界剛好為 $||x||^2$。此不等式的證明要求對內積的性質有進一步掌握,個人認為是很好的練習。

======================
Theorem: Bessel's Inequality
令 $H$ 為 Hilbert Space 且 $x \in H$。若 $\{e_i\} \subset H $ 為一組 orthonormal sequence 則
\[
\sum_{i=1}^\infty | \langle x,e_i \rangle |^2 \leq ||x||^2
\]其中 $\langle \cdot, \cdot \rangle$ 為 $H$ 上的內積運算。
======================


Proof:
令 $\{e_i\} \subset H $ 為一組 orthonormal sequence 且 $x \in H$,現在觀察 $x$ 與 部分和$\sum_{i=1}^n \langle x, e_i \rangle$ 的差異 (透過內積):
\begin{align*}
  0 \leqslant {\left\| {x - \sum\limits_{i = 1}^n {\langle x,{e_i}\rangle {e_i}} } \right\|^2} &= \left\langle {x - \sum\limits_{i = 1}^n {\langle x,{e_i}\rangle } {e_i},x - \sum\limits_{i = 1}^n {\langle x,{e_i}\rangle {e_i}} } \right\rangle  \hfill \\
   &= \left\langle {x,x} \right\rangle  - \left\langle {\sum\limits_{i = 1}^n {\langle x,{e_i}\rangle } {e_i},x} \right\rangle  - \left\langle {x,\sum\limits_{i = 1}^n {\langle x,{e_i}\rangle {e_i}} } \right\rangle  \\
& \hspace{15mm}+ \left\langle {\sum\limits_{i = 1}^n {\langle x,{e_i}\rangle } {e_i},\sum\limits_{j = 1}^n {\langle x,{e_j}\rangle {e_j}} } \right\rangle  \hfill \\
  & = {\left\| x \right\|^2} - \sum\limits_{i = 1}^n {\langle x,{e_i}\rangle } \left\langle {{e_i},x} \right\rangle  - \sum\limits_{i = 1}^n {\overline {\langle x,{e_i}\rangle } } \left\langle {x,{e_i}} \right\rangle  \\
& \hspace{15mm}+ \sum\limits_{i = 1}^n {\langle x,{e_i}\rangle } \sum\limits_{j = 1}^n {\overline {\langle x,{e_j}\rangle } } \left\langle {{e_i},{e_j}} \right\rangle  \hfill \\
   &= {\left\| x \right\|^2} - \sum\limits_{i = 1}^n {{{\left| {\langle x,{e_i}\rangle } \right|}^2}}  - \sum\limits_{i = 1}^n {{{\left| {\langle x,{e_i}\rangle } \right|}^2}}  + \sum\limits_{i = 1}^n {{{\left| {\langle x,{e_i}\rangle } \right|}^2}}  \hfill \\
\end{align*}  上述第三條等式成立是因為應用了 一些內積的 FACT (請參閱下文),另外最後一條等式成立是因為應用了 $\{e_i\}$ 為 orthonormal 的性質 (亦即 $(e_i,e_j) = 1$ 當 $i=j$,且當 $i \neq j$時,$(e_i,e_j)=0$) 故上式最後一項滿足
\[\left( {\sum\limits_{i = 1}^n {\langle x,{e_i}\rangle } \sum\limits_{j = 1}^n {\overline {\langle x,{e_j}\rangle } } } \right)\left\langle {{e_i},{e_j}} \right\rangle  = \sum\limits_{i = 1}^n {{{\left| {\langle x,{e_i}\rangle } \right|}^2}} \]至此,我們得到
\[0 \leqslant {\left\| {x - \sum\limits_{i = 1}^n {\langle x,{e_i}\rangle {e_i}} } \right\|^2} = {\left\| x \right\|^2} - \sum\limits_{i = 1}^n {{{\left| {\langle x,{e_i}\rangle } \right|}^2}} \]亦即
\[\sum\limits_{i = 1}^n {{{\left| {\langle x,{e_i}\rangle } \right|}^2}}  \leqslant {\left\| x \right\|^2}\]注意到此式對任意 $n \in \mathbb{N}$ 皆成立,故取極限亦成立,
\[\sum\limits_{i = 1}^\infty  {{{\left| {\langle x,{e_i}\rangle } \right|}^2}}  \leqslant {\left\| x \right\|^2}\]至此得證。$\square$



Comments:
上述結果保證
\[
\sum_{i=1}^\infty | \langle x,e_i \rangle |^2  <\infty
\]



=======================
FACT: 內積運算的一些常用性質
令 $H$ 為 Hilbert Space,取 $(e_1,e_2,...,e_n)$ 為一組向量 滿足 $e_i \in H$ 且 $x \in H$, $a_i,b_i \in \mathbb{C}$ 則我們有以下一些內積性質:
\[\left\{ \begin{gathered}
  \left\langle {\sum\limits_{i = 1}^n {{b_i}} {e_i},x} \right\rangle  = \sum\limits_{i = 1}^n {{b_i}} \left\langle {x,{e_i}} \right\rangle  \hfill \\
  \left\langle {x,\sum\limits_{i = 1}^n {{b_i}} {e_i}} \right\rangle  = \sum\limits_{i = 1}^n {\overline {{b_i}} } \left\langle {x,{e_i}} \right\rangle  \hfill \\
  \left\langle {\sum\limits_{i = 1}^n {{a_i}} {e_i},\sum\limits_{j = 1}^n {{b_j}} {e_j}} \right\rangle  = \sum\limits_{i = 1}^n {{a_i}} \sum\limits_{j = 1}^n {\overline {{b_j}} } \left\langle {{e_i},{e_j}} \right\rangle  \hfill \\
\end{gathered}  \right.\]其中 $ {\bar a}$ 表示 對 $a$ 取 complex conjugate。
========================

2017年5月7日 星期日

[數理統計] 一致性估計 與 弱大數法則

以下討論 點估計理論 中一些比較重要的性質與應用。

令 $\Theta$ 為某參數空間。

===================
Definition: (Point) Estimator and Estimate
給定 $X_1,...,X_n$ 為 i.i.d. 隨機試驗 來自 pdf $f(x;\theta)$ 其中 $\theta \in \Theta$ 為未知參數,現在定義新的隨機變數 $\widehat \theta$ 為  $X_1,...,X_n$ 的函數,寫為
\[
\widehat \theta := \widehat \theta(X_1,...,X_n)
\]則我們稱此 $\widehat \theta$ 為用以估計參數  $\theta$ 的估計量 (estimator of $\theta$),且若我們能取得隨機樣本的觀察值,比如 $X_1 = x_1,X_2=x_2,...,X_n = x_n$ 則我們稱
\[
\widehat \theta := \widehat \theta(x_1,...,x_n)
\]為參數 $\theta$ 的估計值(estimate)
===================

Comments:
1. 參數 $\theta$ 為分配中的未知常數,但估計量 $\widehat \theta$ 為隨機變數
2. 在數理統計中的 隨機取樣 與 機率論中 iid 隨機變數 視為等價敘述。
3. $\widehat \theta$為  $X_1,...,X_n$ 的函數,但與 $\theta$ 無關!在數理統計中稱此與待估計參數無關的性質為 統計量 (static)。



那麼怎樣的估計量才算是 "好" 的估計量?以下給出一些常見的 "好" 估計量定義。


===================
Definition: What "Good" Properties Should an Estimator Hold? 
令 \[
\widehat \theta := \widehat \theta(X_1,...,X_n)
\] 為某未知參數 $\theta$ 的估計量,則
我們稱 $\widehat \theta $ 為 不偏估計量(unbiased estimator) 若 $E[\widehat \theta] = \theta$
我們稱 $\widehat \theta$ 為 漸進不偏估計量 (asymptotic unbiased estimator) 若 $E[\widehat \theta] \to \theta$ 當 $n \to \infty$
我們稱 $\widehat \theta$ 為 一致估計量 (consistent estimator) 若 \[\widehat \theta \mathop  \to \limits^P \theta
\]亦即給定任意 $\varepsilon >0$ 我們有
\[\mathop {\lim }\limits_{n \to \infty } P\left( {\left| {\hat \theta  - \theta } \right| \geqslant \varepsilon } \right) = 0\]
===================

Comments:
1. 不偏估計在直覺上表示平均而言,估計量 $\widehat \theta$ 不會 高估 或者 低估 參數 $\theta$。
2. 一致性估計則表示隨機試驗夠多之後,對參數的估計量 $\widehat \theta$ 與 真實參數 $\theta$ 幾乎沒有差別。


以下我們簡介一個非常重要的點估計結果定理:若估計量為 不偏  或者漸進不偏,且 估計量 $$的變異收斂到 $0$ 則保證此估計量為一致估計量。

===================
Theorem: Consistency Estimator (Sufficient Condition for Probability Convergent Estimator)
令 $\widehat \theta := \widehat \theta(X_1,...,X_n)$ 為 $\theta \in \Theta$ 的 estimator,若 $E[\widehat \theta] = \theta$ 或 $E[\widehat{\theta}] \to \theta$ 且 $Var(\widehat{\theta}) \to 0$ 當 $n \to \infty$ 則
\[\widehat \theta \mathop  \to \limits^P \theta
\]===================

Proof:
要證明 $\widehat \theta \mathop  \to \limits^P \theta $ 由機率收斂定義可知,給定 $\varepsilon >0$ 我們要證明
\[\mathop {\lim }\limits_{n \to \infty } P\left( {\left| {\hat \theta  - \theta } \right| \geqslant \varepsilon } \right) = 0\]故首先觀察
\[
P\left( {\left| {\hat \theta  - \theta } \right| \geqslant \varepsilon } \right) = P\left( {{{\left( {\hat \theta  - \theta } \right)}^2} \geqslant {\varepsilon ^2}} \right) \leqslant \frac{{E\left[ {{{\left( {\hat \theta  - \theta } \right)}^2}} \right]}}{{{\varepsilon ^2}}}\;\;\;\;\;\;\; (*)
\]上述不等式成立是因為利用了 Generalized Chebyshev's inequality。現在我們觀察
 \begin{align*}
  E \left[ {{{\left( {\widehat \theta  - \theta } \right)}^2}} \right]
&= E\left[ {{{\left( {\widehat \theta  - E\left[ {\widehat \theta } \right] + E\left[ {\widehat \theta } \right] - \theta } \right)}^2}} \right] \hfill \\
  & = E\left[ {{{\left( {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right)}^2} + 2\left( {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right)\left( {E\left[ {\widehat \theta } \right] - \theta } \right) + {{\left( {E\left[ {\widehat \theta } \right] - \theta } \right)}^2}} \right] \hfill \\
  & = \underbrace {E\left[ {{{\left( {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right)}^2}} \right]}_{ = Var\left( {\widehat \theta } \right)} + \underbrace {2E\left[ {\left( {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right)\left( {E\left[ {\widehat \theta } \right] - \theta } \right)} \right]}_{ = 2\left( {E\left[ {\widehat \theta } \right] - \theta } \right)E\left[ {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right]} \\
&\;\;\;\;\; + \underbrace {E\left[ {{{\left( {E\left[ {\widehat \theta } \right] - \theta } \right)}^2}} \right]}_{ = {{\left( {E\left[ {\widehat \theta } \right] - \theta } \right)}^2}} \hfill \\
 &  = Var\left( {\widehat \theta } \right) + 2\left( {E\left[ {\widehat \theta } \right] - \theta } \right)E\left[ {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right] + {\left( {E\left[ {\widehat \theta } \right] - \theta } \right)^2} \hfill \\
\end{align*}
上述結果用到了 $E[\widehat{\theta}] = constant$ 故
\[\left\{ \begin{align*}
  &E\left[ {{{\left( {E\left[ {\widehat \theta } \right] - \theta } \right)}^2}} \right] = {\left( {E\left[ {\widehat \theta } \right] - \theta } \right)^2}; \hfill \\
  &2E\left[ {\left( {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right)\left( {E\left[ {\widehat \theta } \right] - \theta } \right)} \right] = 2\left( {E\left[ {\widehat \theta } \right] - \theta } \right)\underbrace {E\left[ {\widehat \theta  - E\left[ {\widehat \theta } \right]} \right]}_{ = 0} = 0 \hfill \\
\end{align*}  \right.\]故
\[E\left[ {{{\left( {\widehat \theta  - \theta } \right)}^2}} \right] = Var\left( {\widehat \theta } \right) + {\left( {E\left[ {\widehat \theta } \right] - \theta } \right)^2}\]現在將此結果代入 $(*)$ ,可得
\begin{align*}
  P\left( {\left| {\widehat \theta  - \theta } \right| \geqslant \varepsilon } \right) &\leqslant \frac{{E\left[ {{{\left( {\widehat \theta  - \theta } \right)}^2}} \right]}}{{{\varepsilon ^2}}} \hfill \\
  & = \frac{{Var\left( {\widehat \theta } \right) + {{\left( {E\left[ {\widehat \theta } \right] - \theta } \right)}^2}}}{{{\varepsilon ^2}}} \to \frac{1}{\varepsilon^2 }\left( {0 + 0} \right) = 0 \hfill \\
\end{align*}
又因為  機率測度恆正,我們有 $0 \leq P\left( {\left| {\widehat \theta  - \theta } \right| \geqslant \varepsilon } \right)  \to 0 $ 由極限的夾擊定理可知
\[\mathop {\lim }\limits_{n \to \infty } P\left( {\left| {\widehat \theta  - \theta } \right| \geqslant \varepsilon } \right) = 0\]故此得證。 $\square$


Comments:
在數理統計中,一致估計 與 機率論中的 機率收斂為等價。


上述定理可用來非常快速的證明弱大數法則 WLLN:亦即用 一組隨機試驗的 樣本平均數 作為 估計量來估計平均可以保證 機率收斂。

=============
Theorem 2: Weak Law of Large Numbers
令 $X_1,X_2,...$ 為 i.i.d. 隨機變數 數列 且具有共同 $E[X_1] = \mu$ 與變異數 $Var(X_1) = \sigma^2 < \infty$ 現在令
\[
\bar X := \frac{1}{n} \sum_{i=1}^n X_i
\]則
\[
\bar X \mathop  \to \limits^P \mu
\]=============

Proof:
此結果可利用前述 Consistency Estimator Theorem 求證,也就是把 $\widehat \theta := \bar X$ 且 $\theta := \mu$。則我們僅需證明 $\bar X$ 為不偏估計 ($E[\bar X] = \mu$) 或者 漸進不偏估計 ($E[\bar X] \to \mu$),且$\bar X$變異數收斂到 $0$ 即可。故現在觀察
\begin{align*}
  E\left[ {\bar X} \right] &= E\left[ {\frac{1}{n}\sum\limits_{i = 1}^n {{X_i}} } \right] \hfill \\
   &= \frac{1}{n}E\left[ {\sum\limits_{i = 1}^n {{X_i}} } \right] \hfill \\
  & = \frac{1}{n}\sum\limits_{i = 1}^n {E\left[ {{X_i}} \right]}  \hfill \\
  &\mathop  = \limits^{X_i \; i.i.d.} \frac{1}{n}\sum\limits_{i = 1}^n {E\left[ {{X_1}} \right]}  \hfill \\
   &= \frac{1}{n}nE\left[ {{X_1}} \right] = E\left[ {{X_1}} \right] = \mu  \hfill \\
\end{align*} 故我們知道 $\bar X$ 為不偏估計量,現在我們僅需證明 變異數收斂到 $0$ 。觀察
\begin{align*}
  Var\left( {\bar X} \right) &= Var\left( {\frac{1}{n}\sum\limits_{i = 1}^n {{X_i}} } \right) \hfill \\
  &\mathop  = \limits^{X_i \; i.i.d.} \frac{1}{{{n^2}}}\sum\limits_{i = 1}^n {Var\left( {{X_i}} \right)}  \hfill \\
  & = \frac{1}{{{n^2}}}nVar\left[ {{X_1}} \right] \hfill \\
 &  = \frac{1}{{{n^2}}}n{\sigma ^2} = \frac{1}{n}{\sigma ^2} \to 0 \hfill \\
\end{align*} 故由 consistency theorem ,
\[
\bar X \mathop  \to \limits^P \mu
\]至此得證。 $\square$

2017年5月2日 星期二

[機率論] Chebyshev's Inequality 的推廣型

以下介紹一個在機率論中 相當有用的 一條不等式,稱為 Chebyshev inequality,此不等式將 期望值 與 機率測度 做出一定程度的連結 來用以估計 期望值的下界 (或者說 求某機率測度的上界)。以下我們給出此不等式之陳述與證明:讀者可注意要求的假設條件並不多,證明也稍具巧思。

================
Theorem: Generalized Chebyshev's Inequality 
令 $X$ 為 任意 連續型 隨機變數 配備 機率密度函數 $f_X$ ,現在定義 $g(X)$ 為任意非負函數,若 $E[g(X)]$ 存在,則 對任意常數 $c>0$,我們有
\[
P(g(X) \geq c) \leq \frac{E[g(X)]}{c}
\]================


Proof: 假設 $X$ 為連續型隨機變數且 $E[g(X)]$ 存在,令 $c >0$ 為任意正值常數。由於期望值 $E[g(X)]$ 存在,由定義可知我們有
\[
E\left[ {g\left( X \right)} \right] = \int_{ - \infty }^\infty  {g\left( x \right){f_X}\left( x \right)dx}
\]其中 $f_X(x)$ 為 $X$ 的 機率密度函數 (Probability Density Function, pdf)。現在觀察上述右式積分,我們可將其等價寫為
\begin{align*}
  \int_{ - \infty }^\infty  {g\left( x \right){f_X}\left( x \right)dx}  &= \int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {g\left( x \right){f_X}\left( x \right)dx}  + \int_{\left\{ {x:g\left( x \right) < c} \right\}}^{} {g\left( x \right){f_X}\left( x \right)dx}  \hfill \\
 &  \geq \int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {g\left( x \right){f_X}\left( x \right)dx}  \;\;\;\; (*)
\end{align*} 注意到上述不等式成立 是 因為對所有 $x$ 而言, $g(x) \geq 0$ 且 pdf $f_X(x) \geq 0$。現在我們觀察不等式右方的積分式子 $ \int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {g\left( x \right){f_X}\left( x \right)dx} $ 可以發現此積分範圍是對所有的 $x$ 滿足 $g(x) \geq c$,這表示我們可以進一步寫出此積分的下界
\[\int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {g\left( x \right){f_X}\left( x \right)dx}  \geqslant \int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {c{f_X}\left( x \right)dx} \;\;\;\;\; (**)
\]由 $(*)$ 與 $(**)$ 我們可知
\begin{align*}
  E\left[ {g\left( X \right)} \right] &\geq \int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {c{f_X}\left( x \right)dx}  \hfill \\
   &= c\int_{\left\{ {x:g\left( x \right) \geqslant c} \right\}}^{} {{f_X}\left( x \right)dx}  \hfill \\
   &= cP\left( {g\left( x \right) \geqslant c} \right) \hfill \\
\end{align*} 故
\[\frac{{E\left[ {g\left( X \right)} \right]}}{c} \geqslant P\left( {g\left( x \right) \geqslant c} \right)\]至此得證。$\square$

Comments:
1. 上述定理對 離散型隨機變數仍然成立,僅須將證明的積分部分 $\int (\cdot)$ 改成累加 $\sum (\cdot)$ 即可。
2. Chebyshev's inequality 的界的 "鬆緊程度" 依隨機變數情況而定,故拿來做精準上下界估計不一定準確。
3. 若 $g(X):=X$ 則上述 廣義 Chebyshev's inequality 又稱作 Markov's inequality。
4. 上述廣義的 Chebyshev's inequality 應用在於如何"識別"或者 "選取" 適當的 非負函數 $g(X)$,我們會在以下再作進一步說明。



以下我們看個上述定理的應用例:

==================
FACT 1: Standard Chebyshev's inequality 
令 $X$ 為隨機變數具有 有限期望值 與變異數,記作 $ E[X] := \mu$ 且 $E[(X- \mu)^2] =\sigma^2$ 。則對任意 $n>0$ 而言,我們有
\[
P( (X-\mu)^2 \geq n^2 \sigma^2) \leq \frac{1}{n^2}
\]==================

Proof:
給定 $n >0$, 定義 $g(X) := (X-\mu)^2 \geq 0$,則 Generalized Chebyshev's inequality 告訴我們對任意 $c>0$,我們有
\[P({\left( {X - \mu } \right)^2} \geqslant c) \leqslant \frac{{E\left[ {{{\left( {X - \mu } \right)}^2}} \right]}}{c}\]又因為 ${E\left[ {{{\left( {X - \mu } \right)}^2}} \right] = {\sigma ^2}}$ 故上式可改寫為
\[
P({\left( {X - \mu } \right)^2} \geqslant c) \leqslant \frac{{{\sigma ^2}}}{c}
\]現在取 $c:=\sigma^2 n^2 >0$ 則\[
P({\left( {X - \mu } \right)^2} \geqslant n^2 \sigma^2 ) \leqslant \frac{1}{n^2}
\]至此得證。$\square$

==================
FACT 2:
令 $X$ 為隨機變數配備 期望值 $\mu$ 且令 $E[(X-\mu)^{2k}]$ 對任意正整數 $k$ 都存在,則對任意 $c >0$,
\[
P(|X-\mu| \geq c) \leq E[(X-\mu)^{2k}]/c^{2k}
\]==================

Proof: omitted. (取 $g(X):= (X-\mu)^{2k}$ 且 $c = d^{1/2k}, \;\; \forall d>0$ )


以下結果為利用 Chebyshev inequality 與 動差生成函數 Moment Generating Function (mgf) 拉上關係:


==================
FACT 3:
令 $X$ 為隨機變數配備 mgf 滿足下列條件:存在 $\delta>0$ 使得當 $t \in (-\delta, \delta)$,其 mgf $M_X(t)$ 存在,則
\[
P(X \geq c) \leq e^{-ct} M_X(t), \;\;\; t \in (0, \delta)
\]且
\[
P(X \leq c) \leq e^{-ct} M_X(t), \;\;\; t \in (-\delta,0)
\]==================

Proof: omitted (取 $g(X):= e^{tX}$ 且 $c = \frac{\log d}{ t}, \;\; \forall d>0$  )


2017年4月30日 星期日

[機率論] 動差生成函數 的 常見應用 (1)

令 $\Omega$ 為 樣本空間,定義 $X : \Omega \to \mathbb{R}$ 為配備 機率分配函數 $f_X$ 的 隨機變數。

==================
Definition: k-th Moment
隨機變數 $X : \Omega \to \mathbb{R}$  的 k階 動差 (k-th moment) 定為 $E[X^k]$
==================


==================
Definition: Moment Generating Function
令 $X$ 為隨機變數,若存在 $\delta >0$ 使得對 $t \in (-\delta,\delta )$ 而言,期望值 $E[e^{tX}]$ 存在,則 $X$ 的 動差生成函數 (Moment Generating Function, mgf) 存在,且定義為
\[
M_X(t) := E[e^{tX}], \;\;\; t \in (-\delta,\delta )
\]除此之外,若上述條件成立則
\[
D_t^k E[e^{tX}] = E[D_t^k e^{tX}]
\]其中 $D_t^k $ 表示為對 $t$ 微分 $k$次 微分算子。
==================

Comments:
1. 動差生成函數 $M_X(t)$ 是一個以 $t$ 為變數 的函數,目的在於 "產生動差",至於如何產生我們會在下面進行討論。
2. 上述定義僅僅要求 mgf 在 $t = 0$ 附近開區間 $t \in (-\delta, \delta)$ 期望值存在,此條件也保證積分與微分互換性。
3. 若在開區間 $ t \in (-\delta, \delta)$ 期望值存在,立刻可得知 $M_X(0) = E[e^{0}] = 1$


FACT: 上述動差生成函數算是非常便利的工具,我們可以透過其產生各種具有常見分配的隨機變數之一,二階動差。假定某隨機變數之 mgf 存在,則此隨機變數的 k-th 動差表為
\[
D_t^k M_X (0) = E[X^k], \;\; k=1,2,...
\]

Comment:
由上述討論可知,若我們想求期望值則
\[
D_t M_X(0) = E[X]
\]若我們想求變異數則
\begin{align*}
Var(X) &= E[(X - E[X] )^2]\\
& = E[X^2] -(E[X])^2\\
& = D_t^2 M_X(0) - D_t M_X(0)
\end{align*}

==================
Theorem: MGF唯一決定分配函數且 兩 MGF 相等 保證 分配函數相等 
令 $X,Y$ 為兩隨機變數且其各自對應的 mgf $M_X, M_Y$ 在含 $0$ 開區間存在。則此兩隨機變數的分配函數 $$
F_X(z) = F_Y(z), \;\;\; \forall \;\; z \in \mathbb{R}
$$若且唯若 存在 $\delta >0$ 使得 對任意 $t \in (-\delta, \delta)$ 而言
$$M_X(t) = M_Y(t)$$
==================
Proof: omitted


Comments:
1. 上述定理是非常強大的結果,他說明了要決定某隨機變數的分配可以透過 mgf 來唯一決定。也就是說 mgf 與 分配函數為 1-1 對應,且如果當分配函數很困難取得,我們可以透過計算 mgf 來幫助我們確定分配。
2. 機率論中另外有一種函數稱作 特性函數 (characteristic function) 其作用與 mgf相仿,且永遠存在,定義為 $\varphi_X(t) := E[e^{i tX}]$,但是使式子蘊含複數,一般在計算上會稍微比 mgf 複雜一些。


以下我們展示幾個例子來使用動差生成函數求得 期望值 與 變異數。首先我們看幾個離散隨機變數的例子:

==================
Example 1: Bernoulli Random Variable 的 期望值 與 變異數
令 $X$ 為 Bernoulli random variable 滿足 $P(X=1) = p$ 且 $P(X=0)=1-p$。
(a) 試問其 mgf 是否存在?若存在則求其 mgf
(b) 利用 part(a) 證明 $E[X] = p$ 與 $Var(X) = p (1-p)$。
==================

Proof (a): 以下我們使用 mgf 方法來求期望值 以及 變異數 ,首先計算 mgf ,由定義可知
\begin{align*}
  {M_X}(t) &= E\left[ {{e^{tX}}} \right] \hfill \\
   &= \sum\limits_{x \in \left\{ {0,1} \right\}}^{} {{e^{tx}}{f_X}\left( x \right)}  \hfill \\
   &= {e^{t1}}p + {e^{t0}}\left( {1 - p} \right) \hfill \\
   &= {e^t}p + \left( {1 - p} \right) \hfill \\
\end{align*} 注意到上式對任意 $t \in \mathbb{R}$ 成立,故存在 $\delta >0$ 使得對 $t \in (-\delta,\delta )$ 而言,期望值 $E[e^{tX}]$ 存在。

Proof (b): 
以下我們計算 mgf 的微分
\[\left\{ \begin{align*}
  {D_t}{M_X}(t) &= {D_t}\left( {{e^t}p + \left( {1 - p} \right)} \right) = {e^t}p \hfill \\
  D_t^2{M_X}(t) &= {e^t}p \hfill \\
\end{align*}  \right.\] 故 $ E[X] = {D_t}{M_X}(0) = p $ 且
\begin{align*}
 Var(X) &= D_t^2{M_X}(0) - {({D_t}{M_X}(0))^2} \hfill \\
 &= p - {p^2} = p\left( {1 - p} \right) \hfill \\
 \end{align*}至此得證 $\square$





==================
Example 2: Poisson Random Variable 的 期望值 與 變異數
令 $X \sim Poisson(\lambda)$ 為 Poisson random variable  配備機率密度函數
\[
f_X(x) = \frac{\lambda^x e^{-\lambda}}{x!},\;\;\; x=0,1,2,...
\]
(a) 試問其 mgf 是否存在?若存在則求其 mgf
(b) 利用 part(a) 求 $E[X] = \lambda$ 與 $Var(X) = \lambda$。
==================

Proof (a):
由 mgr 定義,我們觀察 \begin{align*}
  {M_X}\left( t \right) &= E\left[ {{e^{tX}}} \right] \hfill \\
   &= \sum\limits_{x = 0}^\infty  {{e^{tx}}{f_X}(x)}  \hfill \\
   &= \sum\limits_{x = 0}^\infty  {{e^{tx}}\frac{{{\lambda ^x}{e^{ - \lambda }}}}{{x!}}}  \hfill \\
   &= {e^{ - \lambda }}\underbrace {\sum\limits_{x = 0}^\infty  {\frac{{{{\left( {{e^t}\lambda } \right)}^x}}}{{x!}}} }_{ = {e^{{e^t}\lambda }}} = {e^{\left( {{e^t} - 1} \right)\lambda }} \hfill \\
\end{align*}
上式最後一條等式成立因為
\[{e^z}: = \sum\limits_{k = 0}^\infty  {\frac{{{z^k}}}{{k!}}},\;\;\; \forall \; z \in \mathbb{C} \]注意到 $M_X(t)$ 對任意 $t \in \mathbb{R}$ 有定義,故存在 $\delta >0$ 使得對 $t \in (-\delta,\delta )$ 而言,期望值 $E[e^{tX}]$ 存在。

Proof (b):
首先對 mgf 求一階 與 二階導數,
\[\left\{ \begin{gathered}
  {D_t}{M_X}\left( t \right) = {D_t}\left( {{e^{\left( {{e^t} - 1} \right)\lambda }}} \right) = {e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t} \hfill \\
  D_t^2{M_X}\left( t \right) = D_t^2\left( {{e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t}} \right) = {e^{\left( {{e^t} - 1} \right)\lambda }}{\left( {\lambda {e^t}} \right)^2} + {e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t} \hfill \\
\end{gathered}  \right.\] 故期望值為 $E\left[ X \right] = {D_t}{M_X}\left( 0 \right) = \lambda $ 且變異數為
\begin{align*}
  Var\left( X \right) &= D_t^2{M_X}\left( 0 \right) - {\left( {{D_t}{M_X}\left( 0 \right)} \right)^2} \hfill \\
   &= {\left. {\left( {{e^{\left( {{e^t} - 1} \right)\lambda }}{{\left( {\lambda {e^t}} \right)}^2} + {e^{\left( {{e^t} - 1} \right)\lambda }}\lambda {e^t}} \right)} \right|_{t = 0}} - {\lambda ^2} = \lambda  \hfill \\
\end{align*} 至此得證 $\square$






下面例子是 連續隨機變數 的情況。

=======================
Example 3: Normal Distribution
令 $X \sim \mathcal{N}(\mu, \sigma^2)$ 配備機率密度函數
\[{f_X}(x): = \frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( { - \frac{{{{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}} \right), \;\;\; x \in \mathbb{R}
\]其中 $\mu \in \mathbb{R}, \sigma >0$
(a) 試求 $X$ 的 mgf
(b) 利用 part(a) 證明 $E[X] =\mu$ 與 $Var(X) =\sigma^2$
=======================

Proof: (a)
由定義出發,觀察
\begin{align*}
  {M_X}\left( t \right) &= E\left[ {{e^{tX}}} \right] \hfill \\
   &= \int_{ - \infty }^\infty  {{e^{tx}}{f_X}(x)} dx \hfill \\
   &= \int_{ - \infty }^\infty  {{e^{tx}}\frac{1}{{\sqrt {2\pi } \sigma }}{e^{ - \frac{{{{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= \frac{1}{{\sqrt {2\pi } \sigma }}\int_{ - \infty }^\infty  {{e^{\frac{{2{\sigma ^2}tx - {{\left( {x - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= \frac{1}{{\sqrt {2\pi } \sigma }}{e^{\frac{{ - {\mu ^2}}}{{2{\sigma ^2}}}}}\int_{ - \infty }^\infty  {{e^{\frac{{ - \left( {{x^2} - 2\left( {{\sigma ^2}t + \mu } \right)x + {{\left( {{\sigma ^2}t + \mu } \right)}^2}} \right)}}{{2{\sigma ^2}}}}}{e^{\frac{{{{\left( {{\sigma ^2}t + \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= \frac{1}{{\sqrt {2\pi } \sigma }}{e^{\frac{{ - {\mu ^2}}}{{2{\sigma ^2}}}}}{e^{\frac{{{{\left( {{\sigma ^2}t + \mu } \right)}^2}}}{{2{\sigma ^2}}}}}\int_{ - \infty }^\infty  {{e^{\frac{{ - {{\left( {x - \left( {{\sigma ^2}t + \mu } \right)} \right)}^2}}}{{2{\sigma ^2}}}}}} dx \hfill \\
   &= {e^{\frac{{{\sigma ^2}{t^2} + 2\mu t}}{2}}}\underbrace {\frac{1}{{\sqrt {2\pi } \sigma }}\int_{ - \infty }^\infty  {{e^{\frac{{ - {{\left( {x - \left( {{\sigma ^2}t + \mu } \right)} \right)}^2}}}{{2{\sigma ^2}}}}}} dx}_{ = 1{\text{ }}\left( {{\text{pdf of normal }}\mathcal{N}\left( {\mu+\sigma^2 t ,{\sigma ^2}} \right)} \right)} \hfill \\
   &= {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}} \hfill \\
\end{align*}

Proof (b):
以下我們計算 mgf 的微分
\[\left\{ \begin{align*}
  {D_t}{M_X}\left( t \right) &= {D_t}\left( {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}} \right) = {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}\left( {\mu  + {\sigma ^2}t} \right) \hfill \\
  D_t^2{M_X}\left( t \right) &= {D_t}\left( {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}\left( {\mu  + {\sigma ^2}t} \right)} \right) = {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{\left( {\mu  + {\sigma ^2}t} \right)^2} + {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{\sigma ^2} \hfill \\
\end{align*}  \right.\]故對應的期望值為
\[E\left[ X \right] = {D_t}{M_X}\left( 0 \right) = {\left. {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}\left( {\mu  + {\sigma ^2}t} \right)} \right|_{t = 0}} = \mu \]
且 變異數 為
\begin{align*}
  Var\left( X \right) &= D_t^2{M_X}\left( 0 \right) - {\left( {D_t^1{M_X}\left( 0 \right)} \right)^2} \hfill \\
   &= {\left. {\left( {{e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{{\left( {\mu  + {\sigma ^2}t} \right)}^2} + {e^{\mu t + \frac{{{\sigma ^2}{t^2}}}{2}}}{\sigma ^2}} \right)} \right|_{t = 0}} - {\mu ^2} \hfill \\
   &= {\sigma ^2} \hfill \\
\end{align*} 至此得證 $\square$

2017年4月26日 星期三

[訊號處理] 2d Convolution & 簡單的影像處理 (利用 MATLAB )

一般在處理影像的時候我們會把 影像(image) 視為 二維離散訊號, 更一般的說法是將其視為矩陣,亦即給定任一張 (灰階) 影像 我們可以將其分割成 $M \times N$ 矩陣 (像素 pixel),並且將其記作 $x[m,n]$,其中 $m=0,...,M-1$ 與 $n = 0,...,N-1$ 。

Comments:
如果要處理的影像是彩色的,一個常用的做法是把該影像以 $R,G,B$ 三色分別存成 三個矩陣,最後在做疊合。在此不做贅述。

以下我們利用 MATLAB 來執行 (灰階)影像處理,我們使用的圖檔是 MATLAB內建的圖檔 cameraman.tif ,當然讀者可自行讀入任何自己想要的圖檔。在MATLAB 輸入

MATLAB Code For Loading the Image
x = imread('cameraman.tif' )
imagesc(x)

則會顯示
圖1a: cameraman.tif 原圖 ($256 \times 256$)

Comments:
1. 上述影像訊號 $x[m,n]$ 為 $M \times N = 256 \times 256 $ 矩陣。
2. 一般而言在影像處理領域,我們將圖片最左上角點視為座標原點,對應的影像訊號為 $x[0,0]$ 。
3. 上述影像透過 MATLAB 2017a 讀圖會顯示有一些藍綠等顏色,若要使此圖顯示為灰階(gray scale) 讀者可鍵入 colormap(gray) 來使其變成灰階影像,亦即

MATLAB Code For Gray Scale Colormap
x = imread('cameraman.tif' )
imagesc(x)
colormap(gray)

則我們會得到如下圖

圖1b: 灰階影像



現在我們可對上述影像進行一些常見的基本處理。令 $x[m,n]$ 為輸入影像訊號,且假設 影像 濾波器 為 線性非時變 (Linear Time Invariant, LTI) 故其對應的 impulse response $h[m,n]$ 可以用來完全描述我們的影像濾波器,最後 我們令 影像處理過後的輸出訊號 $y[m,n]$ 可表為 $x[m,n]$ 與 $h[m,n]$ 的 convolution ,差別僅在此時我們的 convolution運算為二維運算,故我們寫成
\[
y[m,n] = h[m,n]**x[m,n]: = \sum\limits_{k = 0}^{M - 1} {\sum\limits_{l = 0}^{N - 1} {h\left[ {k,l} \right]x\left[ {m - k,n - l} \right]} }
\] 其中 $**$ 表示二維 convolution,在 MATLAB 中我們可以使用 conv2 來執行此運算。下圖顯示了上述討論的觀點:

Comments:
1. 上述討論中提及的 脈衝響應 $h[m,n]$ 在影像處理中又被稱為 點擴散函數 (point spread function, PSF)
2. 給定任意 影像 $x[m,n]$ 其影像尺寸為 $M_1 \times N_1$ 且 影像濾波器  $h[m,n]$ 具有 影像尺寸為$M_2 \times N_2$ 則其經過 2d convolution之後的 輸出 $y[m,n]$ 之影像尺寸會變成
$$
(M_1 + M_2 -1) \times (N_1 + N_2 -1)
$$ 故上述 2d convolution 執行之後原圖尺寸會變成
$$
(256+2-1) \times (256 + 2 -1) = 257 \times 257
$$也就是說透過2d convolution之後 影像邊緣 會多跑出一些額外的像素。讀者可參考下方後續的討論。
3. 上述討論中我們假設 影像濾波器為 LTI ,故輸入與輸出關係為  (time-domain) convolution,那麼熟悉訊號與系統的讀者不難做出以下猜想,對輸入 $x$ 與 輸出 $y$ 與 影像濾波器 $h$ 取 Discrete Fourier Transform (DFT) 可得到 對應的 DFT係數 如 \[\begin{gathered}
  X[k,l] = DFT(x[m,n]); \hfill \\
  H[k,l] = DFT(h[m,n]); \hfill \\
  Y[k,l] = DFT(y[m,n]). \hfill \\
\end{gathered} \] 假設 $N,M$ 夠大,則其在 輸入輸出在頻域上對應的關係為 frequency domain  為相乘
$$
Y[k,l] = H[k,l] X[k,l]
$$
一但計算出 $Y[k,l]$ 在對其取 Inverse Digital Fourier Transform (IDFT) 即可得到 $y[m,n]$。在 MATLAB 中,上述討論可以使用 fft2ifft2 來實現,在此不做贅述。


以下我們探討幾種常見的影像濾波器,亦即幾種常見的 $h[.]$:

Image Filtering:
以下為幾種常見的影像濾波器的例子:
1. 影像模糊 (Image blurring)
\[{h_b}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{1/10} \\
  {1/10}&{1/10}
\end{array}} \right];\]

MATLAB code for Image Blurring
h_b = 1/10 * [1 1; 1 1];
y_b = conv2(x, h_b, 'same')
imagesc( abs(y_b) )

圖2: 模糊效果

讀者可以比較此圖與之前的圖不難發現此圖較原圖為模糊。另外可注意到模糊之後的影像邊緣部分出現額外不屬於原圖的像素。


2. 邊緣偵測:圖像的邊緣可以透過特定的濾波器設計來偵測,比如說我們可以考慮
\[\begin{gathered}
  {h_h}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{1/10} \\
  { - 1/10}&{ - 1/10}
\end{array}} \right]; \hfill \\
  {h_v}: = \left[ {\begin{array}{*{20}{c}}
  {1/10}&{ - 1/10} \\
  {1/10}&{ - 1/10}
\end{array}} \right]. \hfill \\
\end{gathered} \]其中 $h_h$ 用以偵測 影像的水平邊緣,$h_v$ 用以偵測影像的垂直邊緣。

MATLAB Code for Horizontal Edge Detection
h_h = 1/10 * [1 1; -1 -1];
y_h = conv2(x, h_h, 'same')
imagesc( abs(y_h) )

我們得到對於影像水平邊緣的偵測如下圖所示
圖3: 水平邊緣偵測


MATLAB Code for Vertical Edge Detection
h_v = 1/10 * [1 -1; 1 -1];
y_v = conv2(x, h_v, 'same')
imagesc( abs(y_v) )

上述code得到對於垂直邊緣的偵測如下圖
圖4: 垂直邊緣偵測

Comments:
1. 讀者可自行改變 $h_v, h_h, h_b$  的矩陣的係數來看看得到的圖有什麼不同。
2. 影像處理有非常多的有趣的問題可以進一步討論,比如說被 模糊後之後的影像是否可以把它還原?如果可以該怎麼做? 直接用 conv2 或者用 fft2 何者運算較快?影像太大該怎麼壓縮?等等。