推 hwanger : 推一般情況的分析 09/08 12:57
引述《ethan0221 (Ethan)》之銘言：
: Derive the distribution (pdf) of the r-th power y of a normally distributed
: random variable x, that is, y=x^r, x ~ N(μ, σ2). Assume r > 0. What happens
: if r = 0 or r < 0? What is E(1/x) if x is normal? Then find the mean of x^r
: using the approximation by the delta method. What is required for the delta
: method approximation to be valid?
分幾個問題:
(1) r > 0 時
(2) r ≦ 0
(3) E[1/X]
(4) Delta method
(1) r > 0 時
首先, x^r 必須能定義狂所有 x in R, 至少是 almost everywhere.
因此, r 必須是整數, 或特定的有理數, 如 k/3, k/5 之類的, 其中
k 是整數. 一般地說, 如果 r 不是整數, 而其最簡分數 q/p 之分母
p 是奇數, 則 x^r 除了 r < 0 時在 x = 0 無定義之外, 都有定義.
h(x) = x^r, x in R 有可能是偶函數也有可能是奇函數.
當 r 以整數或最簡分數表現, 其分子是偶數時, 則 h(x) 是偶函數;
當 r 是奇整數或其最簡分數的分子是奇數時, 則 h(x) 是奇函數.
限制 r > 0 時, h(x) 在 x ≧ 0 部分是嚴格遞增的, 因此, 如果它
是奇函數, 則它在整個定義域 (R) 上是嚴格遞增的, 所以是一對一.
同果 h(x) 是偶函數, 把 R 分 x≧0 和 x<0 兩部分, 則它在這兩部
分分別是一對一的.
X^r 之機率分配推導, 最基本的有分配函數法和 Jacobian 法.
(a) h(x) 是奇函數時,
設 t ≧ 0, 則
P[X^r ≦ t] = P[X ≦ t^(1/r)] = ∫_(-∞,t^(1/r)] f(x) dx
其中 f(x) 是 X 的 p.d.f.
t < 0 則
P[X^r ≦ t] = P[(-X)^r ≧ (-t)] = P[-X ≧ (-t)^(1/r)]
= P[X ≦ -(-t)^(1/r)] = ∫_(-∞,-(-t)^(1/r)] f(x) dx
Jacobian 法是針對 p.d.f. 直接變換:
h 的反函數是
h^(-1)(y) = y^(1/r) if y ≧ 0,
= -(-y)^(1/r) if y < 0.
故其 Jacobian, 在 1 variable 即導數, 是
J(y) = (h^(-1))'(y) = (1/r)|y|^(1/r -1)
所以, X^r 的 p.d.f. 是
g(t) = f(h^(-1)(t))|J(t)|
= {(1/r)t^(1/r-1)/[√(2π)σ]}exp{-(t^(1/r)-μ)^2/(2σ^2)}
if t ≧ 0;
= {(1/r)(-t)^(1/r-1)/[√(2π)σ]}exp{-[-(-t)^(1/r)-μ]^2/(2σ^2)}
if t < 0.
事實上 h(x) 為奇函數時, 依前面所述可知 1/r = p/q, q 是奇數,
因此 h 的反函數可直接寫 h^(-1)(t) = t^(1/r), 且此也是奇函數,
其導函數是偶函數 (除了 0 以外). 故
P[X^r ≦ t] = P[X ≦ t^(1/r)] = ∫_(-∞,t^(1/r) f(x) dx,
其 p.d.f. 也可簡單寫
g(t) = {(1/r)|t|^(1/r-1)/[√(2π)σ]}exp{-(t^(1/r)-μ)^2/(2σ^2)}
all t in R-{0}.
(b) h(x) = x^r 為偶函數時, h(x) ≧ 0, 也就是說 X^r 非負.
P[X^r ≦ t] = P[-t^(1/r) ≦ X ≦ t^(1/r)] = ∫_[-t,t] f(x) dx,
when t ≧ 0,
= 0 when t < 0.
h 的局部反函數是
h^(-1)(t) = t^(1/r) if t = h(x), x ≧ 0;
= -t^(1/r) if t = h(x), x < 0.
所以兩部分 Jacobian 的絕對值形式都是 (1/r)t^(1/r-1). 故
X^r 之 p.d.f. 為
(1/r)t^(1/r-1)f(t^(1/r)) + (1/r)t^(1/r-1)f(-t^(1/r))
= {(1/r)t^(1/r-1)/[√(2π)σ]} ×
([exp{-(t^(1/r)-μ)^2/(2σ^2)}+exp{-(t^(1/r)+μ)^2/(2σ^2)})
(2) r ≦ 0.
當 r = 0 時, h(x) = 1 for all x≠0,
所以, 除非 X = 0 with probability 則 X^r 限入定義困境, 否則
P[X^0 = 1] = 1.
當 r < 0 時, h(x) = x^r 在 x≠0 只要 x^(-r) 有定義則 h(x) 有
定義, 事實上 h(x) = 1/x^(-r). 這表示:
(i) r 必須是負整數或最簡分數為 -q/p, 其中 p 為奇數.
(ii) h(x) (除 x=0 之外) 可能是偶函數或奇函數.
(iii) 在 x > 0 時, h(x) 是嚴格遞減的. x → +∞ 時 h(x)→0.
當 x → 0+ 時 h(x) → +∞. 至於 x < 0 的圖形特性, 根據
h(x) 是奇或偶, 及 h(x) 在 x > 0 之表現. 易知.
欽求 X^r 之機率分配, Jacobian 法變換 p.d.f. 是比較簡單的.
當 h(x) 是奇函數時, 它是一對一的, 其反函數是
h^(-1)(t) = t^(1/r), t≠0
Jacobian 為 (1/r)t^(1/r-1), 故 X^r 之 p.d.f. 為
g(t) = {(-1/r)|t|^(1/r-1)/[√(2π)σ]}exp{-(t^(1/r)-μ)^2/(2σ^2)}
all t in R-{0}.
當 h(x) 是偶函數時, 仍是分 x≧0 與 x<0 兩部分轉換、加總.
(3)
E[1/X] = ∫_(-∞,∞) (1/x) f(x) dx = ∫_(-∞,∞) t g(t) dt.
從 X 之 p.d.f. f(x) 來看, 在 0 附近 f(x) 接近正數
{1/[√(2π)σ]} e^{-μ^2/(2σ^2)}
但 1/x 在 0 附近的積分發散, 或說不存在.
當然也可先得到 1/X = X^(-1) 的 p.d.f.
g(t) = {1/[√(2π)σt^2] e^{-(1/t -μ)^2/(2σ^2)}
而後討論 ∫_(-∞,∞) t g(t) dt 的歈散問題. 對此積分式, 我們
可以發現當 t → ±∞ 時, (g(t) 的) 指數部分趨近於一個正常數
e^{-μ^2/(2σ^2)}, 但 ∫_[M,+∞) t/t^2 dt 是發散的, 因此
∫_(-∞,∞) t g(t) dt 是發散的. 也就是說 E[1/X] 不存在.
(4)
Delta method 是一種近似法, 通常是用在所謂 "大樣本理論".
簡言之, 依中央極限定理知在適當條件下
Y_n = (X_1+...+X_n)/n asym. d. as N(μ,σ^2/n)
for large n.
考慮一可徹分一對一變換: W_n = h(Y_n), 則
W_n asymptotically distributed as N(h(μ), (h'(μ))^2 σ^2/n)
在省略掉樣本數 n 及 "大樣本" 論述, 可以說:
在 (X 的機率分配) N(μ,σ^2) 中,
若 σ << |μ|, h(x) 是一可微分一對一變換,
則可以用 N(h(μ),(h'(μ))^2 σ^2) 近似
h(X) 的機率分配.
這是因為: 在 μ 鄰近,
h(X) ≒ h(μ) + h'(μ)(X-μ)
依此,
X^r ≒ μ^r + nμ^(r-1) (X-μ)
所以
E[X^r] ≒ μ^r = (E[X])^r
但由
E[X^r] = E[(X^(r/2))^2] > (E[X^(r/2)])^2,
反過來說
(E[X^r])^2 < E[X^(2r)]
可知實際上除了 r = 1 以外, E[X^r] ≠ (E[X])^r.
回顧 delta method 的條件, 最基本的需要 σ << μ. 這個條件是
很模糊的. 事實上, 若 h 可二次微分, 由二階 Taylor's expansion,
h(X) = h(μ) + h'(μ)(X-μ) + 0.5 h"(ζ(X))(X-μ)^2
for some ζ(X) between μ and X.
此處 ζ(X) 的 "for some" 是指 "存在" 這麼一點介於 μ 和 X 之
間, 但並不確知在哪裡. 所以 delta method 採用一階近似, 除非
h"(ζ)(X-μ)^2 保證夠小, 否則計算結果可能沒什麼參考價值. 以
h(x) = x^r 來說, 就是
r(r-1)(ζ(X))^(r-2)(X-μ)^2
保證夠小. 那麼, 我們需要的就是
r(r-1)X^(r-2)(X-μ)^2 , r(r-1)μ^(r-2)(X-μ)^2
都保證夠小. 理論上是不可能的, 因為常態分配在整個 R 上都未消失,
我們只能期望止列二階項 "幾乎" 都很小.
如果對 Taylor's 二階展式的 (X-μ)^2 項取期望值, 則易知至少
r(r-1)μ^(r-2)σ^2 = r(r-1)μ^r(σ/μ)^2
要夠小. 這也是前面說至少 σ << μ 的理由之一.
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 1.165.122.98 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1599536022.A.121.html