2011年12月23日金曜日

survreg と weibreg

Rです。


survreg
> library(survival)
> library(MASS)
#library(MASS)はデータ"gehan"を使うため

> survreg(Surv(time,cens)~treat,dist="weibull", data=gehan) #デフォルトが"weibull"なので dist="weibull" は省略可
Call:
survreg(formula = Surv(time, cens) ~ treat, data = gehan, dist = "weibull")

Coefficients:
 (Intercept) treatcontrol
    3.515687    -1.267335

Scale= 0.7321944

Loglik(model)= -106.6   Loglik(intercept only)= -116.4
 Chisq= 19.65 on 1 degrees of freedom, p= 9.3e-06
n= 42

グラフはこう。
折れ線は survfit による生存関数のカプランマイヤー推定量(ノンパラメトリック)。なめらかな曲線がsurvreg による生存関数のパラメトリックモデル(ワイブル分布を仮定)。

グラフのコード:
> sr <-survreg(Surv(time,cens)~treat,dist="weibull", data=gehan) 
> sf <-survfit(Surv(time,cens)~treat, data=gehan)
> plot(sf)
> curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1])),add=TRUE)
> curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1]+coef(sr)[2])),add=TRUE)
survreg の scale = 1/(pweibull の shape)
survreg の Coefficients = log(pweibull の scale)
である。survreg はなぜか、ふつう shape と呼ぶパラメータが Scale という名前になっている。

weibreg
ehaパッケージをダウンロードしてインストールすると、weibreg という関数が使える。
> library(eha)
> weibreg(Surv(time,cens)~treat, data=gehan)
Call:
weibreg(formula = Surv(time, cens) ~ treat, data = gehan)

Covariate           Mean       Coef Exp(Coef)  se(Coef)    Wald p
treat
            6-MP    0.664     0         1           (reference)
         control    0.336     1.731     5.646     0.413     0.000

log(scale)                    3.516    33.639     0.252     0.000
log(shape)                    0.312     1.366     0.147     0.034

Events                    30
Total time at risk           541
Max. log. likelihood      -106.58
LR test statistic         19.6
Degrees of freedom        1
Overall p-value           9.29141e-06

グラフ。
weibreg による上と同じもの

グラフのコード:
> wr <- weibreg(Surv(time,cens)~treat, data=gehan)
> plot(sf)
> curve(1-pweibull(x,shape=exp(coef(wr)[3]), scale=exp(coef(wr)[2])),add=TRUE)
> curve(1-pweibull(x,shape=exp(coef(wr)[3]), scale=exp(coef(wr)[2]+-coef(wr)[1]/exp(coef(wr)[3]))),add=TRUE)
両者は、パラメータの置き方がちがう。

survreg の係数 = −(weibreg の係数)×(weibreg の shape)
survreg の Scale = 1/(weibreg の shape)

である。ぼくはかなり混乱した。

参考リンク
survreg vs. aftreg (eha)

R: weibreg
R: survreg

eha(イベントヒストリー分析)パッケージ中のオブジェクト一覧
survival(生存分析、ペナルティ付き尤度を含む)パッケージ中のオブジェクト一覧

2011年12月18日日曜日

Twitterウィルス(スパム)について

追記 2012/05/30

本家ヘルプセンターにちゃんとした記事が書かれていた↓



ぼくも前に一度ひっかかったのでメモっときます。


注意
ツイッターでこんな感じのダイレクトメッセージが来て、

誰かがブログでお前のこと悪く書いてるぞ、みたいな感じ

URLをクリックするとこんな感じにログイン画面に行った


としたら、それはスパムです。よくみるとURLが違う。twitter.comじゃなくて、twvitterとかになってる。

ここにIDとパスワードを入れると、今度は自分がスパムダイレクトメッセージをばらまく側に回ってしまいます。


対処
もしまちがってIDとパスワードを入れてしまったら、ツイッターの「設定」から、「アプリ連携」を選んで、許可したおぼえのないヤツとか、よく分からないヤツを全部消してけば(たぶん)大丈夫です。


おれのツイッターのアカウント名が写ってるが…

ついでにパスワードも変えとくといいかも。

スパムメッセージの内容は他にも「あのときの写真ができたよー、みてみて」みたいなのとか色々あります。

R: zero-inflated(ゼロ過剰)モデル

ほとんどが0でときどきピョコっと数字が入っているようなカウントデータ(なにかの数を数えたデータ)に対して回帰分析をするときは、zero-inflated(ゼロ過剰)モデルというやつを使うといいらしい。
zero-inflated モデルは、R で pscl パッケージをインストールしてzeloinflという関数を使うとできる。

例)
zeroinfl(y ~ x, dist = "poisson")


分布は "poisson"(ポアソン分布), "negbin"(負の二項分布), "geometric"(幾何分布)が使えるようである。

R: glmのオプション

R の glm 関数で


イテレーション(繰り返しの計算)回数を増やすには
glm() の引数 maxit に適当な数字をいれる。「アルゴリズムは収束しませんでした」みたいなメッセージが出ても、maxit を増やせばとりあえず収束はする(場合もある)。


イテレーションの収束判定(許容差)を変えるには
glm() の引数 epsilon に適当な数字をいれる。これを変えたいと思ったことはあんまりない。


例)
>  x <- c(18,17,15,20,10,20,25,13,12)
>  y <- c(0, 0, 0, 0, 1, 1, 0, 1, 0)
> glm(y ~ x, family=binomial(logit), maxit=1000, epsilon=0.0001) #この計算においてイテレーション回数とか変える必要は別にないが、いい例が思いつかなかった。
 
Call:  glm(formula = y ~ x, family = binomial(logit), maxit = 1000, 
    epsilon = 1e-04)

Coefficients:
(Intercept)            x  
     2.6828      -0.2111  

Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
Null Deviance:     11.46 
Residual Deviance: 10.09  AIC: 14.09 


また trace = TRUE にすると各イテレーションごとに出力してくる。
例)

> glm(y ~ x, family=binomial(logit), maxit=1000, epsilon=0.0001, trace=TRUE)
Deviance = 10.09889 Iterations - 1
Deviance = 10.09096 Iterations - 2
Deviance = 10.09096 Iterations - 3

Call:  glm(formula = y ~ x, family = binomial(logit), maxit = 1000,
    epsilon = 1e-04, trace = TRUE)

Coefficients:
(Intercept)            x
     2.6828      -0.2111

Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
Null Deviance:    11.46
Residual Deviance: 10.09 AIC: 14.09


参考URL: http://rss.acs.unt.edu/Rdoc/library/stats/html/glm.control.html

2011年12月3日土曜日

R:& と &&の違い

最近ようやくちょっと R に慣れてきた感がある

& と &&の違い
どちらも論理積なんだけど、

& だとベクトルの要素それぞれに対して
&& だとベクトルの最初の要素に対して

の演算になるみたいです。

例)
> c(TRUE,FALSE,FALSE) & TRUE
[1]  TRUE FALSE FALSE
> c(TRUE,FALSE,FALSE) && TRUE
[1] TRUE
> c(FALSE,TRUE,FALSE) & TRUE
[1] FALSE  TRUE FALSE
> c(FALSE,TRUE,FALSE) && TRUE
[1] FALSE

&& を使わなきゃいけない場面って少ない気がするので、とりあえず & をつかっておけばオッケーだと思います。| と || についても同様。

その他一般原則
  • R はとにかくベクトル‐行列が基本
  • なので for 文はなるべく使わないほうが計算が速くなる
  • このことはちょうど、数学でベクトル‐行列の計算をするとき、成分をいちいち書き下すよりもベクトル‐行列表記のままの方がかんたんだったりすることと対応している
  • for 文を使わないで代わりに何を使うのかというと apply とか merge とか subset とかである
  • ベクトルについての操作は大抵専用の速くてかんたんな関数があるので、そっちを使う
  • あと、一度書いたコードを後から見なおすのはけっこう大変なので最初からていねいに書くといい

2011年12月2日金曜日

養老孟司『バカの壁』への批判

動機
なんで今さらこんなエントリを書こうと思ったかというと、養老先生の「y=ax」を褒めてる人を最近になって読書メーターで見たからです。以下に養老先生の「y=ax」がなんでダメなのかを書きます。

y=ax
養老先生の「y=ax」というのはこんな感じです。

  • 人間の脳みそは式「y=ax」で説明できる
  • x は脳への入力(五感とか何らかの情報)、y は出力
  • a(の絶対値)は関心が大きい入力に対しては大きい値、関心が小さいことには小さい値をとる
まず、この「y=ax」っていう書き方が良くないです。入力に対していろんな値をとるなら a は定数じゃないことになります。 a を関数にして、y=a(x) とかって書くべきです。まあそれはいいや。この説明でなにが分かるかを見ていくことにします。

入力:関心大 → a:大
入力:関心小 → a:小

 でした。a が大きかったり小さかったりするとどうなるかというと、

入力:関心大 → a:大 → ax:大(= y:大)
入力:関心小 → a:小 → ax:小(= y:小)

そして、出力 y が大きいってことは、その入力に対しいての反応が大きいってことを表わしているんだそうです。なるほど。つまり、

入力:関心大 → a:大 → ax:大(= y:大) → 反応:大
入力:関心小 → a:小 → ax:小(= y:小) → 反応:小

人間の脳は、関心が大きいことには大きく反応するのか! でもこの「反応の大小」っていうのは何を表している量なんでしょう。どうも「反応の大きさ」はその入力に対する関心の大きさを表しているみたいです。

入力:関心大 → a:大 → ax:大(= y:大) → 反応:大 → 関心:大
入力:関心小 → a:小 → ax:小(= y:小) → 反応:小 → 関心:小

途中の y とか a とかを省略して結果をみてみます。

入力:関心大 → 関心:大
入力:関心小 → 関心:小

これですっきりしました。実に論理的。もっとも、この説明は一体なにを説明できているのかというと……たぶん何も説明できてないと思います。

(余談:数学をやってるとたまにこういうことがある。「よし、計算できそうだ。解けた! 0 = 0 ……ってこりゃ当たり前だよ。何やってんだおれ……」みたいな。)

この「y=ax」はもちろん褒めてる人だけじゃなくて貶してる人もいっぱいいたのですが、「あまりに単純化がひどい」みたいな趣旨のものが多く、私が言ってるようなことを書いてる人は少ないように見えました。「単純化がひどい」っていう批判も正しいのでしょうが、それ以前にこのモデル自身がまったく無意味だと思います。

ほかにも、
  • フロイトの無意識について述べた部分とかもおかしいです。「最近、おれの授業で寝てる学生が多い」 → 「学生は睡眠時間のことをちゃんと考えてない」 → 「都会っ子は無意識を意識してないのだ! だからフロイトは無意識を再発見しなければならかったのだ!」みたいな。フロイトの言う無意識って寝てる時間のことじゃないんですけど……。フロイトの「無意識」は意識がない状態を指す「無意識」とはまったく別もの。内田樹もなかよく対談なんかしてないで、つっこんでやればいいのに。まちがってることにまちがってると言うのも友情だと思うな、ぼくは。
  • ピカソについて述べた部分もおかしいです。「ピカソの(キュビズム時代の)絵はめちゃくちゃなようでいて、顔のパーツ部分など細部はリアルにかいてある」 → 「ピカソは見た景色を脳の中で独自に再構成できたのだ!」みたいなの。顔のパーツ部分、リアルか? まあ、ものの感じ方は人それぞれですからね……みたいな解釈もできるが、たぶん養老先生、ピカソの絵見てないだろ。それに、パーツ部分がリアルだったしても、ピカソの脳自体が特殊だと考えるのはかなり飛躍してると思います。
  • あとなんだっけ、共同体の再構築みたいな話にも不満があったような気がするが、忘れた。
  • それと、なんていうか、みんな「人それぞれいろんな意見がある」とか、「物ごとには良い面と悪い面がある」みたいな話好きだよなー。いや、それ自体はいいんだけど。それは議論の出発点であって結論じゃないよー。


おわりに
ぼくが「バカの壁」を読んだのはかなり前、中学生(いや、高校生だったっけ? 忘れた)のころだったので、このエントリは記憶をたよりに書いてます。なので、引用などは不正確です。でも大筋では間違っていないはずです。


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 推定量ができた.


参考文献

2011年10月28日金曜日

リンク:TeX 関連

以前に ↓このページへのリンクを貼ったことがあったが、
日本語 LaTeX を使うときに注意するべきこと

↓このページもおもしろい。
数学の常識・非常識—由緒正しい TEX 入力法

「滑稽である」「見苦しい」「不格好である」などの力強い断定でつづられた名文であり、わたしのようなボンクラ学生は居住まいを正さざるを得ない。

(ぼくの文章は冷笑的に見えたり、バカにしているように感じられたりするかもしれませんが表面的に読み取れる意味だけを読み取って欲しいです)

ちなみに、後者のリンクでは定理環境で英数字をイタリックにしないように \upshape を使うことが推奨されているが、この点については、前者のリンクで紹介されている方法の方がシンプルで良い。

そして、なぜか \usepackage{amsthm} を使うと

\begin{theorem}[定理の名前] 文

と書いたとき、

定理 1 (定理の名前) 文
(かっこは半角)

みたいな感じで「定理の名前」の部分が太字になってくれないので、ぼくは \usepackage{theorem}  の方を使ってる。こっちだと

定理 1(定理の名前)
(かっこは全角)

みたいにプリントされる。

関連リンク:theorem.sty: TeX パッケージ

あと、「重要だけど『当然の常識』や『暗黙のルール』扱いになっててあまり他人が教えてくれないこと」つながりで
数学の教科書や授業に出てくる文字と用語
このページもおもしろかった。大学の数学科三、四年生でも「定義」と「定理」の区別がついてないっぽい人はけっこういる。もったいない

2011年10月19日水曜日

TeX の数式用フォント

TeX で数式を書くときにときどき使うフォント一覧.


boldsymbol, pmb, mathfrak, mathb, を使うには \usepackage{amsmath,amssymb}
mathscr (花文字)には \usepackage{mathrsfs} が必要. 

今までベクトルとかを書く太字には
¥boldmath{¥mbox $x$}
を使っていましたが,明らかに
¥boldsymbol{x}
の方が楽ですな.

以下,上の画像のソース
デフォルト

$abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ$

Y\hspace{-.75em}=boldsymbol\{\}

$\boldsymbol{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=rm\{\}

$\rm{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathbf\{\}

$\mathbf{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathsf\{\}

$\mathsf{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathit\{\}

$\mathit{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathtt\{\}

$\mathtt{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathcal\{\}

$\mathcal{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=pmb\{\}

$\pmb{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathfrak\{\}

$\mathfrak{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathbb\{\}

$\mathbb{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

Y\hspace{-.75em}=mathscr\{\}

$\mathscr{abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ}$

2011年9月9日金曜日

カプラン‐マイヤー推定量の手短な導出


  • カプラン‐マイヤー推定量(Kaplan-Meier estimator)の何たるかについては他のサイトをご覧ください.
  • この投稿は,Thomas R. Fleming, David P. Harrington "Counting Processes and Survival Analysis" の p.3-5 あたりを参考にしています.
  • 以下はちゃんとした厳密な導出ではないです.
  • 上記の本には「3章と6章で計数過程とマルチンゲール表現によって,カプラン‐マイヤー推定量をきちんと表現する」みたいなことが書いてある.
$T$ は非負,連続型の確率変数で故障(=死亡)時間を表すのにつかう.分布関数は $F(t)$ ,密度関数を $f(t)=\frac{d}{dt}F(t)$ ,生存関数は $S(x)=\{T>t\}=1-F(t)$ ,ハザード関数を
\begin{align*}
\lambda (t) & = \lim _{\Delta t \to 0+0} \frac{1}{\Delta t}P\{t \le T<t+\Delta t|T \ge t\} \\
 & = -[\frac{d}{dt}S(t)]/S(t) \\
& = f(t)/S(t).
\end{align*}
とする.関数 $\Lambda (t)=\int^t_0 \lambda (u) du$ は $T$ の累積ハザード関数(cumulative hazard function)と呼ばれる.

また,$S(t)=\exp\{-\Lambda (t)\}$である.なぜなら
\begin{align*}
&\frac{d}{dt}\log\{1-F(t)\}=-\frac{f(t)}{1-F(t)}\\
&\Leftrightarrow \log S(t)=-\int^t_0\Lambda (u) du \\
&\Leftrightarrow S(t)=\exp\{-\Lambda (t)\}
\end{align*}
だからである.

$0=t_0<t_1<\cdots < t_m =t$ は区間 $[0,t]$ の分割, $d_l$ を $[t_{l-1},t_l)$ の時点での故障したモノの数の合計, $y_l$ を時間 $t_{l-1}$ 以前に故障していないモノ(これから故障が起きるかもしれない集団なので『リスクセット』などと呼ばれる)の数の合計とする.

十分小さい $\Delta t$ に対して,
\begin{align*}
\Lambda (t + \Delta t)-\Lambda (t) & \approx \lambda (t) \Delta t \\
& = \lim _{\Delta t \to 0+0} \frac{1}{\Delta t}P(t \le T<t+\Delta t|T \ge t) \Delta t \\
& \approx P(t \le T<t+\Delta t|T \ge t)
\end{align*}
上式の1行目

が成り立つ.これより, $\Lambda (t_l) -\Lambda (t_{l-1}) \approx P(t_{l-1} \le T \le t_l |T \ge t_{l-1})$ の推定量は,愚直に考えると $d_l/y_l$ であり,
\[\hat \Lambda (t)=\sum_{l:t_l \le t}d_l/y_l .\]
この推定量は Nelson (1986) が最初に提案したので,ネルソン(Nelson)推定量,ネルソン‐アーラン(Nelson-Aalen)推定量などと呼ばれる.

さて, $S(t) = \exp\{-\Lambda(t)\}$ であったので $\hat S(t) = \exp\{- \hat \Lambda(t)\}$ .これより,$d_l/y_l \approx 0$ のとき,
\begin{align*}
 \hat S(t) & =\prod_{l:t_l \le t} \exp(- d_l/y_l ) \\
& \approx \prod_{l:t_l \le t}(1-d_l/ y_l)
\end{align*}
2行目は一次近似(マクローリン展開の第一項目をとること)による.最後の式がカプラン‐マイヤー推定量(Kapla-Meier (1958) product limit estimator)である.

http://aaaaushisan.blogspot.jp/2011/11/blog-post.html につづく)

2011年8月12日金曜日

カリエラ型の婚姻規則とクラインの四元群(橋爪大三郎『はじめての構造主義』への部分的批判)

なんと! ここに書かれている話題については、ぼくよりちゃんとした人がちゃんとした記事をかいてくれました!
犬Q日記
しかもぼくの群の説明はまるっきりまちがっていたので訂正しました。知らないで勝手なこと書いてすみませんでした!
(9/25)

準備(置換)

置換というものがある.たとえば,

\[
\alpha=\left(
\begin{array}{cccc}
a & b & c & d \\
b & a & d & c
\end{array}
\right)
\]

と書いてあったら,αは a b c d を b a d c と並び替えるような操作ですよ,という意味になる.この操作αのことを置換という.
置換はあみだくじで表せる.αという置換をあみだくじで表すと図1のようになる.


あみだくじを通すと a b c d が b a d c に入れ替わる. また,

\[
\beta=\left(
\begin{array}{cccc}
a & b & c & d \\
c & d & a & b
\end{array}
\right)
\]

このように書いてあったらβは a b c d を c d a b と並び替える置換だ.あみだくじで表すと図2のようになる.


置換にはかけ算(かけ算は・で表す.×はあまり使わない)が定義できる.α・αはαの操作を2回繰り返すこと,α・βはαをしてからさらにβをすることだ.あみだくじでやってみる.


αを2回繰り返すと元の a b c d に戻った.このことをα・α=Iと書いたりする.Iは何もしないという置換だ(図4).

\[
I=\left(
\begin{array}{cccc}
a & b & c & d \\
a & b & c & d
\end{array}
\right)
\]


さて,α・β,ついでにβ・αは図5のようになる.



α・βもβ・αも結果は同じ d c b a になった.α・βをあらためてγと書くことにすると,


\[
 \gamma= \left(
\begin{array}{cccc}
a & b & c & d \\
d & c & b & a
\end{array}
\right)
\]
\[\alpha \cdot \beta =\gamma=\beta \cdot \alpha\]

である.
勘のいい人は気づいているかもしれない.βもγも2回繰り返すともとのa b c d に戻る.


α,β,γの関係をまとめると図8のようになる.


または,こんな表でまとめることもできる.

・ I α β γ 
II α β γ 
αα I γ β 
ββ γ I α 
γγ β α I 

(γ・βとかγ・αはあみだくじで確かめてはいないけれども,γ・β=(α・β)・β=α・(β・β)=α・I =αという具合に計算してやればこの表が完成する.)
この表みたいな構造を持ったものを『クラインの四元群』という(らしい).


本題


さて(やっと),本題に入る.
オーストラリア大陸の北部砂漠にもともと住んでいた人たちには,結婚に関して複雑な規則があった.
その内の一つ『カリエラ型』というのはこんなふうだ.
「この社会の人びとは、運動会かなにかのときみたいに、全員が四つの、だいたい人数の等しいグループ(婚姻クラス)に分かれている。(中略)父親がA1で母親がB2なら、あなたはB1。(中略)そしてあなたは、A2の異性としか、結婚できない。」(橋爪大三郎『はじめての構造主義』(講談社現代新書) p.81より)
図で書くと以下のようになる.何回か図をなぞってみると意味がわかります.



例)父親がA2で母親がB1なら? あなたはB2.結婚できるのは,A1の異性だけ.
例2)父親がB2で母親がA1なら? あなたはA2.結婚できるのは,B1の異性だけ.

図9は図8と同じ形をしている.同じ構造を持っていると言ってもいい.『カリエラ型』の婚姻規則は『クラインの四元群』と同型である,という言いかたもできるだろう.
しかし,
要するに、オーストラリアの原住民の結婚のルールは、抽象代数学の、群の構造とまったく同じものなのだ!

これは、なかなかのことではないだろうか。
ヨーロッパ世界が、えっちらおっちら数学をやって、「クラインの四元群」にたどりつくまでに、短くみても二千年はかかった。つい最近まで、だれもそんなもの、知らなかったのである。(中略)先端的な現代数学の成果とみえたものが、なんのことはない、「未開」と見下していた人びとの思考に、先回りされていたのだ。
(同 p.181より)
と言われると,うーん,そうかなー?
自分で長々と説明しておいて言うのもなんだけど,この婚姻規則を理解するのに,『クラインの四元群』や置換の話を理解する必要はまったくない.図9だけをいきなり見れば十分だ.
そして「代数学の群」というのは群の公理を満たすものを言う.
先の「置換」と「置換のかけ算」の世界はこれらが全部成り立っていた(単位元はI,αの逆元はαそのもの)から,群だけれども,ほかにも群はたくさんある.「分数」と「分数のかけ算」の世界も群だ.
「ヨーロッパ世界の数学」の成果は,「群という考え方によって,置換などの操作も,数字を扱うのと同じように数学で扱えるようにしたこと」だと思う.『クラインの四元群』は群のなかでも基本的なものだからわざわざ名前が付けられているだけで,ものを入れ替えたり対応させたり……っていうこと自体はヨーロッパでもどこでも昔からやっていた.日本でもあみだくじは昔からあった.数学の成果は置換そのもののことではない.

橋爪先生はぜんぜん水準のちがうものどうしを比較していると思う.
「二足歩行ロボットを作れるようになったのは最近だ.しかし,オーストラリア人は昔から二足歩行している.先端的な科学の成果が,「未開」と見下していた人びとに,先回りされていた.」
と言っているのに近いと思う.

この規則を発見した人はすごいと思うし,この規則が『クラインの四元群』だという指摘もおもしろいとは思う.しかし,この規則を説明するのに『クラインの四元群』を持ち出す必要があるかというと……ない.
もっと言うと,なんでわざわざ『クラインの四元群』なんて難しい言葉を使う必要があるんだ! 『知の欺瞞』じゃないか!


しかし,一方で,

ぼくは「群」とかそういう言葉を使いたくなる気持ちも分かる気がする.
神話とか婚姻規則とかを,なんというか,そのー,博物学的に羅列するだけじゃなくて,もっとこう,神話とか婚姻規則とかの構成要素の空間みたいなものを定義して,規則や神話の実際の表れ方と,その変化の仕方みたいなものを整理していったら,人間の思考パターンの本質みたいなものが,群だか環だか体だか分からないけど,そういう集合として表現できそうな気がするじゃないか! それができたらすごくおもしろいじゃないか!
たぶん,「構造主義」が目指していたのはそういうことだったんじゃないだろうか.
けれども……(参考リンク

参考文献:橋爪大三郎『はじめての構造主義』(講談社現代新書)
合わせておすすめ:内田樹『寝ながら学べる構造主義』(文春新書)

2011年8月4日木曜日

発表資料などのフォント

  • プログラムのコードなど:“Courier New”(等幅なのでインデントがずれない)
  • 数式中のアルファベット:“Times New Roman” の斜字
  • ギリシャ文字:“Symbol”

を使うとよい.

あまり関連はないがおもしろかったので
参考リンク:日本語 LaTeX を使うときに注意するべきこと(黒木玄さん)
(初学者がいきなりこれを見てしまったら「うわ、TeXってめんどくせー」って思いそうだ……)

2011年8月2日火曜日

R による生存時間分析(指数分布を仮定したパラメトリックモデル)

(あきらかに間違ったこと書いてたので訂正;8月2日)

参考リンク
[連載]フリーソフトによるデータ解析・マイニング第37回 Rと生存時間分析(2)
たとえば R で以下のようにするとき
> library(survival);library(MASS)
> ge.pasf <- survreg(Surv(time, cens)~treat,data=gehan,dist="exponential")
> ge.pasf
Call:
survreg(formula = Surv(time, cens) ~ treat, data = gehan, dist = "exponential")

Coefficients:
 (Intercept) treatcontrol 
    3.686098    -1.526614 

Scale fixed at 1 

Loglik(model)= -108.5   Loglik(intercept only)= -116.8
 Chisq= 16.49 on 1 degrees of freedom, p= 4.9e-05 
n= 42 
Coefficients がなにを意味しているのか,最初よく分からなかったが,生存曲線が \[S(t)= e^{-\lambda t} \] ただし \[ \lambda = 1/(\exp(3.686098 + (-1.526614)x)) \] ここで x は treat が control のとき1,その他は0をとる変数,
であるという意味のようだ.
(つまり,
“6-MP”のとき,$1/\exp(3.686098)=0.02506963$
“control”のとき,$1/\exp(3.686098 + (-1.526614))=0.1153846$)

ノンパラメトリックモデルと重ねてグラフを書いてみる.
(グラフのコマンド)
library(survival);library(MASS)
ge.pasf <- survreg(Surv(time, cens)~treat,data=gehan,dist="exponential")
plot(survfit(Surv(time,cens)~treat, data=gehan),lty=1:2)
curve(1-pexp(x,1/exp(coef(ge.pasf)[1])),lty=1,add=TRUE) 
curve(1-pexp( x,1/exp(sum(coef(ge.pasf))) ),lty=2,add=TRUE)
確かめるために,自分でパラメータλを推定してみる.
ある一つの生存時間の観測値が得られる確率(尤度)は \[(f(t)dt)^{d_i} (1-F(t))^{1-d_i}\] $t$は生存時間(本当は$t_i$だがしばらく$i$を省略), $d_i$は打ち切りのとき 0,打ち切りでないとき 1.
データの件数を$n$とすると,尤度関数は \[L=\prod^{n}_{i=1}(f(t)dt)^{d_i} (1-F(t))^{1-d_i}\] \[\log L = \sum^{n}_{i=1}(d_i \log(f(t)dt) + (1-d_i)\log(1-F(t)) \] いま,指数分布を仮定しているので \[F(t) = 1-e^{-\lambda t} \] \[f(t) =\lambda e^{-\lambda t} \] を代入して, \[ \frac{\partial}{\partial \lambda }\log L =0 \] を解くと \[ \hat \lambda= \left(\sum^{n}_{i=1} d_i \right)/ \left(\sum^{n}_{i=1} t_i\right)\] これをR のコードで書くと
>  sum(gehan[gehan$treat=="6-MP",3])/sum(gehan[gehan$treat=="6-MP",2]) 
[1] 0.02506964
> sum(gehan[gehan$treat=="control",3])/sum(gehan[gehan$treat=="control",2]) 
[1] 0.1153846
上と一致していることがわかる.
(補1) R の関数が \[ \lambda = 1/\exp(\mathbf{\beta x}) \] と指数をとっているのは,制約$\beta > 0$をなくして数値計算をしやすくするためと思われます.
(補2) グラフのコマンドはこれでも良い
lower.tail=FALSE
とすると上側確率を求めてくれる.
ge.pasf <- survreg(Surv(time, cens)~treat,data=gehan,dist="exponential")
plot(survfit(Surv(time,cens)~treat, data=gehan),lty=1:2)
curve(pexp(x, 1/exp(coef(ge.pasf)[1]),lower.tail=FALSE),from=0, to=max(gehan$time),lty=1,add=TRUE) 
curve(pexp(x, 1/exp(sum(coef(ge.pasf))),lower.tail=FALSE),from=0, to=max(gehan$time),lty=2,add=TRUE)

2011年7月17日日曜日

一見かんたんに計算できそうだが新しい関数になってしまう積分

これらは一見知っている関数(初等関数)になりそうにみえて新しい関数になってしまう。

通称数学科のパチンコ。

ちょっと工夫して(置換積分や部分積分で)計算しようとするのは時間の無駄なのでやめよう!


名前
\[ \int \frac{\sin x}{x} dx \] ディリクレ積分
\[ \int \frac{\cos x}{x} dx \]
\[ \int \frac{e^{-x}}{x} dx \]
\[ \int e^{-x^2} dx \] ガウス積分
\[ \int \sin (x^2) dx \] フレネル積分
\[ \int \cos (x^2) dx \] フレネル積分
\[ \int \frac{dx}{\log x} dx \] 対数積分
\[ \int \frac{dx}{\sqrt{1-x^3}} dx \] \[\int \frac{dx}{\sqrt{1-x^4}} dx \]: 楕円積分

参考文献:難波誠『数学シリーズ 微分積分学』(裳華房)

2011年7月7日木曜日

積分をかんたんに計算する工夫

前置き

数学や統計学、その他を勉強しているとこの手の計算はたまにでてくる。そこでつまづいているのはもったいないので、これくらいは暗記しておくのも悪くない。

公式

\[
\int e^x f(x) dx = e^x (f(x)-f'(x)+f''(x)-f'''(x)+ \cdots)
\]
\[
\int e^{-x} f(x) dx = -e^{-x} (f(x)+f'(x)+f''(x)+f'''(x)+ \cdots)
\]
\[
\int f(x) \sin(x) dx = f(x)\sin(x) + f'(x)\cos(x) +f''(x)(- \sin(x)) + f'''(x)(-\cos(x))+ \cdots
\]
\[
\int f(x) \cos(x) dx = f(x)(-\cos(x)) + f'(x)(\sin(x)) +f''(x)\cos(x) + f'''(x)(-\sin(x))+ \cdots
\]

証明
部分積分の公式を繰り返し適用する。



\begin{align*}
\int xe^{-x} dx & = -e^{-x}(x + x') \\
& = -(x + 1)e^{-x} .
\end{align*}
\begin{align*}
\int x^2 e^{-x} dx & = -e^{-x}(x^2 + (x^2)'+(x^2)'') \\
& = -(x^2 + 2x + 2)e^{-x} .
\end{align*}

参考
藤田岳彦「弱点克服 大学生の確率・統計」(東京図書)

余談
数学科には暗記を嫌う人が多い気がする。その理由はだいたい以下の2通りに分類できる。

  1. 理由はわからないけどとにかくいやだ。嫌いなものは嫌い。どうしても暗記できない。
  2. 数学は暗記よりも考える力のほうが大事だから
前者はその人が背負ってしまった業なのでしようがない。生きていくしかない。

後者はたぶん間違ってる。どこが間違ってるかというと暗記に頼らないほうが考える力がつくっていうようなことはない。

考えるっていうのは脳みそのワーキングメモリの上でオブジェクト(言葉やモノ)をなにかの規則通りに操作することであるので、記憶したり、紙に書いたりして、ワーキングメモリをいっぱいあけておいたほうがいろいろ考えられる。

参考リンク:あらゆる勉強に通じるコツ

2011年6月7日火曜日

Windows XP に LaTeX beamer をいれる

Beamer をいれる方法.

↓このページをみれば分かります.
よしいずの雑記帳
けれども,これだけだとTeXファイルをコンパイルするときによくわからないエラーが出てしまう.

そこで,↓このページに書いてあるように
Triad sou.
ここ
から beamerbasecompatibility.sty をダウンロードしてきて,今入っている beamerbasecompatibility.sty を探して上書きする.

文字を14ポイント以上にあげるときはさらにパッチが必要で,
LaTeX Beamer 入門
にあるように,
ここ
から extsizes.zip をダウンロード&解凍して,texmf フォルダに移す(一番最初にやったのと同じ).

これでセットアップは終了.

その他メモ

上記リンクに beamer で作ったスライドのお手本がある.

デフォルトだと,本文だけでなく数式もゴシックになってしまうので
\usefonttheme[onlymath]{serif}
を使うと良い.

命題や主張を枠で囲うのには,
¥begin{block}{ここにブロックのタイトル}

¥end{block}
が便利.

(S.M. 君に感謝)

2011年3月26日土曜日