2011年11月30日水曜日

R: ヒストグラムに密度関数のグラフを重ねてかく

R でヒストグラムの上に密度関数を重ねた図がかきたいときは、hist() の引数 freq に FALSE を与えてやればいい。
例)
> r =rnorm(100)
> hist(r,freq=FALSE)
> curve(dnorm(x),add=TRUE,lty=2) 

ヒストグラムの縦軸が density(密度)というのは分かりにくいかもしれないが、
面積=たて×よこ=割合
である。
この例の場合、ヒストグラムの幅は0.5刻みなので、たとえば、0.5〜1の範囲に入る標本の個数は
(約0.3) × 0.5 = 約0.15
全体の約15%=15個くらい、と分かる。

参考リンク:Rと統計学―その4

2011年11月11日金曜日

リンク:GLM

「指数型分布族のナチュラルパラメータは一般化線形モデル(R の glm)の自然リンク関数である」
という話をメモっておこうかと思ったのですが,このページが分かりやすかったのでリンクを貼っておくことにしました(そしてこのページはおそらく下記の本を参照していると思われます).


これはたとえば,二項分布の確率関数が
\begin{align}
f(y;n,p) &= \left( \begin{array}{l}
n \\
y \\
\end{array}\right) p^n(1-p)^{n-y} \\
&= \exp( \log( \left( \begin{array}{l}
n \\
y \\
\end{array}\right) p^n(1-p)^{n-y} ))\\
&= \exp (y \log \frac{p}{1-p} + n \log(1-p) + \log\left( \begin{array}{l}
n \\
y \\
\end{array}\right) )
\end{align}
と表わせ,これの $y$ とかかっている部分( $\log \frac{p}{1-p}$ )がロジットになっているということが,二値の説明変数のモデリングにロジスティック回帰がよく使われることの,一つの理論的根拠になっているということのようです.なんか不思議な感じがする…….

ネルソン推定量からカプラン・マイヤー推定量を導出する

本エントリは カプラン‐マイヤー推定量の手短な導出 の続き的な内容です.


準備

確率変数 $\overline{N}(t)$ は時刻 $t$ までに故障していないモノ(これから故障が起きるかもしれない集団『リスクセット』)の数の合計, $\overline{Y}(t)$ は時刻 $t$ までに起こったイベント数(故障したモノの数の合計)にそれぞれ対応していると考えることにする.
また,$f(t) - f(t-) $ を表すのに $\Delta f(t)$ と書くする.


導出

累積ハザードを
\begin{align}
\Lambda (t) &= \int^{t}_{0}\{1-F(s-)\}^{-1}dF(s),\\
(&= \int^{t}_{0}\frac{f(s)}{1-F(s-)}\, ds )
\end{align}
と定義する.これより,
\begin{eqnarray*}
d\Lambda(s) &=& \frac{dF(s)}{1-F(s-)},\\
dF(s) &=& \{1-F(s-)\}d\Lambda(s).
\end{eqnarray*}
$F(0)=0$ より,
\begin{eqnarray}
F(t) &=& \int^{t}_{0}dF(s) \nonumber\\
&=& \int^{t}_{0} \{ 1- F(s-) \}d\Lambda (s).
\tag{1}
\end{eqnarray}


$\Lambda$ の推定量,ネルソン (Nelson) 推定量 $\hat{\Lambda} = \int \overline{Y}^{-1} d\overline{N}$ は使っていいことにする.


$S$ の推定のために,$\Lambda$ のネルソン推定量 $\hat{\Lambda} = \int \overline{Y}^{-1} d\overline{N}$ を式(1)に代入する.
つまり,式(1)は
\begin{align*}
1-F(t) &= 1-\int^{t}_{0}\{1-F(s-)\}d\Lambda(s), \\
S(t) &=1- \int^{t}_{0}S(s-)d\Lambda(s),
\end{align*}
を含意するので,これに $\hat{\Lambda}$ を代入する.
$\hat{S}$ は再帰的に定義され,
\begin{eqnarray}
\hat{S}(t) = 1- \int^{t}_{0} \hat{S}(s-)d\hat{\Lambda}(s).
\end{eqnarray}
これより
\begin{align*}
\Delta \hat{S}(t) &= \hat{S}(t)-\hat{S}(t-), \\
\hat{S}(t-)-\hat{S}(t) &= -\Delta \hat{S}(t) \\
&= \hat{S}(t-) \Delta \hat{\Lambda}(t) \\
&= \hat{S}(t-) \Delta \left( \int \overline{Y}^{-1}(t) d\overline{N}(t) \right) \\
&= \hat{S}(t-) \frac{\Delta \overline{N}(t)}{\overline{Y}(t)} ,\\
\hat{S}(t) &= \hat{S}(t-)- \hat{S}(t-)\frac{\Delta \overline{N}(t)}{\overline{Y}(t)}\\
 &= \hat{S}(t-) \left\{ 1- \frac{\Delta \overline{N}(t)}{\overline{Y}(t)} \right\},
\end{align*}
であるから,
\begin{align*}
\hat{S}(t-) &= \hat{S}((t-)-)\left\{1-\frac{\Delta \bar{N}(t-)}{\overline{Y}(t-)}\right\}, \\
\hat{S}((t-)-) &= \hat{S}(((t-)-)-)\left\{1-\frac{\Delta \bar{N}((t-)-)}{\overline{Y}((t-)-)}\right\}, \\
\vdots
\end{align*}
よって
\begin{align*}
\hat{S}(t) = \prod _{s \le t} \left\{ 1- \frac{\Delta \overline{N}(s)}{\overline{Y}(s)} \right\},
\end{align*}
これで以前に紹介した Kaplan-Meier 推定量ができた.


参考文献