2012年11月18日日曜日

僕の新しいブログ

廿TT

2012年11月11日日曜日

R による大量データの散布図:あるいは iplots for mac

要約:R で大量データの散布図を描くとき、普通に plot を使うと点が重なりあってなんだかよくわからないことになるかもしれない。そういう場合は iplots パッケージとか使って透明度をあげたりすればいい。

 iplots で
とりあえず winows ユーザーの方はRコンソールをたちあげて

x <- rnorm(10000)
y <- rnorm(10000)
install.packages("iplots")
library(iplots)
iplot(x, y)
と入力してみてください。 矢印キーで操作できます。

mac で library(iplots) とやると、

> library(iplots)
 要求されたパッケージ rJava をロード中です
Note: On Mac OS X we strongly recommend using iplots from within JGR.
Proceed at your own risk as iplots cannot resolve potential ev.loop deadlocks.
'Yes' is assumed for all dialogs as they cannot be shown without a deadlock,
also ievent.wait() is disabled.
Error :  .onLoadはloadNamespace()(iplotsに対する)の中で失敗しました、詳細は:
 call: .jnew("org/rosuda/iplots/Framework")
 error: java.lang.InternalError: Can't start the AWT because Java was started on the first
thread.  Make sure StartOnFirstThread is not specified in your application's Info.plist or on the command line
 エラー:  '‘iplots’' に対するパッケージもしくは名前空間のロードが失敗しました
というエラーが出る。しかし JGR をインストールすれば mac でも動く。

JGRのコンソール

その他の方法
hexbin パッケージではこんなこともできる

install.packages("hexbin")
library(hexbin)
x <- rnorm(10000)
y <- rnorm(10000)
bin<-hexbin(x, y, xbins=50)
plot(bin, main="Hexagonal Binning")
「それほど大量のデータではないが重なりがある」みたいな場合は、jitter 関数でノイズを加えて点の座標をずらしてやればいい。
x <- rep(1:3,3)
y <- rep(1:3,3)
par(mfrow=c(1,2))
plot(x, y)
plot(jitter(x,0.1), jitter(y,0.1))

リンク集

Y先生、並びにK先生に感謝いたします。ありがとうございました。
https://twitter.com/yoshiroy/status/250231399150395392

2012年10月11日木曜日

クリップボード経由で R とExcel を使う

追記(11/17):せっかくなので関数化した。
マトリックス or データフレームならそのまま、それ以外は転置をとってクリップボードに書き出す。

#mac の場合
cpy <- function(a){
    if(is.data.frame(a) | is.matrix(a)){ write.table(a, pipe("pbcopy"), sep="\t", row.names=FALSE, quote=FALSE)
    }else write.table(t(a), pipe("pbcopy"), sep="\t", row.names=FALSE, quote=FALSE)
}
#windows の場合
cpy <- function(a){
    if(is.data.frame(a) | is.matrix(a)){ write.table(a, file="clipboard", sep="\t", row.names=FALSE, quote=FALSE)
    }else write.table(t(a), file="clipboard", sep="\t", row.names=FALSE, quote=FALSE)
}

ちなみに .Rprofile というファイル名で作業ディレクトリに保存しておくと、R を起動したとき自動的に読み込まれる。mac の場合は .Rprofile の改行コードを LF にする必要があるようだ。



ちょっとした表を作るときなどに使う小技をメモ


ウィンドウズの場合:
クリップボードから読み込み;
a <- read.table("clipboard", header=TRUE) 
クリップボードに書き出し;
write.table(a, file="clipboard",sep="\t" ,row.names=FALSE)
# "\t" はタブ

マックの場合:
クリップボードから読み込み;
a <- read.table(pipe("pbpaste"))
クリップボードに書き出し;
write.table(a, pipe("pbcopy"), sep="\t", row.names=FALSE)


例えばこんな Excel の表があるとして、

コピーして
read.table で R に読み込むと

こんな感じ。

あるいはwrite.tableで R から書きだして

ペーストして、横に並べたりして
行名を書き出すのは
write.table(row.names(a), pipe("pbcopy"), sep="\t",row.names=FALSE, col.names=FALSE)
(ウィンドウズだと write.table(row.names(a), file="clipboard", sep="\t", row.names=FALSE, col.names=FALSE) )
で、

 こんな感じだ。


参考&関連リンク:
R の初歩
ちょっとした R のこと
表計算シートをクリップボード経由で読み込みクロス表(分割表)を作成・保存する手順

2012年8月10日金曜日

円グラフやめろ


要約:円グラフはやめて、代わりにドットプロットや棒グラフを使いましょう。

 いままで見落としていたが、統計ソフト R のヘルプに "Pie charts are a very bad" なんて記述があっておどろいた。円グラフが悪いってどういうこと? 以下拙訳。
注意
 円グラフ(パイチャート)は情報を提示する方法として、とても悪いもののひとつです。(人間の)目は長さを判別するのは得意ですが、相対的な面積を判別するのは苦手です。この手のデータを図示するのには、棒グラフ(柱状グラフ)やドットプロットが使えます。
 Cleveland (1985)、264 ページには『円グラフで表せるようなデータは、常にドットプロットで表せる。この意味は、(人間が不得意な)正確な角度の判別の代わりに、(人間が得意とする)共通の尺度の下での位置の判別で行えるということである』とあります。この立場は、Cleveland の実証的研究および、McGill の知覚心理学者としての調査に基づいています。 
参考文献
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.
Cleveland, W. S. (1985) The Elements of Graphing Data. Wadsworth: Monterey, CA, USA.

ということだそうです。はい。わかりました。これからは円グラフをやめて、ドットプロットや棒グラフを使います。

円グラフ。0.3 は 0.2 よりだいぶ大きく見える。0.3 と 0.4 は大差ない感じがする(私見)。

棒グラフ

ドットプロット

 上図のコード(R):
X <- c(0.1,0.2, 0.3, 0.4)
pie(X, labels=X)
barplot(X, horiz=T, space=0)
dotchart(X, pch=16)

参考リンク:
 奥村先生のページは全体的に参考になる。
など。
 ウィキペディアにも円グラフについての批判が書かれている。知らなかった……。

2012年6月28日木曜日

普通の日記:Chrome no longer supports Mac OS X 10.5

追記(8/18):その後、けっきょく Firefox にのりかえた。


chrome の「よくみるページ」が突然表示されなくなった。


そしてこんなメッセージが
詳しく見る を押すと


"Chrome no longer supports Mac OS X 10.5"


"Upgrading your Mac to 10.6 (Snow Leopard) or 10.7 (Lion)"

だそうです。 はい。けどおれはとりあえずこの黄色いメッセージが消えてくれればいい。

カチカチやっても消えない


どうすればこのメッセージなくなるの……。

解決策:New Tab Redirect! という拡張機能を入れて My new tab page を about:blank にする。


これですっきりしました。ぺんぺん草いっぽん生えない快適な画面。これで明日からまたらくらくインターネット生活だ!

2012年6月19日火曜日

母数とはなにか

はじめに


「母数とは」で検索すると以下のような記事が出てきます。

母数の意味 - 数学 - 教えて!goo (現時点では google 検索で上から3番目)
分母と母数の違いを教えてください。 - livedoor ナレッジ (4 番目)

これらの「ベストアンサー」にはまちがいが書いてあります。

どうやら「『母数と分母はちがう!』という話はどこかで聞いたことがあるけど、けっきょく『母数ってなんなのか』がよくわからないので、やはり微妙に間違った知識だけが頭に残ってしまった」という人がけっこういるようです。

そうなってくると、単に 「母数」を「分母」の意味で使うのはやめろ と言っただけでは無責任なのかもしれません。

そこでちょっとぼくなりに「母数」を説明してみます。


とりあえず


  • 「母数」は「分母の数」や「全数」とは、まったく別ものと思ってください。
  • 母数はパラメーターparametor)の訳語です。「母数」というより「パラメーター」といった方がわかりやすいかもしれません。
  • 統計学の重要な目的の一つは、母数推定することです。
  • 統計学においては、母数推定されるものなのです。


母数を推定してみる


例えばコインを4回投げて、

裏、表、表、裏

というデータが得られたとします。
「表が出る確率」を p で表すことにします。さて、p はいくつでしょうか?
こう質問されたら多くの人が、「4回中2回が表になってるから、『表が出る確率』は半々、p =1/2 じゃない?」と推測すると思います。

しかし、同じコインをさらに4回投げて、

裏、表、表、裏、裏、裏、裏、裏

というデータが得られたとします。

このとき「表が出る確率」p はいくつでしょうか?
「このコインは裏が出やすいようだ。8回中2回しか表がないから、p2/8 = 1/4 だ」
と考える人もいるでしょうし、
「裏が続けてでたのは偶然だろう。いくらなんでも p1/4 ってことはないだろう
と考える人もいるでしょう。

では、正解は? 本当の p はいくつなのでしょうか?
それはわかりません。神のみぞ知る世界です。我々はデータから p を推測するしかないのです。
この「本当の p 」のことを母数というのです。

この例では母数は p だけでしたが、母数が何種類もある場合もあります。
例えば、コインの裏表でなく、ルーレットのデータだったら? あるいは複雑なアンケートのデータだったら? 母数が何種類もないと困りそうです。


統計学の世界観



統計学のわかりにくさは、データ(標本)から推定した p = 1/2 や p1/4 が、母数の p そのものではないというところにあるような気がします。
統計学では、例えばデータの平均をとっても、それが平均そのものであるとは考えません。データの平均は真の平均の推定値であると考えるのです。「真の平均」というのも場合によっては母数の一種だったりします。

統計学の世界観

プラトンは「我々が見ているこの現実は、洞窟に映る影のようのものだ。我々は真の実体(イデア)そのものではなく、イデアの影を見ているにすぎない」と考えました。

プラトンの世界観

この考え方は統計学に近いかもしれません。統計学では、影のような現実からデータを集めて真実に近づくための方法を探しているのです。

またオタク諸氏はフィギュアを眺めたり、プラモデルを組み立てたりしながら、二次元の美少女を夢想します。

オタクの世界観

この考え方も統計学に近いかもしれません。統計学では、プラモデルの代わりに「統計モデル」を組み立てたり眺めたりして、本質を理解しようとするのです。
先ほどの「表が出る確率は1/4」とかいうのも、「統計モデル」の一種です。


最後に


母数というのは、wikipedia にある通り「確率分布を特徴付ける数」のことです。なので「確率分布」のことを勉強すれば「母数」のことは自然と理解できます。逆に「確率分布」の勉強をしないで一足飛びに「母数」だけを理解しようとしてもあまり意味がありません、ながながと書いておいて最後にこんなことをいうのもアレですが。
先ほどの「表が出る確率は1/4」とかいうのも、「確率分布」の一種で「ベルヌーイ分布」という立派な名前がついています。


参考文献:
森 博嗣『数奇にして模型』(講談社ノベルス)
本田 透『喪男(モダン)の哲学史』(講談社)

2012年6月17日日曜日

R: 不完全ガンマ関数の計算

pgamma(x, a) で \[ \frac{\gamma(x, a)}{\Gamma(a)}, \] gamma(a) で Γ(a) が求められるので,

  • 第1種不完全ガンマ関数(下側) \[ \gamma(x, a) = \int_0^x t^{a-1}\,e^{-t}\,dt \] は,pgamma(x, a) * gamma(a)

  • 第2種不完全ガンマ関数(上側) \[ \Gamma(x, a) = \int_x^{\infty} t^{a-1}\,e^{-t}\,dt \] は,pgamma(x, a, lower=FALSE) * gamma(a)
とやれば計算できる.

参考リンク;
ガンマ分布(Wikipedia)
不完全ガンマ関数(Wikipedia)
R: The Gamma Distribution
不完全ガンマ関数 - 高精度計算サイト
オーバーフロー、アンダーフロー - さかもとたつお

2012年6月5日火曜日

計算メモ: 生存関数と平均の関係

\begin{align} \mu = \int^{\infty}_{0}x f(x) \, dx \end{align} とし(もう少していねいに言うと,確率密度関数が $f(x)$, support が $[0,\infty)$, 平均が $\mu < \infty$ の分布を仮定している), かつ $\displaystyle \lim _{x \to \infty} f(x) = 0$ とする.
このとき, \begin{align} \int^{\infty}_{0} [1 - F(x) ]\, dx = \mu. \end{align}

なぜなら, 部分積分により, \begin{align} \mu &= \int^{\infty}_{0}x f(x) \, dx \\ &= \int^{\infty}_{0}x \{-[1 - F(x)] \}' \, dx \\ &= \biggl[ x \{-[1 - F(x)] \}\biggr]^{\infty}_{0} - \int^{\infty}_{0} \{-[1 - F(x)] \} \, dx \\ &= \int^{\infty}_{0} [1 - F(x) ]\, dx \end{align} だから.
ここで 3 行目の第一項が 0 になるのは, ロピタルの定理を用いて, \begin{align} & \lim _{x \to \infty} \left[ x \{-[1 - F(x)] \} \right] \\ &= \lim _{x \to \infty} \left[ - \frac{x}{\{1-F(x)\}^{-1}} \right] \\ &= \lim _{x \to \infty} \left[- \frac{1}{\{-f(x)\}^{-2}} \right] \\ &= \lim _{x \to \infty} \left[ -\{f(x)\}^{2} \right] \\ &= 0 \end{align} のように計算した.

この性質はときどき当たり前のように使われている.
たとえばこの本↓ の 235 ページに出てくる.

 

2012年5月19日土曜日

『計算統計学の方法 』読書ノート2;EMアルゴリズムの練習問題(生存時間分析)

追記
グラフを入れた (6/1)



この本↓ の事例を勉強のために実装してみる. R を使用.


小西 貞則, 越智 義道, 大森 裕浩(2008)『シリーズ〈予測と発見の科学〉 5 計算統計学の方法 ―ブートストラップ,EMアルゴリズム,MCMC―』(朝倉書店)


「6.4 中途打ち切りデータと単回帰」より.
前回は一応,本を読んでなくても理解できるように文章を書いたけど,なんか帯に短したすきに長し的な感じになってしまったので,今回は計算メモとコードだけを羅列することにした

計算メモ

  • p.98 (6.26) 式 $z^{(k)}_{i} = E_{\boldsymbol{\theta}^{(k)}}[Z_i | \boldsymbol{Y}=\boldsymbol{y}]$ の計算

いま,計算したいのは $i=m+1,\ldots ,n$ の場合なので, \begin{align} z^{(k)}_{i} &= E_{\boldsymbol{\theta}^{(k)}}[Z_i | \boldsymbol{Y}=\boldsymbol{y}] \\ &= E_{\boldsymbol{\theta}^{(k)}}[X_i | X_i > c_i] \end{align} 連続型の分布の時,条件付き期待値$E[X | X > c]$は,
\begin{align} E[X | X > c] &= \frac{\int^{\infty}_{c}x \cdot f(x)dx}{P(X>c)}\\ &= \frac{\int^{\infty}_{c}x \cdot f(x)dx}{1-F(c)} \end{align} なので,p.98 の \begin{align} \int^{\infty}_{c} x \phi \left(\frac{x-\mu}{\sigma} \right) \, dx= \mu\left\{1-\Phi \left(\frac{c-\mu}{\sigma} \right) \right\} +\sigma \phi \left(\frac{c-\mu}{\sigma} \right) \end{align} より \begin{align} z^{(k)}_{i} &= \frac{\int^{\infty}_{c_i} x_i \phi \left(\frac{x_i-\mu^{(k)}_i}{\sigma^{(k)}} \right)\, dx_i}{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}} \right)} \\ &=\frac{ \mu^{(k)}_{i}\left\{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}}{\sigma^{(k)}} \right) \right\} +\sigma^{(k)} \phi \left(\frac{c_i-\mu^{(k)}_{i}}{\sigma^{(k)}} \right)}{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}} \right)} \\ &= \mu^{(k)}_{i}+\sigma^{(k)} \frac{\phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}}\right) }{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}}\right)}. \end{align}

  •  p.99-100 Mステップの計算
 \begin{align} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)})=& - \frac{n}{2} \log2\pi - \frac{n}{2} \log \sigma^2 -\frac{1}{ \sigma^2}\sum^{n}_{i=1}(\mu_i)^2 \\ &- \frac{1}{ \sigma^2} \left\{ \sum^{m}_{i=1} x^2_1 -2 \sum^{m}_{i=1} \mu_i x_i + \sum^{n}_{i=m+1} (z^2_i)^{(k)} -2 \sum^{n}_{i=m+1} \mu_i z^{(k)}_{i}\right\} \end{align} ベクトルはベクトルのまま微分する方法があったような気がするが、よく分からないので一つずつ計算する. \begin{align} \mu_i &= \boldsymbol{v}^{T}\boldsymbol{\beta}\\ &=\beta_0+\beta_1v_{1i}+\beta_2v_{2i} \end{align} なので, \begin{align} 0 &=\frac{\partial}{\partial \beta _0} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &=- \frac{1}{2 \sigma^2} \sum^{n}_{i=1}2\mu _i - \frac{1}{ 2\sigma^2}\left\{ -2 \sum^{m}_{i=1} x_i + -2 \sum^{n}_{i=m+1} z^{(k)}_{i} \right\} \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1}(\mu_i -x^{\ast}_{i})\right\}, \end{align} \begin{align} 0 &=\frac{\partial}{\partial \beta _1} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &=- \frac{1}{2 \sigma^2} \sum^{n}_{i=1}2\mu _i v_{1i}- \frac{1}{ 2\sigma^2}\left\{ -2 \sum^{m}_{i=1} v_{1i} x_i + -2 \sum^{n}_{i=m+1} v_{1i} z^{(k)}_{i} \right\} \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1}v_{1i} (\mu_i -x^{\ast}_{i})\right\}, \end{align} 同様, \begin{align} 0 &=\frac{\partial}{\partial \beta _2} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1}v_{2i} (\mu_i -x^{\ast}_{i})\right\}, \end{align} ここで, \begin{align} \boldsymbol{x}^{\ast} =(x^{\ast}_{1},\cdots , x^{\ast}_{n})^T = (c_{1}, \cdots ,c_m, z^{(k)}_{m+1}, \cdots, z^{(k)}_{n})^T . \end{align} まとめると, \begin{align} \boldsymbol{0} &=\frac{\partial}{\partial \boldsymbol{\beta}} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1} \boldsymbol{v}_i (\mu_i -x^{\ast}_{i})\right\}. \end{align} これを解いて, \begin{align} \boldsymbol{0}&= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1} \boldsymbol{v}_i (\mu_i -x^{\ast}_{i})\right\} \\ &= \sum^{n}_{i=1} \boldsymbol{v}_i (\boldsymbol{v}^T_i \boldsymbol{\beta} -x^{\ast}_{i}) \\ &= \sum^{n}_{i=1} \boldsymbol{v}_i \boldsymbol{v}^T_i \boldsymbol{\beta} - \sum^{n}_{i=1} \boldsymbol{v}_i x^{\ast}_{i} \\ \boldsymbol{\beta}^{(k+1)} &= \left( \sum^{n}_{i=1} \boldsymbol{v}_i \boldsymbol{v}^T_i \right)^{-1} \left( \sum^{n}_{i=1} \boldsymbol{v}_i x^{\ast}_{i} \right) . \end{align} その他の計算は本にある通り.

R コード
aml <-data.frame(
  d =c( 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0),
  t =c(18, 9,28,31,39,19,45, 6, 8,15,23,28, 7,12, 9, 8, 2,26,10, 4, 3, 4,18, 8, 3,14, 3,13,13,35),
  v1=c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
  v2=c( 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0))
aml <-aml[order(aml$d,decreasing =T),]
aml <-transform(aml,x=log(t))
attach(aml)
v <- t(cbind(1, v1, v2))
#初期値
beta <-c(beta_0 =2.43, beta_1=-0.59, beta_2=-0.11)
sigma2 =0.57
theta<-matrix(0, nrow=19, ncol=4)
c <-x[d==1] #打ち切りでない部分
z <-x[d==0] #中途打ち切りデータの部分
n <-length(d)
for(k in 1:19){
  mu <-t(v) %*% beta
  h= dnorm((z-mu[d==0])/sqrt(sigma2))/(1-pnorm((z-mu[d==0])/sqrt(sigma2)))
  zk=mu[d==0]+sqrt(sigma2)*h
  z2k=mu[d==0]^2+sigma2+sqrt(sigma2)*(z+mu[d==0])*h
  xast <-c(c,zk)
  beta2 <- solve(v %*% t(v)) %*% (v %*% xast)
  mu2 =t(v) %*% beta2
  sigma22 <- (sum(c^2-2*mu2[d==1]*c+mu2[d==1]^2)+sum(z2k-2*mu2[d==0]*zk+mu2[d==0]^2))/n
  mu =mu2
  beta =beta2
  sigma2 =sigma22
  theta[k,] =c(t(beta), sigma2)
}
detach(aml)
#library(xtable)
#print(xtable(theta,digits=7),type="html")

実行結果

収束の様子はこんな感じ



$k$ $\beta^{(k)}_{0}$ $\beta^{(k)}_{1}$ $\beta^{(k)}_{2}$ $(\sigma^2)^{(k)}$
1 3.1489857 -0.9346340 -0.2141950 0.6825710
2 3.2484719 -1.0066516 -0.2489479 0.7693717
3 3.2811932 -1.0289521 -0.2613609 0.8030753
4 3.2930759 -1.0369361 -0.2659463 0.8158525
5 3.2974980 -1.0398975 -0.2676593 0.8206739
6 3.2991564 -1.0410071 -0.2683023 0.8224909
7 3.2997799 -1.0414242 -0.2685441 0.8231754
8 3.3000147 -1.0415812 -0.2686352 0.8234332
9 3.3001030 -1.0416404 -0.2686695 0.8235303
10 3.3001363 -1.0416626 -0.2686824 0.8235669
11 3.3001489 -1.0416710 -0.2686872 0.8235806
12 3.3001536 -1.0416742 -0.2686891 0.8235858
13 3.3001554 -1.0416753 -0.2686898 0.8235878
14 3.3001560 -1.0416758 -0.2686900 0.8235885
15 3.3001563 -1.0416760 -0.2686901 0.8235888
16 3.3001564 -1.0416760 -0.2686902 0.8235889
17 3.3001564 -1.0416761 -0.2686902 0.8235889
18 3.3001564 -1.0416761 -0.2686902 0.8235890
19 3.3001564 -1.0416761 -0.2686902 0.8235890


答え合わせ
> sr <- survreg(Surv(x,d)~v1+v2,dist="gaussian", data=aml)
> sr
Call:
survreg(formula = Surv(x, d) ~ v1 + v2, data = aml, dist = "gaussian")

Coefficients:
(Intercept)          v1          v2 
  3.3001564  -1.0416761  -0.2686902 

Scale= 0.907518 

Loglik(model)= -36.8   Loglik(intercept only)= -40.7
 Chisq= 7.71 on 2 degrees of freedom, p= 0.021 
n= 30 
#survreg は分散でなく標準偏差で表示しているので,
> (sr$scale)^2
[1] 0.823589
> sr$it
[1] 4
survreg は 4回のイテレーションで収束している.優秀.

参考文献(孫引き)
$\int^{\infty}_{c} x \phi \left(\frac{x-\mu}{\sigma} \right) \, dx$ と $\int^{\infty}_{c} x^2 \phi \left(\frac{x-\mu}{\sigma} \right) \, dx$ の計算はこの本,解析に使用したAMLデータはこの本に出ているらしい.

2012年5月6日日曜日

『計算統計学の方法 』読書ノート;EMアルゴリズムの練習問題(入門編)

追記
グラフを入れた (6/1)


この本↓ の事例を勉強のために実装してみる. R を使用.


小西 貞則, 越智 義道, 大森 裕浩(2008)『シリーズ〈予測と発見の科学〉 5 計算統計学の方法 ―ブートストラップ,EMアルゴリズム,MCMC―』(朝倉書店)


「6.1 1変量正規分布の場合」より.


これは, 欠測データがある場合に, 正規分布のパラメータ $\mu$ と $\sigma ^2$ を推定する問題である.
以下の様なデータがあるとする.
54.62, 57.98, 49.41, 51.77, 46.72, 43.09, NA, 51.68, 40.33, 50.78, 43.63, NA, 44.89, 52.34, 47.28, 45.73, 49.02, 49.52, 46.41, 45.44, 46.50, 53.57, NA
NA の部分は何らかの理由で値が分からない. この NA の部分を欠測データと呼ぶ.
こういう場合には, EM アルゴリズムが使える.

準備
標記の都合上, 上のデータを並べかえる. そして, 前半の l 個は観測の行われたデータ, 後半の m 個が欠測データであることにする. また, 完全な観測に対応する確率変数ベクトルを $\boldsymbol{X}=(X_1, \ldots,X_{l},X_{l+1}, \ldots, X_n)^T$, 実際に観測の行われたデータに対応する確率変数ベクトルを $\boldsymbol{Y}=(Y_1, \ldots, Y_l)^T$, 欠測データに対応する確率変数ベクトルを $\boldsymbol{Z}=(Z_1,\ldots,Z_m)^T$として,その実現値を小文字で表すことにする(T は転置).
こんな風に.
$x_1$ $x_2$ $x_3$ ...... $x_{18}$ $x_{19}$ $x_{20}$ $x_{21}$ $x_{22}$ $x_{23}$
$y_1$ $y_2$ $y_3$ ...... $y_{18}$ $y_{19}$ $y_{20}$ $z_{1}$ $z_{2}$ $z_{3}$
54.62 57.98 49.41 ...... 45.44 46.50 53.57 ? ? ?

$X_i$ ($i =1, \ldots ,n$) は正規分布に従うと仮定する.
このとき, (仮に欠測部分のデータが観測されたとして)完全データに対する尤度関数 $L^C$ は, \begin{align} L^C (\boldsymbol{\theta},\boldsymbol{x}) &= \prod ^{n}_{i=1} f(x_i) \\ &= \prod ^{n}_{i=1} \frac{1}{\sqrt{2 \pi \sigma ^2}} \exp\left(- \frac{(x_i - \mu)^2}{2 \sigma ^2} \right) \end{align} 対数尤度関数 $l^C$ は \begin{align} l^C (\boldsymbol{\theta},\boldsymbol{x}) &=\log \left( \prod ^{n}_{i=1} \frac{1}{\sqrt{2 \pi \sigma ^2}} \exp\left(- \frac{(x_i - \mu)^2}{2 \sigma ^2} \right)\right) \\ &=\sum^{n}_{i=1}\left( -\frac{1}{2} \log (2 \pi \sigma^2) - \frac{(x_i - \mu)^2}{2 \sigma ^2} \right) \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\sum^{n}_{i=1} (x_i - \mu)^2 \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (x_i - \mu)^2 + \sum^{n}_{i=l+1} (x_i - \mu)^2 \right\} \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (y_i - \mu)^2 + \sum^{m}_{j=1} (z_j - \mu)^2 \right\} \end{align} である.ただし,ここで $\boldsymbol{\theta}= (\mu, \sigma^2)$ である.

E ステップ
E ステップでは,$Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(k)}) = E_{\boldsymbol{\theta}^{(k)}}[l^C (\boldsymbol{\theta},\boldsymbol{X})|\boldsymbol{Y}=\boldsymbol{y} ]$ を計算する.
\begin{align} Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(k)})&= E_{\boldsymbol{\theta}^{(k)}}[l^c (\boldsymbol{\theta},\boldsymbol{X})|\boldsymbol{Y}=\boldsymbol{y} ] \\ &= E_{\boldsymbol{\theta}^{(k)}}[ -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (Y_i - \mu)^2 + \sum^{m}_{j=1} (Z_j - \mu)^2 |\boldsymbol{Y}=\boldsymbol{y} ]\right\}\\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} E_{\boldsymbol{\theta}^{(k)}}[ (Y_i - \mu)^2|\boldsymbol{Y}=\boldsymbol{y}] + \sum^{m}_{j=1} E_{\boldsymbol{\theta}^{(k)}}[ (Z_j - \mu)^2 |\boldsymbol{Y}=\boldsymbol{y} ]\right\} \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (y_i - \mu)^2 + \sum^{m}_{j=1} E_{\boldsymbol{\theta}^{(k)}}[ (Z_j - \mu)^2]\right\} \end{align} ところで, $E_{\boldsymbol{\theta}^{(k)}}[Z]=\mu^{(k)}$, $E_{\boldsymbol{\theta}^{(k)}}[(Z - \mu^{k})^2]=(\sigma ^2)^{(k)}$ であるから, \begin{align} E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu)^2]&= E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu^{(k)}+\mu^{(k)}-\mu)^2] \\ &= E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu^{(k)})^2 +2(Z_j -\mu^{(k)})(\mu^{(k)} -\mu)+ (\mu^{(k)} -\mu)^2] \\ &= (\sigma ^2)^{(k)} +2E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu^{(k)})(\mu^{(k)} -\mu)]+ (\mu^{(k)} -\mu)^2 \\ &= (\sigma ^2)^{(k)} +2E_{\boldsymbol{\theta}^{(k)}}[(Z_j \mu^{(k)} -Z_j \mu -(\mu^{(k)})^2 -\mu\mu^{(k)})]+ (\mu^{(k)} -\mu)^2 \\ &= (\sigma ^2)^{(k)} +2( (\mu^{(k)})^2 - \mu^{(k)} \mu -(\mu^{(k)})^2 -\mu\mu^{(k)})+ (\mu^{(k)} -\mu)^2 \\ &= (\sigma ^2)^{(k)} + (\mu^{(k)} -\mu)^2 \end{align} よって, \[ Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(k)})= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (y_i - \mu)^2 + m (\mu^{(k)} - \mu)^2 + m(\sigma^2)^{(k)} \right\} \]
M ステップ
M ステップでは,$Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)})$ を最大化する $\boldsymbol{\theta}$ を求める.したがって, 以下の方程式をとけば良い. \begin{align} \boldsymbol{0}&= \frac{\partial}{\partial \boldsymbol{\theta}} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)})\\ &= \left( \begin{array}{c} \frac{\partial}{\partial \mu} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ \frac{\partial}{\partial \sigma ^2} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \end{array} \right) \\ &= \left( \begin{array}{c} \frac{1}{\sigma^2}\left\{\sum^{m}_{i=1}(y_i - \mu) + m(\mu^{(k)} -\mu) \right\} \\ -\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2} \left\{ \sum^m_{i=1}(y_i -\mu)^2+m(\mu^{(k)} -\mu)^2 + m(\sigma)^{(k)} \right\} \end{array} \right) \end{align} よって, \begin{align} \boldsymbol{\theta}^{(k+1)} &= \left( \begin{array}{c} \mu ^{(k+1)} \\ (\sigma ^2) ^{(k+1)} \end{array} \right) \\ &= \left( \begin{array}{c} \frac{1}{n}\left\{ \sum^{l}_{i=1}y_i + m \mu ^{(k)}\right\} \\ \frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2+m(\mu ^{(k)}- \mu ^{(k+1)} )^2 +m(\sigma ^2)^{(k)}\right\} \end{array} \right) \end{align} これが EM アルゴリズムの更新式となる.

実装
#データ
x <-c(54.62, 57.98, 49.41, 51.77, 46.72, 43.09, NA, 51.68, 40.33, 50.78, 43.63, NA, 44.89, 52.34, 47.28, 45.73, 49.02, 49.52, 46.41, 45.44, 46.50, 53.57, NA)

y <- x[!is.na(x)]
m <- sum(is.na(x))
n=length(x)

mu1=0; sigma1=1 #初期値
for(k in 1:100){
  mu2= (sum(y)+m*mu1)/n ;sigma2= (sum((y-mu2)^2)+m*(mu1-mu2)^2+m*(sigma1))/n #更新式
  if(abs(mu1-mu2) + abs(sigma1 -sigma2) < 1e-05) break #収束判定
  mu1=mu2
  sigma1=sigma2
}
結果,
> k
[1] 11
> mu2
[1] 48.5355
> sigma2
[1] 17.93601
$\mu=48.5355$, $\sigma^2 =17.9360$ となり, (k=11 で break しているということは,10回目と11回目がほぼ同じ値だったということなので)10回の反復で収束したことが分かる.実はこれらは欠測データを除いたデータの標本平均,標本分散と変わらない.
> mean(y)
[1] 48.5355
#R の関数 var はデフォルトだと不偏分散を求めるので…
> var(y)*(length(y)-1)/length(y)
[1] 17.936
ついでに,欠測観測数 m が 0,1,2,3 と変化するとき,収束状況がどのようになるか見てみよう.
theta =matrix(0, nrow=10, ncol=8)
for(m in 0:3){
  mu1=0; sigma1=1  #初期値
  for(k in 1:10){
    mu2= (sum(y)+m*mu1)/n
    sigma2= (sum((y-mu2)^2)+m*(mu1-mu2)^2+m*(sigma1))/n
    mu1=mu2
    sigma1=sigma2
    theta[k,(2*m+1):(2*m+2)]=c(mu2,sigma2)
  }
}
#print(xtable(theta),type="html")
結果,



m=0 m=1 m=2 m=3
k $\mu^{(k)}$ $(\sigma^2)^{(k)}$ $\mu^{(k)}$ $(\sigma^2)^{(k)}$ $\mu^{(k)}$ $(\sigma^2)^{(k)}$ $\mu^{(k)}$ $(\sigma^2)^{(k)}$
1 48.54 17.94 46.22 123.96 44.12 211.08 42.20 282.91
2 48.54 17.94 48.43 23.23 48.13 37.10 47.71 57.04
3 48.54 17.94 48.53 18.19 48.50 19.69 48.43 23.11
4 48.54 17.94 48.54 17.95 48.53 18.10 48.52 18.61
5 48.54 17.94 48.54 17.94 48.54 17.95 48.53 18.02
6 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.95
7 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94
8 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94
9 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94
10 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94

となり, m が大きくなるにつれて,収束が遅くなる傾向が見て取れる.
また,m = 0 のときは一回で収束している.これは m=0 のときの更新式が, \begin{align} \boldsymbol{\theta}^{(k+1)} &= \left( \begin{array}{c} \mu ^{(k+1)} \\ (\sigma ^2) ^{(k+1)} \end{array} \right) \\ &= \left( \begin{array}{c} \frac{1}{n} \sum^{l}_{i=1}y_i\\ \frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2 \right\} \end{array} \right) \\ &= \left( \begin{array}{c} \bar{y}\\ s^2_y \end{array} \right) \end{align} となるのだから,考えてみれば当たり前のことであった.

その他
更新式をみると,EM アルゴリズムとは,欠測 $\boldsymbol{z}$ の値を, $\boldsymbol{Z}$ の期待値 $E_{\boldsymbol{\theta}^{(k)}}[\boldsymbol{Z}|\boldsymbol{Y}=\boldsymbol{y}]$ で置き換えて計算する方法ではないことがわかる.
$\mu$ についての更新式 \[\mu^{(k+1)} = \frac{1}{n}\left\{ \sum^{l}_{i=1}y_i + m \mu ^{(k)}\right\}\] は, そのように考えたときと一致している. つまり $\mu^{(k+1)}$ と \[\frac{1}{n}\left\{ \sum^{l}_{i=1}y_i + \sum^{m}_{j=1} E_{\boldsymbol{\theta}^{(k)}}[Z_j|\boldsymbol{Y}=\boldsymbol{y}] \right\}\] は等しい. しかし, $\sigma^2$ についての更新式 \[(\sigma^2)^{(k+1)} = \frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2+m(\mu ^{(k)}- \mu ^{(k+1)} )^2 +m(\sigma ^2)^{(k)}\right\} \] は, \[\frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2 + \sum^{m}_{j=1}(E_{\boldsymbol{\theta}^{(k)}}[Z_j|\boldsymbol{Y}=\boldsymbol{y}] - \mu ^{(k+1)})^2\right\}\] とは別物である.

2012年4月28日土曜日

R: ヒストグラムの一部に色をつける

追記:書き直しました(2013/05/25)


ヒストグラムの一部に色をつける関数を書いた。
あまりエレガントなやり口とはいえないが……


使い方)
R のコンソールに以下をコピーして貼り付ける↓
 histcol <-function(data,part,haba, iro="#0080ff", honsu=20){
  xv0 =seq(part[1],part[2],haba)
  xv =sort(rep(xv0,2))
  xv =xv[-c(1,length(xv))]
  yv=numeric(0)
  for(k in  1:(length(xv0)-1)){
    yv  =c(yv, rep(length(data[data >= xv0[k] & data < xv0[k+1]]),2))
  }
polygon(c(xv, rev(xv)) ,c(rep(0,length(yv)), rev(yv)), col=iro ,density=honsu, border = NA)
}
引数 part に色をつける範囲を x 座標で入力、引数 haba にヒストグラムの幅を入力する。
ヒストグラムの幅はヒストグラムを描いてみて目で見て確認する。
iro は線の色、honsu は線の本数。


使用例)

x <- rnorm(100)
#png("survreg.png")
hist(x)
histcol(data=x,part=c(-1,1),haba=0.5)
#dev.off()
 


2012年4月21日土曜日

R: survreg の分布一覧

このページ(とこの本)では、パラメトリック生存モデルを当てはめる関数 survreg(survival パッケージ)が紹介されている。
そして指定できる確率分布として、ワイブル分布(デフォルト)、指数分布(exponential)、ガウス分布(gaussian) 、対数正規分布 (log-normal) 、ロジスティック分布(logistic)、対数ロジスティック分布(log-logistic)が挙げられている。
しかし実際には、以下の通り、極値分布(extreme)、レイリー分布(rayleigh)、t分布も使える。
> library(survival)
> names(survreg.distributions)
 [1] "extreme"     "logistic"    "gaussian"    "weibull"     "exponential"
 [6] "rayleigh"    "loggaussian" "lognormal"   "loglogistic" "t"     
これらの分布を使う必要がある場面は、たぶんあまり多くはないが、以下に使い方をメモして置く。

例)"Gehan" データにパラメトリックモデルを当てはめる。

パッケージ Fadist を読み込むと対数ロジスティック分布などが使える。
library(survival)
library(FAdist)
モデル当てはめ。
sfs =survfit(Surv(time, cens)~treat,data=gehan)

srs =survreg(Surv(time, cens) ~ treat, dist="exponential", data = gehan)
srs_r =survreg(Surv(time, cens) ~ treat, dist="rayleigh", data = gehan)
srs_w =survreg(Surv(time, cens) ~ treat, dist="weibull", data = gehan)
srs_g =survreg(Surv(time, cens) ~ treat, dist="gaussian", data = gehan)
srs_l =survreg(Surv(time, cens) ~ treat, dist="logistic", data = gehan)
srs_lg =survreg(Surv(time, cens) ~ treat, dist="loggaussian", data = gehan)
srs_ll =survreg(Surv(time, cens) ~ treat, dist="loglogistic", data = gehan)
srs_ext =survreg(Surv(time, cens) ~ treat, dist="extreme", data = gehan)
survreg の "extreme(極値)" はグンベル分布のことらしい。FAdist にもグンベル分布は入っているのだが、パラメータの置き方がよく分からなかったので、ここでは自分で定義することにする。
my_ext = function(t,a=1,l=0){
  1-exp(-exp((t-l)/a))
  }
プロット。
#png("survreg.png")
par(mfrow= c(2,4))

plot(sfs, main="exponential",lty=1:2)
curve(1-pexp(x, rate=1/exp(coef(srs)[1])) , add=TRUE, lty=1)
curve(1-pexp(x, rate=1/exp(coef(srs)[1] + coef(srs)[2])) , add=TRUE, lty=2)

plot(sfs, main="rayleigh",lty=1:2)
curve(1-pweibull(x,shape=1/0.5, scale=exp(coef(srs_r)[1]) ), add=TRUE, lty=1)
curve(1-pweibull(x, shape=1/0.5, scale=exp(coef(srs_r)[1] + coef(srs_r)[2])) , add=TRUE, lty=2)

plot(sfs,  main="Weibull")
curve(1-pweibull(x,shape=1/srs_w$scale, scale=exp(coef(srs_w)[1])) ,add=TRUE, lty=1)
curve(1-pweibull(x,shape=1/srs_w$scale, scale=exp(coef(srs_w)[1] + coef(srs_w)[2])) ,add=TRUE, lty=2)

plot(sfs,lty=1:2,  main="Gaussian")
curve(1-pnorm(x,sd=srs_g$scale, mean=exp(coef(srs_lg)[1])) ,lty=1, add=TRUE)
curve(1-pnorm(x,sd=srs_g$scale, mean=exp(coef(srs_lg)[1] + coef(srs_lg)[2])) ,lty=2, add=TRUE)

plot(sfs,lty=1:2,  main="log-Gaussian")
curve(1-plnorm(x,sdlog=srs_lg$scale, meanlog=(coef(srs_lg)[1])) ,lty=1, add=TRUE)
curve(1-plnorm(x,sdlog=srs_lg$scale, meanlog=(coef(srs_lg)[1] + coef(srs_lg)[2])) ,lty=2, add=TRUE)
 
plot(sfs,lty=1:2,  main="logistic")
curve(1-plogis(x,scale=srs_l$scale, location=(coef(srs_l)[1])),lty=1, add=TRUE)
curve(1-plogis(x,scale=srs_l$scale, location=(coef(srs_l)[1] + coef(srs_l)[2])),lty=2, add=TRUE)
 
plot(sfs,lty=1:2,  main="log-logistic")
curve(1-pllog(x,shape=srs_ll$scale, scale=(coef(srs_ll)[1])),lty=1, add=TRUE)
curve(1-pllog(x,shape=srs_ll$scale, scale=(coef(srs_ll)[1] + coef(srs_ll)[2])),lty=2, add=TRUE)

plot(sfs,lty=1:2,  main="extreme")
curve(1-my_ext(x,a=srs_ext$scale, l=(coef(srs_ext)[1])),lty=1, add=TRUE)
curve(1-my_ext(x,a=srs_ext$scale, l=(coef(srs_ext)[1] + coef(srs_ext)[2])),lty=2, add=TRUE)
#dev.off()
t分布については、よく分からなかった。

2012年4月13日金曜日

R: グラフの一部を塗りつぶす

いつもながら R-Tips からの孫引き

例)
標準正規分布の上側 5% 部分に色を付けたい時は以下のようにする。
 plot(dnorm, -4, 4)
 xvals <- seq(qnorm(0.95), 4, length=10)   # 10等分
 dvals <- dnorm(xvals)  # 対応するグラフの高さ
 polygon(c(xvals,rev(xvals)),c(rep(0,10),rev(dvals)),col="gray") 

polygon の引数 density に正の整数を与えると、塗りつぶしが斜線になる。
density=1 だと 1 本、density=2 だと 2 本の斜線が引かれる。
angle で斜線の傾きを指定できる(省略可。デフォルトだと 45°)。
 polygon(c(xvals,rev(xvals)),  c(rep(0,10),rev(dvals)),col="magenta",angle=60,density=20)


2012年4月3日火曜日

R: 空のリストを作る

追記(1/18/2013)
以下ぜんぶなしで、やっぱり単に l <-list() でやりましょう

追記(10/7)
最近はこれより lapply(1:4, function(i){numeric()}) でやることが多い

かんたんに空のリストを作る方法


例)
長さ4のリストを作る
> l <-mapply(rep,1:4,0)
> l
[[1]]
integer(0)

[[2]]
integer(0)

[[3]]
integer(0)

[[4]]
integer(0)

文字列や行列も代入できる
> l[[1]]="hello"
> l[[2]]= matrix(1:6, nrow=2, ncol=3) 
> l
[[1]]
[1] "hello"

[[2]]
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

[[3]]
integer(0)

[[4]]
integer(0)

2012年3月14日水曜日

R: sort と order の違い

sort() は、小さい順に整列させる関数、
order() は、『整列した各要素の元の位置』を返す関数。

データフレームに対して:order
ベクトルに対して:sort

だと思っておけば、間違えることは少ないはず……


例)データフレーム d を、列 a について整列させたい。

間違い
> d <-data.frame(ID =c("a","b","c","d","e","f"), a =c(1,4,3,6,2,5)) 
> d
   ID a
1  a 1
2  b 4
3  c 3
4  d 6
5  e 2
6  f 5
> sort(d$a)
[1] 1 2 3 4 5 6
> d[sort(d$a),]
   ID a
1  a 1
2  b 4
3  c 3
4  d 6
5  e 2
6  f 5
sort が返すのは 1, 2, 3, 4, 5, 6 という値なので、1行目、2行目、3行目、4行目、5行目、6行目、というふうに取り出されている。

  正解
> order(d$a)
[1] 1 5 3 2 6 4
> d[order(d$a),]
   ID a
1  a 1
5  e 2
3  c 3
2  b 4
6  f 5
4  d 6
(私は間違えたのです)


参考リンク:14. ベクトル計算

2012年2月25日土曜日

「母数」を「分母」の意味で使うのはやめろ

追記:「母数とはなにか」を書いた(6/19)


asahi.com(朝日新聞社):朝日新聞×河合塾 共同調査「ひらく 日本の大学」 - 教育 

今春卒業した大学生の就職率について、文部科学省は1日、「就職氷河期」と呼ばれた2000年春を下回り、過去最低の91・0%だったと発表している。文科省の調査は母数を「就職希望者」としているのに対し、今回の調査は卒業生全体を母数に取っている。
http://www.asahi.com/edu/hiraku/article01.html

難しい専門用語の使い方をまちがえることはよくある。これはしようがないと思う。
なるべく専門用語をつかわないようにして説明するとどうしても厳密さが失われる。これもしようがないと思う。
だけど、簡単に言えることをわざわざ難しく言い換えて、しかもまちがってるっていうのはだめだ。これはしようがなくないんだ。


だから「母数」を「分母」の意味で使うのはやめろ!!!!!!


これは「母数」の概念を理解しろとか、「サンプルサイズ」、「標本の大きさ」などの専門用語を使えとかいう主張ではありません。

上の記事を書いた方はおそらく、「母数」という言葉をどこかで見て、「分母」と同じような意味だと解釈したのだろう。だったら「分母」でいいじゃん。「分母」のほうが多くの人に通じるよ、という主張です。


参考文献
呉智英 『封建主義者かく語りき』 史輝出版


呉智英は「すべからく」を「すべて」の意味で使うことを「民主主義に由来する誤り」としてめちゃくちゃ厳しく批判している。こういう揚げ足取り、子どもっぽくて大好き。

2012年2月10日金曜日

Sweave 入門

文字コードについての記述を訂正 (2/11)

いきなり始めるエスウィーヴ(for windows)

口上
R で作業して TeX でレポートや発表資料をつくってる人は Sweave を使うと楽ができるかもしれない。Sweave は R に最初から入ってる関数。

入れるもの
"Sweave.sty" と "ae.sty" 。これを "Sweave.sty", "ae.sty" という名前で保存して、使っている TeX 環境の "なんとか.sty" っていうファイルがいっぱい入ってるフォルダ(私の場合は "C:¥w32tex¥share¥texmf¥tex¥platex¥base")に置く。

使用例
  1. ↓これを "test.Rnw" という名前で保存する

  2. \documentclass{jsarticle}
    \usepackage[dviout]{graphicx} 
    \usepackage{ascmac}
    \title{タイトル}
    \begin{document}
    \maketitle
    これはテスト
    
    <>=
     data(cars)
     plot(cars)
    @
    
    \end{document}
  3. R を起動して作業ディレクトリを "test.Rnw" のあるところに変更
  4. R コンソールに Sweave("test.Rnw",encoding = "cp932") と入力
これで以下のようなメッセージが出て来て、
> Sweave("test.Rnw",encoding = "cp932")
Writing to file test.tex
Processing code chunks ...
 1 : echo term verbatim eps
You can now run (pdf)latex on ‘test.tex’
('test.tex' に対して LaTeX を実行できます)
作業ディレクトリに "test.tex", "test-001.eps" というファイルができるはず。 "test.tex" の中身はこうなっている。
\documentclass{jsarticle}
\usepackage[dviout]{graphicx} 
\usepackage{ascmac}
\title{タイトル}
\usepackage{Sweave}
\begin{document}
\maketitle
これはテスト

\begin{Schunk}
\begin{Sinput}
> data(cars)
> plot(cars)
\end{Sinput}
\end{Schunk}
\includegraphics{test-001}

\end{document}
 
普通にコンパイルできる。するとこんな感じになる↓


以上です。

蛇足
  • "test.Rnw" の文字コード、UTF-8 にしたら日本語の部分が全部 "NA" になってた。Shift_JIS にしたら通った(Mac だとこの逆だったりするのかな?)
    「R-2.13.1から,関数Sweave()の引数にencordingを明示する必要があります.Shift-JISなら次のようにします.Sweave("test.Rnw", encoding = "cp932") 」ということだったようです。
  • 最初、"Sweave.sty" と "ae.sty" をブラウザから保存したら、html のタグもくっついてきてて、それに気付かずずっとエラーになってた。
  • "test.Rnw" の9行目、<< >> で囲まれた部分にいろいろオプションを付けられる。デフォルトだと pdf=TRUE, eps=FALSE になってる(グラフが pdf で出てきて、eps では出てこない)

参考・関連リンク
R日記の作り方 Sweaveを使う - ryamadaの遺伝学・遺伝統計学メモ
LaTeX - RjpWiki:R から LaTeX の表を出力したりする方法が載ってる。これでいちいち ... & ... & ... ¥¥ とか打たなくていい。

2012年1月4日水曜日

R のグラフ

グラフの見難さに定評があるわたくし
#png("colors.png")
cols = colorRampPalette(c("#0068b7","white","#f39800"))
iro10 =cols(10) 
curve(x^2,col=iro10[1],lwd=2.5,xlab="", ylab="")
for(i in 2:10)curve(i*x^2,col=iro10[i],add=TRUE,lwd=2.5,xlab="", ylab="")
#dev.off()

これで、少しはましになるかしら

参考リンク:
統計グラフの色
53. グラフィックスパラメータ(弐)

ggplot2 よさそう
これとか使いこなせたらすごいよさそう