2007年11月28日

gnuplotで日本語epsを作るための俺用メモ

この2日間の1/3くらいの時間をgnuplotの勉強に費やしてしまったが、 簡単なデータプロットならExcelかRで描いたほうがいい気もしてきた。まあ、2日程度で何が分かるというものでもないが。

ともあれ、思ったより面倒だった日本語epsの出力の仕方をメモっておこう。以下、 使用しているgnuplotはMaximaについてきたvar 4.2。

1. ターミナルの設定

set terminal postscript eps enhanced color "GothicBBB-Medium-EUC-H"

カラー&ゴシックで出力。フォントサイズも大きくした方が良いかもしれない。

2. 保存先の設定

set output "hoge.eps"

分かる場所に分かる名前で。

3. プロット

plot "datafile"

もしくは

load "scriptfile"

loadするスクリプトファイルのエンコードをEUCにしておくのを忘れない。

参考

posted by Rion778 at 01:03 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2007年11月24日

PCを使った数値積分(台形公式、シンプソンの公式)

20071124math

この定積分を原始関数を求めず行うことを 考えてみよう。 グラフで示すと、

20071124graph1

この水色の部分の面積を求めることに相当する。

まずは積分記号の意味を考える。

∫の下と上についた数字が積分の区間だが、この中のある数字を中身の数式に与える。帰ってくるのはy、つまり高さ。それにdx、 つまりほんのちょっとの幅を掛ける。そうすると、縦×横で面積が求まる。dxを無限に細い幅に設定し、 指定された区間のあらゆる値について面積を計算し、これをすべて足し合わせれば指定された区間の面積が分かるだろう。

で、原始関数さえ分かればこの果てしない作業が簡単になる。

が、原始関数というのはいつでもどこでも求められるとは限らないので、これに頼らず近時解を求める方法が必要となる。

まず思いつくのは、dxにそれなりに小さい適当な値を与え、長方形の面積の和として近似的に解を求める方法だろう。 イメージとしてはこう。

 20071124graph2

なんとなく近い値が出そうだが、誤差のあることも分かる。ともかく、Rで計算してみよう。区間分割の幅は0.1とする。

a <- 0
x <- 2
while(x < 5.9){
 a <- a + f(x)*0.1
 x <- x+0.1
 }

aの中に解が入る。値は 8.253492 だった。

次に考えるのは、長方形で無い多角形の面積の和として解を求める方法。例えば台形なら公式から簡単に面積が求められるので、 台形を採用することにする。イメージとしてはこう。

20071124graph3

見るからにさっきよりも誤差が少なくなりそうな感じ。さっきと同じ区間分割で計算してみる。 具体的にはaにおけるf(x)の値とa+dxにおけるf(x)の値を上底・下底とし、区間幅を高さとし、「(上底+下底)×高さ/2」 という例の公式を使う。

b <- 0
x <- 2
while(x < 5.9){
 b <- b + (f(x)+f(x+0.1))*0.05
 x <- x+0.1
 }

結果は 8.411346。

しかしグラフにはまだ少し隙間がある。もっと精度を上げられないだろうか。例えば、多角形の「頭」に「丸み」を持たせられたら、 もっとぴったりさせられるだろう。イメージとしては…めんどいのでパス。まあ、とにかくさっきよりさらにぴったりすると。

で、「丸み」を持たせるには、「頭」を関数で表現する必要がある。前回説明したラグランジュ補間を使うと、 3つの点から二次関数が作れるので、これを利用して「頭」を二次関数にする。

こうすると縦×横みたいな簡単な計算では面積が求められないので定積分を利用する必要がある。しかし、 区間幅はこちらで設定するので既知であり、二次関数の定積分ということも分かっているので、先に計算をある程度終わらせておくことができる。

a,m,bという三点(a<b、mはaとbの中点)と、それに対応するf(x)の値、つまりf(a)、f(m)、 f(b)が分かっている場合、区間aからbまでの積分は次の式(シンプソンの公式)で近似される。

20071124math2

ここでh=m-a=b-mである。ラグランジュ補間して定積分するとちゃんと出てくるので暇な人は確認を。

ともかくRで計算してみよう。用いるf(x)の値の数はさっきまでと同じだが、 3点を使って区間を定めるので分割された区間の幅は倍になる。

c <- 0
x <- 2
while(x < 6){
 c <- c + (f(x) + 4*f(x+0.1) + f(x+0.2))*0.1/3
 x <- x + 0.2
 }

結果は 8.414278 。

で、これまで三つの値が出た。

  1. 長方形による近時 …8.253492
  2. 台形公式             …8.411346
  3. シンプソンの公式  …8.414278

本当に真の値に近づいていっているのだろうか。

手計算をして確認すればいいのだが、面倒くさい。

実はRには数値積分をしてくれる関数があるのでそいつを使おう。

f <- function(x) sin(x+2.2) - cos(x^2.3/factorial(3)) +2
integrate(f, 2, 6)

integrate(関数, 積分範囲の下限, 積分範囲の上限)とするだけで答えが返ってくるのだから簡単なものだ。で、 返ってくるのは

8.414256 with absolute error < 4.2e-08

後ろの英語は4.2×10-8以下の小さな誤差が含まれることを言っている。先ほどの3つの値と比較すると、 1→2→3の順に精度が良くなっていることが分かるだろう。特にシンプソンの公式は非常に精度が高いことが分かる。

posted by Rion778 at 12:16 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月23日

Rによるラグランジュ補間

今日は数値積分を調べてたら一日終わってしまった。

区間を細かく分割し、分割された区間の面積をそれぞれ計算して全部足し合わせる、というのがPCを使う場合の数値積分の基本的な方法。 要は力技だ。

それで、区間を細かく分割するほど精度は高まるのだけど、 細かく分割された区間の面積をどう計算するのかによっても数値積分の精度は異なってくる。

例えば、関数がいくらかでも斜めの部分を持つならば長方形で計算するよりも台形で計算したほうが精度は良くなるだろう(参考: Wikipedia 「台形公式」)

そして、関数がいくらかでも滑らかさを持つならば同じように滑らかなもの、 例えば二次関数を利用してやると精度はさらに良くなるだろう。その方法にシンプソンの公式(参考:Wikipedia 「シンプソンの公式」)と呼ばれるものがある。こいつを理解するには、任意の3点から二次関数を導く方法を理解しなければいけない。

そこで用いられるのがラグランジュ補間という方法。(n+1)個の点があったとき、それらを通るn次関数を求める。

手法の解説は以下を参考に。簡潔で分かりやすいです。

で、上のページにはC言語によるラグランジュ補間のソースコードが置いてある。Cのが早いんだろうとは思いつつ、 どうにかRに移植できないか考えてみた。しかしそのままRに書き写してみたところで上手くいかない。ところどころ値が飛んでしまう。

仕方ないので数式を見つつ最初から組み立てる。

ソースの前に適当な「点」を用意。

set.seed(77)
x <- c(1:7)
y <- runif(7)*20-10

今回は乱数の「種」を指定したので同じようにRに打ち込めば同じ結果が返ってくると思う。

で、以下がラグランジュ補間。

X <- numeric(0) 
Px <- numeric(0)
n   <- 0
for(t in seq(min(x),max(x),0.01)){
n <- n+1
px <- 0   
for(i in 1:length(x)){
Ux <- 1
Dx <- 1
for(j in 1:length(x)){
 if (i==j) Ux <- Ux
 else Ux <- Ux*(t-x[j])
 }
for(j in 1:length(x)){
 if (i==j) Dx <- Dx
 else Dx <- Dx*(x[i]-x[j])
 }
 px <- px + y[i]*Ux/Dx
 }
 X[n] <- t
 Px[n] <- px
}

これで、Pxの中にラグランジュ補間で得られた値が格納される。Xと共にプロットしてやればいい。 ついでに元になったyの値を赤丸でプロット。

plot(X,Px,type="l",ylim=c(-20,20))
points(x,y,col="red")
grid()

ここまでのスクリプトをすべて実行すると次のグラフが得られる。

lag

要するにこういう風に点と点の間を繋いでくれる。この場合曲線は6次関数である。 ちなみにxはソートしてあれば等間隔じゃなくてもいい。

そして3点ならば二次関数で繋ぐ。例えば

x <- c(1:3)
y <- c(5,3,9)

としてさっきのスクリプトを動かし、プロットすれば、

となる。

こんな感じで3点から二次関数を作り、その間を定積分するのがシンプソンの公式と呼ばれる方法なんだけどそれはまた次回にでも。

posted by Rion778 at 20:19 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月19日

中心極限定理(2)

前回より

平均μ、分散σ2である何らかの分布からサンプリングされたデータの平均値は、平均μ、 分散σ2/(サンプリング数)の正規分布に従う。

これが中心極限定理であった。

これをさらに詳しく確認するためには、 実際に分散がσ2の母集団からn個のサンプリングを繰り返し、平均をいくつも計算し、 その分散をσ2/nと比較すればいい。

n=1の場合、n=2の場合、 n=3の場合とひとつずつnを増やしながら平均値の分散をシミュレーションによって生成していき、 理論値であるσ2/nとの比を逐一計算していき、nを横軸に、 理論値と実測値の比を縦軸にプロットして行ったらより視覚的に中心極限定理を確認できるだろう。もちろん、 値は1を中心に集まるはずだ。

では、例によってRでスクリプトを書いて確認。

まず、母集団を適当に決める。分散と平均値が決まっていればなんでもいいので、 今回は自由度3のカイ二乗分布としよう。確率密度は次のグラフのようなものである。

chisq

で、この分布に従う乱数を発生させるには、rchisq関数を使う。 たとえば10個生成するなら

> rchisq(10,3)
 [1] 0.6557817 6.0603239 0.4138968 1.1510150 1.7536676 1.8107585 4.7430857
 [8] 1.3513523 3.9005390 7.6347070

詳しくはヘルプを参照のこと。

それで、 やりたいことはサンプリング数を増やしながら平均値の分散を調べることだった。そのために次のようなスクリプトを書いてみた。

mean.var <- numeric(0)
for(a in 1:1000){
 y <- numeric(0)
 for(i in 1:300){
 z <- rchisq(a,3)
 y[i] <- mean(z)
 }
 mean.var[a] <- var(y)
 }

for文の中でfor文を回している。

サンプル数1で300回サンプリングし、 それぞれの平均(サンプル数1の場合はそのままの値だが)を計算、その平均値の分散を計算して、 その値をmean.varのa番目に格納。そうしたら今度はサンプル数を1つ増やして2にして300回サンプリングし、… 以下サンプリング数1000まで繰り返し。というようなことをやっている。

当然だが馬鹿みたいに重いので注意。 私が大学に持っていっている旧世代のVAIOではメモリが足りなくなるかもしれない。

それで、このスクリプトを回して得られるのは要するに何か。それは、 標本の平均値が真の平均値に対してどの程度ばらつくのかという指標と、 その指標がサンプル数に応じてどう変化するのかという2つの情報である。とりあえずプロットしてみよう。

 mean_var

x軸がサンプル数、y軸が平均値のばらつきの程度(分散)である。 サンプル数を増やすと、急激にばらつきが減少する、つまり真の平均値へ近づきやすくなっているということがわかるだろう。ちなみに、 母集団、つまり自由度3のカイ二乗分布における平均値は3、分散は6である。そして、 たとえばサンプル数200の時の平均値の分散は0.02936468であり、 6/200=0.03に近いことからも中心極限定理が確認できる。

なお、 確率密度関数のわかっている確率変数は積分によって平均と分散を正確に知ることができる。Rでやるなら次のようにする。

f <- function(x){
 y <- x*dchisq(x,3)
 return(y)
 }
integrate(f,-Inf,+Inf)
g <- function(x){
 y <- ((x-3)^2)*dchisq(x,3)
 return(y)
 }
integrate(g,-Inf,+Inf)

さて、 やりたいことは平均値の分散が母集団の分散÷サンプル数と等しくなるのかどうかの確認であった。だから、 先ほど計算した各nについての分散の値を6/nと比較してやればいい。つまり、

mean.var/(6/1:1000)

を計算するのだ。この計算では1000個の答えが返ってくるわけだが、 中心極限定理が正しいなら、そのすべての値が1を中心に分布するはずである。それではプロットしてみよう。

var

中心に1の直線を赤で引いた。1を中心に値が集中していることが分かる。

まとめ

さて、今日分かったことは何だったか。

まず一つ、サンプルから計算される平均値は、 サンプル数を増やすと急激に真の平均値へ近づくということである。それは、今日の2つ目のグラフから分かる。 あのグラフではサンプル数1000の場合までプロットしているので分かりにくいが、 100の場合までにしてプロットしなおすと次のような形になる。

mean_var2

このグラフから分かるように、 サンプル数を20程度まで増やすうちは平均値の推定値の精度は急激に高められる。しかし、 40以上に増やしたところでさほど意味のないこともわかる。 この関係は実際に実験でサンプル数を決定する際の指標となるので重要である。

次に、サンプルから計算される平均値のばらつきの程度は、 母分散とサンプリング数から計算することが可能であるという点である。母分散とサンプル数から計算された平均値の分散が、 実際に観測された平均値の分散の値とよく一致することは3つ目のグラフからわかる。なお、母分散はふつう未知であるため、 その推定値である不偏分散(標本から計算される)を用いることが多い。

以下関係のない話。

僕の知りたいのは調和平均のバラツキの仕方なのだ。 サンプル数を増やすと真の値に近づくということ、分散が小さくなるということは分かっている。でも、 算術平均の時のような簡単な関係が見つからない。「母分散/サンプル数」では実際の分散より小さくなってしまう。 どうも調和平均はサンプル数が平均値へ影響を与えてしまったりしてややこしい。 シミュレーションによる力技では有効な推定値を与えられそうにない。誰かすでに答えを持ってたりしないだろうか。

posted by Rion778 at 22:29 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月14日

Rを使った1階常微分方程式の数値解の求め方(Euler法)

1階の常微分方程式 ってーのは例えばこういうもの。

dydt

この式の場合、yの変化量がyそのものによって定義される。

このままだと、「tが1増えるときにyはyだけ増加する」と言っているに過ぎない。実際に数値解を得るときは、例えば 「t=0の時y=1である」みたいな初期条件を与え、その上で「t=10のとき、yはいくつになっているか?」と問うわけだ。

t=0でy=1、増加量はyの値なのだから、y=1という条件より、単純に「y=t×1+1」として計算していいのだろうか? 答えは11なのだろうか?

そんなのでいいはずがない。

例えば、いきなりt=10の時のyを計算せず、t=5の時のyを計算し、そのyを使ってt=10の時のyを計算したらどうなるか。

まず、t=5の時を考える。とりあえず「y=t×1+1」として計算すると、yは6になる。

次に、t=10の時を考える。最初の微分方程式より、変化量はyの値であるので、今度は「t×6」を使わなければいけない。 tはさっきの5から5だけ増えているので、増加量は30。さっきの6と足してyの値は36だ。 間に1ステップ置くだけで値は3倍以上になってしまった。

同様に、tを細かくとればとるほどに値は増えていく。それでも、何かしらの値には収束していく。tの区間を「無限に」 細かくとれば真のyの値が分かるだろう。

が、まずは「力技」で近似解を求めてみよう。方法としてはこうだ。

  1. yとtについて初期値を与える
  2. tをほんの少しだけ増やす。
  3. tをほんの少し増やした時のyの増分を計算し、最初のyに足す。
  4. 再度tをほんの少しだけ増やす。
  5. tをほんの少し増やした時のyの増分をさっきのyを使って計算し、さっきのyに足す。
  6. 再度tをほんの少しだけ増やす。
  7. 以下好きなだけ繰り返し

では、この方法でもって最初に示した微分方程式を解くプログラムをRで書いてみよう。初期値はt=0、y=1とする。

y <- 1
t <- 0
i <- 1
Y <- numeric(0)
T <- numeric(0)
while(t <= 10){
 Y[i] <- y
 T[i] <- t
 y <- y + y*0.01
 t <- t+0.01
 i <- i+1
 }
plot(T,Y,type="l")

ここでの「ほんの少し」は0.01。とりあえずt=10まで計算を続け、プロットしてみた。結果はこう。

Euler

最初にやっていた10やら36どころの数字ではない。t=10におけるyの値は20959.16である。

この「力技」のことをEuler法とも呼ぶ。やたらいろんなものに出てきますね、オイラー先生の名前は。

なお、方法は他にもある。例えばRunge-Kutta法だとか。それらに比べるとEuler法は精度がボチボチなんだけど、 簡単な方法なのでこうした説明にしばしば使われるのです。

ところで、最初の微分方程式は変数を分離して積分することにより割と簡単に解ける。つまり、y=f(t)の形にできる。 手順は省略するけど、t=t1の時のyをy1、t=t2の時のyをy2とする(t1<t2)と、

solvedydt

例えばコイツにy1=1、t1=0を初期値として与え、plot関数でプロットしてみると次のようなグラフが描ける。

fx

さっきのグラフとほぼ同じ。

んがっ、t=10の時の値は22026.47である。

1000くらいの誤差があったわけです。

ちなみに、「tをほんの少し増やす」の「ほんの少し」をさっきの1/10、0.001にして計算すると21916.68となる。 この辺りはマシンパワーとの相談。

参考文献

posted by Rion778 at 23:02 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月13日

チロルチョコタワー

1+2+3+4+...+99+100=?

この問題を解くには、

100+99+...4+3+2+1

という式を足してやると簡単。つまり、101が百個の半分、答えは10100/2=5050。

これを一般化すれば、1から任意の自然数nまでの数字を一つずつ増やしながら全て足し合わせた時の答えはn(n+1)/2。

そしてこれは、三角に積み上げられたブロックの個数を求めるのにも使える。

例えば、底辺から10個、9個、8個…2個、1個とブロックが積み上げてあったなら、 その三角のタワーは10×11/2=55個のブロックから作られていることになる。

つまり、すでに出来上がっている三角タワーのブロック数を知るのは容易い。

ここまでは普通の算数。

では、任意の数のブロックを用いたとき、最大で何段のタワーが作れるのか、という問題を作ってみよう。「余った」 ブロックは如何様にも使えないものとする。

そうすると少し問題は複雑になる。

とりあえず、Rでプログラムを書いて力技で解いてみよう。考え方としては、こう。

  1. まずブロックを一つ取る。
  2. 残ったブロックが、次の段に必要なブロックの数以上であるか判断。
  3. 足りなければ終了。出来上がった段数を調査。
  4. 足りていたら、次の段に必要なだけブロックを取る。
  5. 2に戻って終わるまで繰り返し。

では、Rのスクリプト。

my.pyramid <- function(x){
 b <- 0
 while(x > b){
b <- b+1
x <- x-b
}
 return(b)
 }

繰り返しの条件が「次の段に必要なブロック数以上であるか」ではなく、「今の段に必要なブロック数より大きいか」になっているけど、 結果は同じ。実行結果はこちら。

> my.pyramid(1)
[1] 1
> my.pyramid(3)
[1] 2
> my.pyramid(5)
[1] 2
> my.pyramid(7)
[1] 3
> my.pyramid(9)
[1] 3
> my.pyramid(10)
[1] 4

うん。ちゃんと結果が出ている。

でもこれだとプロットするのに都合が悪い。ちゃんと方程式を解いて一般的な求め方を出そう。

まず、段数をx、ブロック数をnとすると、n=x(x+1)/2であることは最初に示したとおり。とりあえずこれをxについて解く。 めんどいのでMaximaで。

2007y11m13d_224532620

うーん、便利だな。Maximaは。

で、xは当然0より大きいので、必要なのは2つ目の解。こいつをRで関数にしてやる。

my.pyramid2 <- function(x){
 n <- (sqrt(8*x+1)-1)/2
 return(floor(n))
 }

returnのところで使っているfloor関数は小数点以下を切り捨ててくれる。 切り捨てないと2.3段とか変な答えが返ってきてしまうので少し細工をしました。

実行結果はさっきの関数と同じになるので省略。気になる人は各自確認を。

で、こいつをプロットしてみよう。ブロックの個数が増えるに従い、どんな変化をするのだろうか。100個くらいまで見てみよう。

plot(my.pyramid2,0,100)

pyramid

こんな感じ。100個つかったって12段にしかならない。

ところで、なんでこんなことをしていたかと言いますと、チロルチョコがたくさんあったんですね。で、それを積み上げるアホウが一人。 そのアホウが言うわけです。

「このチロルチョコ全部つかったら何段のピラミッドが作れるかなぁ」

posted by Rion778 at 22:48 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月12日

四分位数とヒンジ

データを大きさ順に並べたときに、小さいほうからデータ全体の1/4が含まれるような順位の値を第1四分位数、 同様に3/4が含まれるような順位の値を第3四分位と言う。第1四分位数は下側四分位、第3四分位数は上側四分位とも言う。なお、 第2四分位数は中央値に等しい。

一方、中央値以下のデータの中央値を下側ヒンジ、中央値以上のデータの中央値を上側ヒンジという。

下側四分位数と下側ヒンジ、上側四分位数と上側ヒンジは大体同じで、 特に区別されることなく用いられている場合もあるが実際は異なるものであり、対象となるデータによっては値がいくらか異なる。

例えば、次のような1から10までの数値が一つずつ含まれるデータ

1, 2, 3, 4, 5, 6, 7, 8, 9, 10

について下側四分位、上側四分位を計算するとそれぞれ3.25、7.75となるが、下側ヒンジ、上側ヒンジの値はそれぞれ3、 8である。

一方、これに0を加えて次のようなデータにすると、

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

四分位数とヒンジの値はそれぞれ2.5、7.5で等しくなる。

ヒンジの定義は中央値以下(以上)のデータの中央値であるので、対象となるデータに中央値が含まれない場合、 中央値からヒンジまでの距離は最小値(最大値)からヒンジまでの距離より「遠く」なる。

一方、四分位数の場合は下側(上側)四分位数から中央値までの距離と最小値(最大値)は常に等しい。

ところで四分位数はデータを並べたとき下から25%、75%の位置にあるデータの値を述べている。 他のパーセントについて述べるときは、そのパーセントをqとしてqパーセンタイルなどと言う。

では、任意のqパーセンタイルを計算するのは具体的にどうするのか。 ヒンジは中央値から最大(最小)値までに含まれるデータでもって再度中央値を計算するだけだから簡単だが、 厳密にq%の位置にあるデータというのはちょっと計算が面倒に思える。

とりあえずWikipediaを覗いてみる(数式の画像はこちらで新たに用意した)。

n 個のデータ x に対する q 分位数 Qq は、 昇順にソートしたデータを xs とすると、
qq
x 

     _, ._
  ( ゚ Д゚)   …………

えー。まずは上から。

(1-q+qn)というのは、(n-1)q+1を整理したものかと。つまり、 (n-1)で全ての個数を0をはじめとする順序(0,1,2...n-1)に対応させ、 qを掛ける事でその中の何番目がq%点に対応するかを計算、再度1を足すことで1をはじめとする順序(1,2,3...n)に戻していると。

で、下の式。

もしも(1-q+qn)が自然数であるならば、 その順位に対応したデータxが存在するということなのでそれをそのままQqとして使える。

しかし、自然数で無い場合、その順位をはさんでいる2つの数字から計算すると。 具体的に計算したほうが分かりやすい。

例えば、

1, 2, 3, 4, 5, 6, 7, 8, 9, 10

という10個の数字からなるデータがあったとして、このデータの25%点を求めるとする。

まずは25%点がいったい何番目のデータに該当するのかを計算。(1-q+qn)に代入して

1-0.25+0.25*10 = 3.25

よって3.25番目のデータが25%点である。

で、3.25は自然数ではないので3番目のデータである3と4番目のデータである4の間を0.25: 0.75で分ける。

(4-3.25)*3 + (3.25-3)*4 = 3.25

よって、25%点は3.25である。

ちなみに、Rには分位数を求める関数のquantileがある。データベクトル以外の引数を特に与えなければ、 最小値、下側四分位、中央値、上側四分位、最大値を返してくれる。特定のパーセンタイルを知りたければ、その値を引数として与える。

> x <- 1:10
> x
 [1]  1  2  3  4  5  6  7  8  9 10
> quantile(x)
 0%   25%   50%   75%  100% 
 1.00  3.25  5.50  7.75 10.00 
> quantile(x,p=0.05)
5%
1.45

ヒンジを知りたいときはfivenum関数が便利。最小値、下側ヒンジ、中央値、上側ヒンジ、 最大値をベクトルとして返してくれる。

> fivenum(x)
[1]  1.0  3.0  5.5  8.0 10.0

ヒンジと四分位数はこのように似ているようで少し違うので、データ数が少なく、 奇数個であるような場合は少し注意したほうがいい気がしないでもないでも。どっちがいいかと言われても知らないけど。

ちなみに、五数要約や箱ヒゲに使うのはヒンジ。なんでだろ。

参考文献

posted by Rion778 at 21:24 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月11日

サケフレークとタマネギのパスタ

DSC 食べ物が無くなって来たころに適当に作ったものって意外とおいしかったりするよね。

材料(1人前)

  • パスタ         …1人分100g
  • サケフレーク     …ビンの1/3(大さじ2くらい?)
  • タマネギ         …中1/4個
  • マヨネーズ       …大さじ3
  • ピエトロドレッシング …大さじ1
  • サラダ油        …大さじ1
  • 塩、コショウ      …適量

プロトコル

  1. パスタを茹でる。アルデンテ。
  2. 茹でている間にタマネギをスライス。
  3. パスタが茹で上がる1分ほど前に熱したフライパンへサラダ油(大さじ1)を入れ、タマネギのスライスを炒めはじめる。
  4. 少ししんなりしてきたところでパスタを投入。塩、コショウで軽く味付け。
  5. さっさとサケフレーク投入。よく混ぜる。
  6. ある程度タマネギに火が通ったら火を止める。
  7. マヨネーズ、ピエトロドレッシングを絡め、完成。
posted by Rion778 at 15:59 | Comment(2) | TrackBack(1) | レシピ | このブログの読者になる | 更新情報をチェックする

分散分析

分散分析とは、データのバラつきの仕方、すなわち分散を調べることにより、グループ間に平均値の差があるのかどうかを調べる手法だ。

誤差によるデータのバラつき(グループの平均値から個々のデータはどの程度バラついているのか?)と、 各グループの平均値のバラつきを比較し、グループ間におけるバラつきの方が非常に大きいときは「グループ間には偶然によらない差がある」 との判断を下す。

具体的には何をどうやっているのか。例によってRを用いて確認していこう。

まずはデータの準備。

分散分析の場合はデータが正規分布に従っており、分散が等しいという前提が必要なので(つまり、 現実の事例に適用するには注意が必要)、分散が1、平均値のみが異なる2種のグループを用意する。正規分布に従った乱数を発生させる関数、 rnormを使用し、次のようにそれぞれ10のデータを用意する。

> x <- rnorm(10,5)
> y <- rnorm(10,1)
> x
 [1] 3.698867 3.440333 3.651131
 [4] 3.882890 5.221521 7.689116
 [7] 5.145143 6.630719 4.197294
[10] 4.613909
> y
 [1]  1.6828048  1.6255364  0.8231488
 [4]  2.5112012 -0.2590962 -0.5415252
 [7]  0.8268164  0.6663882  1.7391940
[10]  0.7674629

グループxは平均5、グループyは平均1。分散分析によって差があると判断されることを期待しよう。なお、 あくまで乱数なので以下の具体的な計算結果などは必ずしも再現性のあるものではない(ちなみに、 乱数の再現性を保障するためには乱数の種を固定してやればいい、けど書き直すの面倒なんで今回はパス)。

まずは各グループの平均値が、互いにどのくらいバラついているのかを計算する。バラつきの程度を判断するには、 平均値からの距離を計算すればいい。つまり、偏差「グループの平均−全体の平均」を計算し、データの個数だけ足し合わせればいい。 しかしそのままだと0になってしまうので、「(グループの平均−全体の平均)^2」を足し合わせたもの、 つまり偏差平方和をバラつきの程度の判断指標として採用し、これを計算することにする。

> a <- sum((mean(x)-mean(c(x,y)))^2*length(x) + (mean(y)-mean(c(x,y)))^2*length(y))
> a
[1] 73.45558

mean(x),mean(y)がそれぞれのグループの平均値、mean(c(x,y))が全体の平均値。 length関数はベクトルの長さを返す。つまり、length(x)はグループx内のデータの個数である10を返すといった具合。

次に、同じような要領で誤差によるバラつきの程度を計算しよう。ここでいう「誤差によるバラつき」というのは、 個々のデータが自分の所属するグループの平均値からどの程度離れているのか、ということを意味している。よって、「(個々のデータ− グループの平均)^2」を足し合わせる。

> b <- sum((x-mean(x))^2 + (y-mean(y))^2)
> b
[1] 25.48589

これで、グループ間、誤差、2種類のバラつきをあらわす偏差平方和が得られた。しかしこのままではデータの個数、 グループ数の影響が残っているので、データの個数を反映した「何か」で割ってやる必要がある。その「何か」は「自由度」と呼ばれるもので、 基本はデータ数なのだけれど、平均値を計算するたびに1つずつ減っていく。

なんで平均値を計算すると減るのか。それは、平均値が決まってしまうと、自由に動かせるデータが一つ減るからだ。例えば、 3つの数字からなるグループがあるとしよう。平均値が何なのか決まっていなければ、その3つにどんな数字が入っていようが問題はない。 しかし、「平均値は2である」という情報が与えられると、2つまでは自由な数字を入れられるけど、 3つ目は平均値が2になるような数字を選ばなければならなくなる。つまり、平均値という情報がデータの自由度を一つ奪ったのである。

少し話がそれたが、偏差平方和をこの自由度で割ったものを不偏分散という。不偏分散はデータの個数に影響されることなく、 データのバラつきを表してくれる。だから、この不偏分散を比べればバラつきかたの違いを知ることができる。それでは、 グループ間の不偏分散を誤差の不偏分散で割って(つまり、 グループ間は誤差に比べて何倍バラついているのかを計算する)分散比を求めてみよう。

> (a/(2-1)) / (b/(10+10-2))
[1] 51.87971

a、つまりグループ間の平方和は、xの平均、 yの平均という2つのデータから全体の平均を引いているので自由度は2-1で1となる(ここでのxの平均、 yの平均は自由度から引く個数に加えないことに注意。「xの平均」「yの平均」はそれぞれ「データ」なのであって、この「データ」 の自由度が「全体の平均」を計算することで奪われる)。そしてb、誤差平方和の場合は、 すべてのデータ(20個)から2つの平均値を用いて算出した値であるので、自由度は20-2で18となる。

で、実際計算すると52倍くらいとなった。やはり相当違うようだ。グループ間がこれだけバラついているということは、 グループ間に差があると言えそうだ。

でも、もしかしたら偶然に52倍の差が出ているのかもしれない。偶然なのかどうかを確認するために、 「F分布」を利用する。分散がまったく同じ母集団からサンプリングしたデータでも、 分散の比をとってみると多少はバラつく。そして、分散が同じ母集団からサンプリングされたデータの分散の比は、 F分布に従ってバラつくことが知られている。そして、F分布は分散比を計算するときに使われた2つの分散の自由度によって形が決定される。 今回の場合、グループ間の分散を計算したときの自由度は1、誤差の分散を計算したときの自由度は18なので、「第一自由度1、 第二自由度18のF分布」を調べれば、51.87971という分散比が偶然によるものか否かを判断することができる。

まずはそのF分布を描いてみよう。

curve(df(x,1,18),0,6,ylim=c(0,1))

df

4倍くらいまでの分散比は割と出そうだが、それ以上の分散比はあまり出そうでないことがグラフより見て取れる。 先の分散比は52倍もあるのだから、偶然と言えそうも無いことはもはや明らかだ。

しかし念のため、F分 布の上側95%点、つまり、それ以上の分散比は5%以下の確率でしか出ないという点を計算し、 その値よりも先ほどの分散比が大きい場合は、そんなことはめったに起こらないのだからグループ間はバラついている、 差があるのだということにしよう。

> qf(0.95,1,18)
[1] 4.413873

4.4程度だった。先ほどの分散比はこれよりも圧倒的に大きいので、グループ間に差があるという判断を下そう。 これで分散分析は終了である。

ちなみに、分散分析を簡単にやってくれる関数も当然ある。oneway.test関数を使うなら、

> xg <- rep(1,10)
> yg <- rep(2,10)
> xy <- data.frame(data=c(x,y), group=as.factor(c(xg,yg)))
> oneway.test(data~group, xy)
One-way analysis of means (not assuming
equal variances)
data:  data and group 
F = 51.8797, num df = 1.000, denom df = 15.706, p-value = 2.344e-06

となる。p-valueというのは分散比51.87971が偶然に出てくる確率である。非常に小さいので、 これは偶然とは言いにくい=グループ間に差があると判断して差し支えないことが分かる。

ところで、分散比は本当にあのような分布でバラつくのだろうか。シミュレーションをしてみよう。平均値、 分散ともに1の標準正規分布から一つは2個、一つは19個のデータをサンプリング、不偏分散を計算(つまり、自由度は1と18)し、 その比を計算するという作業を5000回ほど繰り返してデータを集め、ヒストグラムを描いてみよう。そして、 その上に自由度(1,18)のF分布を重ね描きしてみる。(truehist関数を使っているのでMASSパッケージをお忘れなく)。

z <- numeric(0)
for(i in 1:5000){
 x <- rnorm(2)
 y <- rnorm(19)
 z[i] <- var(x)/var(y)
 }
truehist(z,xlim=c(0,6),ylim=c(0,1))
par(new=T)
curve(df(x,1,18),0,6,ylim=c(0,1))
quantile(z,0.95)
qf(0.95,1,18)

df2

シミュレーションの結果はきれいにF分布の密度関数に重なった。めでたしめでたし。 

参考文献

posted by Rion778 at 04:09 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2007年11月07日

中心極限定理

中心極限定理によれば、元の分布が何であれ、 そこからサンプリングされた標本の平均値は正規分布に従って分布する。

にわかには信じがたい話だ。例えば、こんな分布(ちなみに、自由度2のカイ二乗分布)

shisq

のデータからサンプリングを繰り返し、平均値をいくつも計算する。そうすると、計算された平均値たちはこんな分布

norm

にしたがってバラツクと主張しているのだから。

ともあれ、Rを使って実際に確かめてみよう。

まずは母集団を何にするか。Rを使えばカイ二乗分布などに従う乱数を簡単に発生させられるけど、 そんなのはみんなやってるからもっと変なものがいい。

というわけで、いま僕が思いついたランダムでもなんでもない適当な数字を寄せ集めを母集団とすることにしよう。

> x <- c(1,1,2,3,4,5,5,5,6,7,7,8,8,8,8,9,9,0,0,0,0,0,-1,-2,-4,-5)

こいつの分布を、hist関数でヒストグラムにしてみる

> hist(x)

すると、

myhist

こんな感じにバラついた分布が出来上がった。平均値と分散を計算しておこう。

> mean(x)
[1] 3.230769
> var(x)
[1] 17.30462

不偏分散だけど…まあいいか。それではこの適当な母集団からサンプリングと平均値の計算を繰り返し、 ヒストグラムを描いていくことにしよう。結局やることは前のブートストラップとほとんど同じ。ソースの雛形はこんな感じで。

z <- numeric(0)
for(i in 1:5000){
 y <- sample(x,「サンプリング数」,replace=T)
 z[i] <- mean(y)
 }
truehist(x)

「サンプリング数」のところにサンプリングするデータの個数を入れる。なお、重複を許すことにより、数が限られた母集団を「分布」 とみなしている。ちなみに、truehist関数はちょっといい感じのヒストグラムを書く関数。MASSパッケージに入っているので、 使う前にはパッケージの読み込みをお忘れなく。

まずはサンプリング数1。平均も何も無い。そのままのデータ。

sanp1

当然ながら、元の分布とほぼ同じ。

それじゃサンプリングを2つにしてみよう。

samp2

真ん中に偏った。これは、2つのデータを無作為に取った場合、 近いところにあるデータが選択される確率よりも互いに離れたところにあるデータが選択される確率が高いことによる。 離れた2つのデータの平均を取れば、当然その値は中心寄りになるという訳である。

では、少し飛ばしてサンプリング数を10に増やしてみる。すると、

samp10

もはやほとんど正規分布だ。試しに標準正規分布を重ねてみよう(単位は無視)

samp10plusnorm

ところで、最初の分布は-5〜10までの間でヒストグラムが描かれていたのに、 10個のサンプリングデータの平均値は0〜6くらいの間に散らばっている。

サンプリング数を一気に増やして1000にしてみよう。

samp1000

形は正規分布だが、幅が2.8〜3.6くらいとかなり狭まっている。 ここでこの分布の平均と分散を計算してみよう。

> mean(z)
[1] 3.228538
> var(z)
[1] 0.01700635

元の分布の平均と分散はこうだった

> mean(x)
[1] 3.230769
> var(x)
[1] 17.30462

平均がほぼ同じなのは見ての通りだし、納得もしやすい。

では、分散はどういう関係にあるのか。比を計算してみよう。

> var(x)/var(z)
[1] 1017.538

サンプリング平均の分散は、元の分散の約1000分の1になっている。 1000という数字はサンプリング数と等しい。つまり、 こういうこと。

平均μ、分散σ2である何らかの分布からサンプリングされたデータの平均値は、平均μ、 分散σ2/(サンプリング数)の正規分布に従う。

これが、中心極限定理。

サンプルが多ければ分散が小さくなり、平均値の良い推定が得られる(大数の法則)ということでもある。 サンプリング数をどんどん増やしていけば、分散は0に近づいていくので、平均値の分布はδ関数のような形へ収束していき、 真の平均値が明らかとなる。

とはいえ、ふつうそのためには莫大な(というか無限大の)サンプルが必要なので、 ある程度のサンプル数から平均値が少なくとも存在するであろう区間を推定する。区間推定に関しては専門書を参考にしてもらうとして、 この話はこれでおしまい。お疲れ様でした(自分へ)。

参考文献

posted by Rion778 at 23:55 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。