2009年04月27日

回帰係数(傾き)の多重比較

*お決まりの「作法」として"回帰係数の多重比較"というモノがあるわけではない,と思います.そのまま鵜呑みにしてえらい人に怒られても知りませんよ.

共分散分析

共分散分析を使うと,共変量でもって調整したうえで分散分析をすることができる.たとえば次のようなデータ.

放牧がその後の種子生産に与える影響を調べたデータで,それぞれ初期の根の直径と種子生産量によって散布図が描いてある(統計学:Rを用いた入門書より.データはココ).

このとき,単純に種子生産量の平均値のみを2つの処理間で比べると,「放牧した方が種子生産量が多くなる」という結果が得られる.

しかし散布図から判断すると,それは放牧区に初期の個体サイズが大きなサンプルが集中していた結果のように見える.たとえば回帰直線を2本引いたら,明らかに放牧区の線は「下」にくるだろう.

それで,2本の回帰直線は違うのかとかそういうことをやるのが共分散分析.

これをやるにはそれぞれの回帰直線が平行という仮定が必要になる.要するにその回帰直線を平均値とし,回帰直線からの偏差を分散として分散分析をするんだから.

なので共分散分析は「平行性の検定(回帰の同質性の検定)」「要因効果の検定(共分散分析)」という2つのステップを踏む.

Rで

Rでは次のように行う.まず平行性の検定.

summary.aov(lm(Fruit~Root*Grazing))

ここでFruitは種子生産量,Rootは根の直径,Grazingはグループを表すベクトル.これにより次のような分散分析表が得られる.

             Df  Sum Sq Mean Sq  F value    Pr(>F)    
Root          1 16795.0 16795.0 359.9681 < 2.2e-16 ***
Grazing       1  5264.4  5264.4 112.8316 1.209e-12 ***
Root:Grazing  1     4.8     4.8   0.1031      0.75    
Residuals    36  1679.6    46.7                 

ここでRoot:Grazingとなっている部分は交互作用を表しており,ここが有意でなければ各回帰直線は平行であるということがいえる.この例ではp値が0.75なので問題ない.

共分散分析は次のように行う.

summary.lm(lm(Fruit~Root + Grazing))

モデル式の書き方がちょっと違う.これにより次のような結果が得られる.

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -127.829      9.664  -13.23 1.35e-15 ***
Root              23.560      1.149   20.51  < 2e-16 ***
GrazingUngrazed   36.103      3.357   10.75 6.11e-13 ***

ここでGrazingUngrazedとなっている部分のEstimateが2つの回帰式がどれだけ離れているかを表しており,ここのp値が有意なら2つの回帰式は有意に離れている(つまり切片が有意に異なる)ということが言える.

(もうちょっと説明をしたいところだけど正直イメージがうまく描けてないのでパスです.青木先生のところなんかをよーく読んだら分かるかもしれません.僕は流し読んだのでわかりません><)

回帰式の多重比較

共分散分析では平行性の検定を行う.で,ここで棄却されなければ次へ進める.

逆に言うと,平行性の検定が棄却されるということは切片はともかくとして2つの回帰式の傾きが有意に異なるということ.

というわけで,群が2つしかない場合にそれらの回帰式の傾きを比べる方法については解決だ.上に書いた方法で平行性を調べたらいい.また,「どこかに差がある」ということが分かればいいんだったらやっぱ同じ方法でいい.

でもやっぱ「どことどこに差があるの?」ってことが気になる.

というわけで,群の中から2つを選び総当たりで平行性の検定を行うことにしよう.そのとき問題になるのは検定の多重性というやつ.

Holmの方法による多重比較

たとえば群が3つあったら対比較は3回やらないといけない.危険率が5%として,すべての群が等しいのにどこかに差がある判定をしてしまう確率は1-0.95^3で0.14くらいになる(実際はちょっと違うらしいけど).

それで,全体として有意水準5%で検定をしたいときに,対比較の回数に応じて1回1回の対比較における有意水準を低く調整してやることで何とかしようという多重比較法がある.ボンフェローニ(Bonferroni)の方法をはじめとする有意水準調整法というやつだ.

今回はホルム(Holm)の方法というやつにする.これはそれぞれの対比較で得られたp値を小さい順に並べて,i番目のp値に対する有意水準をα/(k-i+1)にするという方法.kは対比較の回数,αは全体としてそれ以下におさめたい有意水準.ボンフェローニの方法だとすべての有意水準をα/kに調整してしまうので過剰に保守的になってしまうけど,ホルムの方法ではもうちょっとゆるい感じになっている.ゆるいといっても全体としての有意水準はα以下に保たれる.

で,やることは決まった.対比較を繰り返して,Holmの方法で有意水準を調整するんだ.

Rで

関数にしてみました.

pairwise.covar <- function(
                           data.f,       #データフレーム
                           g.i=1,        #グループ変数の列
                           i.i=2,        #独立変数の列
                           d.i=3         #従属変数の列
                           )
{
  group <- data.f[[g.i]]                  #グループ変数
  g.name <- levels(group)               #水準名
  independent <- data.f[[i.i]]            #従属変数
  dependent <- data.f[[d.i]]              #独立変数
  g.n <- length(levels(group))          #グループ数
  k <- g.n*(g.n-1)/2                    #組み合わせパターン数
  pvals <- numeric(k)                   #p値用ベクトル
  cmb <- combn(g.name,2)                #全組み合わせ
  cmb.name <- character(k)              #組み合わせの名前(出力用)
  sign <- character(k)                  #アスタリスク用
  for(i in 1:k){
    pvals[i] <- round(anova(
                            lm(dependent[group==cmb[1,i] | group==cmb[2,i]]~
                               independent[group==cmb[1,i]|group==cmb[2,i]]*
                               group[group==cmb[1,i]|group==cmb[2,i]]
                               )
                            )$Pr[3],7)          #p値の計算と格納
    cmb.name[i] <- paste(cmb[1,i],":",cmb[2,i]) #組み合わせ名格納
    if(pvals[i] < 0.001){
      sign[i] <- "***"                  #0.1%以下は***
    }else if(pvals[i] < 0.01){
      sign[i] <- "**"                   #1%以下は**
    }else if(pvals[i] < 0.05){
      sign[i] <- "*"                    #5%以下は*
    }else if(pvals[i] < 0.1){
      sign[i] <- "."                    #10%以下は.
    }else{
      sign[i] <- " "                    #それ以外は" "
    }
  }
  o <- order(pvals)
  cbind(cmb.name[o],
        p.value=(pvals*(k-order(pvals)+1))[o],
        sign[o])                        #結果のベクトル
}

アスタリスクとか飾りだけどさ,ほら,見やすさのために…

使い方は

pairwise.covar(データフレーム, 
				グループ変数の列,
				独立変数の列,
				従属変数の列)

列の指定は省略可.デフォルトでそれぞれ1, 2, 3.

使ってみる.まずデータの作成.

## 仮データの作成
set.seed(5)
x <- 1:20
y1 <- x+rnorm(20)
y2 <- x*2+rnorm(20)
y3 <- x+10+rnorm(20)
group <- as.factor(rep(1:3, rep(20,3)))
independent <- rep(x, 3)
dependent <- c(y1, y2, y3)
mydata <- data.frame(group, independent, dependent)

xが独立変数,y1~3が従属変数.プロットするとこうなる.

凡例忘れた...○がy1,△がy2,+がy3.

見て分かるようにy1とy3は切片が異なるけど傾きは同じ.y2だけ傾きが違う.

これに対してさっきのを使うと

> pairwise.covar(mydata)
             p.value         
[1,] "1 : 2" "0"        "***"
[2,] "2 : 3" "0"        "***"
[3,] "1 : 3" "0.639056" " "  

となる.ちゃんと計算できてる?

ちなみに有意水準を減らす代わりにp値を増やしているので場合によってはp値が1を超えるとかいうことになるかもしれない.あとp値が考えている有意水準より高いからって傾きが等しい=平行というのにはちょっと無理があるかも.ボンフェローニよりマシとはいえホルムの方法も保守的なのには違いがないので.

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

上限と下限[位相のこころ p.31〜]

A⊆Qについて,

  1. すべてのx∈Aにたいしてx≦a
  2. すべてのx∈Aにたいしてx≦u ならばa≦u

というようなaは,uをxにいくらでも近づけることができるために,x≦aというようなaの中で最小のもの, Aのすべてを抑えるぎりぎりのものとなる.このとき,これを∪Aと書く.上からギュッと押えるイメージ.これをAの上限と呼ぶ.

先の1.,2.において不等号を逆にしたもの,つまりAを下から押えるぎりぎりのものを∩Aと書き,Aの下限と呼ぶ.

特にQの場合は上限をsupA, 下限をinfAなどとも書く.

上限,下限は必ずしも存在するとは限らない.上限・下限が存在しない例として次のようなものが挙げられる.

  • Q  …  有理数全体を「抑える」 ことはできない
  • ∪φ  …  抑える対象が存在しない
  • A={x∈Q|x2≦2}    …  √2は有理数の中に存在していない(Aの中と外の間には途切れ,gapが存在している)

また,P(E)の場合, AP(E)という部分集合の集まりについて∪Aは全てのX∈Aを含むもの, つまり合併集合となる. そして∩Aというのは全てのX∈Aに含まれるもの, つまり共通集合となる.また、 Eの部分集合をひとつも取らないΦP(E)については普通次のようにする.

  • Φ
  • Φ=E

これは最初の1.,2.に当てはめてみると,1.についてはそもそも判断すべきXが存在しないので自動的に満たされ、 2.についてはUというのが任意のU∈P(E)になるため、「任意のUに含まれるX=φ」 「任意のUを含むX=E」というようなことになる。ただ、なじめなかったら定義と思ってしまってもいいそうだ。

また,QP(E)について次の性質が挙げられる.

  • ∪{x|x<a} = a
  • ∪{x|x>a} = a

ただし,自然数の集合Nでは次のようになって成立しない.

  • ∪{x∈N|x<n+1} = n

関数 f:Λ → Q について, たとえば∪f(Λ),∩f(Λ)というようなものを表したいとき,次のような書き方が使われる.

一般のxnについては同様に次の記法が使われる.

(ちなみにLaTeXで∪のでかいやつを描くコマンドが\bigcupだと調べるのにちょっと苦労した. 考えたらわかりそうな気もしたけど.)

また,A⊆QについてAをQに埋め込む単射 i:A → Qではxとi(x)を同一視して次のような書き方がされる.

また,増大列(x0 ≦ x1 ≦ x2...)について,

は同値である.

なぜなら,b<a<cについて,n≧Nならb≦xn≦cとなるようなNがとれる.いま考えているのは増大列なので, xn≦cはすべてのnで成立する.そこでcの下限を考えると, xn≦aとなる.

----

まだ大丈夫だけどこれは他の本(高木貞二の解析概論とか.ほんの1mmくらいかじっただけだけど.)で見たからだな. どうもやっぱある程度高等な予備知識を要求されている気がする.

このままだとすぐ付いていけなくなりそうなのでなんか別にもう一冊くらい探さないと.などと思って集合・ 位相入門 を頼んだ.ちょっと頑張ろう。

教科書

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

順序構造[位相のこころ p.28〜]

有理数全体: Q
空間Eの部分集合全体: P(E)

順序構造を考える対象

(O1) x ≦ x

(O2) x ≦ y, y ≦ z なら x ≦ z

というような「公理」が成立するもの。Qでいうところの大小関係、 P(E)でいうところの包含関係(≦の代わりに⊆を使うアレ)など。

・「公理」を満たす極端な例

x ≦ y なら x = y …  自分以外と比較できない。比較ができるってことは、そいつが自分っていうこと。 離散順序とでも呼ぼう。

任意のx, yで x ≦ y … どいつとでも比較していいけど、 どっちがでかい(あるいは同じ)ことにしてもいい混沌順序とでも呼ぼう。

こういう変なのは避けたい…

(O1)、(O2)だけでは不十分(擬順序という)。「順序」 というときは次を加える。

(O3) x ≦ y, y ≦ x なら x = y

比較のときにxとyを入れ替えてもいいなら、そいつは「同じ」なんだ、という縛りを入れてやる。 別のもの同士を比較するときは入れ替えちゃだめというわけで、混沌順序を避けられる。

また、昔はこれも加えていた

(TO) 任意の x, y は x ≦ y または y ≦ x

けどこれはなくてもいい。 あるとQは大丈夫だけどP(E)とかが順序の仲間に入れない (二つの有理数はいつでも大小関係を定義できるけど、集合XとYはいつも包含関係にあるとは限らない)。 これを満たすものは全順序と呼んでもうちょっとイイモノだとしておく。

次に離散順序を避けるには…

(DO) 任意の x, y について u ≦ x, y ≦v となる u, v が存在する。

集合から適当に2つの元をとってきたとき、 その両方より小さいuとその両方より大きいvが存在するということ。 xとyが別物だったとき、 その両方と比較できるu, vが存在するということだから、 離散順序みたいに比較できる集まりがポツポツと飛んでいるものは避けられるという仕組み(と思う)。

これは有向順序という。Qみたいな全順序はもちろん有向順序。 P(E)もその元Xについてφ⊆X⊆E(少なくとも空集合は含むし、空間には含まれる) なのでやっぱり有向順序。

ついでにもうちょっと厳しくする。

(LO) 任意の x, y について

a ≦ x, y ≦ b

u ≦ x, y ≦ v なら u ≦ a ≦ b ≦ v

となる a, b が存在する。

a, b は次のように書く。

a = x ∩ y

b = x ∪ y

具体的に言うと、

aはxとyの両方より小さい、または両方に含まれる、 といったようなもののなかで一番デカイやつ

bはxとyの両方より大きい、または両方を含む、 といったようなもののなかで一番ちっこいやつ

aは下から押さえる感じ(だから、「デカイやつ」とか言っても実態は普通小さい。小さい中で一番デカイ)。bは上から押さえる感じ (同じく「ちっこい」といっても実態はでかい)。こういうものを束順序という (ちなみにこの束(laticce)はファイバー束とかの「束(bundle)」とは違う。違う概念に同じ名前。ややこしい!)。

ただの言い直しだけど、

Qの中の話だったらaはxとyのうちで小さいほう、bは大きいほうのこと。

P(E)の中の話だったらaは集合XとYの共通集合(共通部分)、bは合併集合(和集合)。

…というように、「順序」 と一言に言っても向きもへったくれもないものや自分と違うものとは比較を拒むようなものからいろんな段階のものが考えられるんだけど、 適当にx, yと2個持ってきたときに、 それらをまとめて押さえられる何かを定義できるようにしておくと便利ですよと、そんな話(だろうと思う) 。それらをまとめて押さえられる何かを定義するのは次のパートですね。

----

位相とか集合とか一度ちゃんとやらないと、とか張り切って読み始めた。

1時間で3ページしか進んでないとかどういうことなの…。

a ≦ x, y ≦ b とか書いてあるときに「a ≦ x」, 「y ≦ b」と勘違いしてて30分くらいロス。 xとyまとめてるだけなのに。x ≦ y, y ≦ zと定規でもって比べたら、a ≦ x, y ≦ bの方はカンマの後のスペースが0.5ミリくらい短かったよ!

本の最初のほうはその人の記号の使い方とか説明の仕方とかに慣れてないからどうでもいいとこで躓いてしまう。 この辺の話だって他の本で見たはずの内容だけどあらためて丁寧に追ってみると今みたいなつまらない躓きが…。 まあでも慣れたからペース上がるはず。きっと。おそらく。多分…。

ところで(TO)とか(LO)とかいう書き方は何なんだろう。人名とかの略かな?

教科書

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

メンデルvsフィッシャー?

(※2016/8追記)Fisherより指摘されたデータのバイアスは確認されなかったということで決着が付いています(cf. https://en.m.wikipedia.org/wiki/Gregor_Mendel#Controversy However, reproduction of the experiments has demonstrated that there is no real bias towards Mendel's data.以下を参照)。

メンデルの法則といえば高校生物でおなじみの優勢・分離・独立の3つの法則のことだが、これの元になった実験については 「データがあまりに都合のいい数値である」という批判があることは有名だ。

要するに、メンデルの法則が正しいかどうかはさておいて、メンデルの実験には捏造疑惑があるのだ(もちろん、 メンデルの法則がいくつかの例外を除いて、 あるいは例外についても適当な説明方法によって正しいということが主張できるのは周知のとおりだが)。

しかし、日本植物整理学会の「みんなのひろば」という一般の人からの質問に答えるコーナーの質問と回答をまとめた本、 「これでナットク! 植物の謎 」 という本で、次のような記述があった。

しかしメンデルにならって交雑の際に慎重に個体を選んだ三人の研究者が、それぞれ独自に出したデータも、メンデルと同様に 「きれい」なものだったので、フィッシャーの批判はあたらないとされています。

本当か?

ちなみに、Webページはこちら。 編集の都合だろうが書籍よりも記述がいくらか多い。

おそらく、「三人の研究者」というのはメンデルの法則を再発見したド・フリース、コレンス、 チェルマクの3人のことだろう。「工学のためのデータサイエンス入門」 にメンデル、コレンス、チェルマクのデータが載っているので引用する。

実験者  黄色個体  緑色個体
メンデル    6022     2001
コレンス    1394      453
チェルマク   3580      1190

ここでそれぞれのp値は0.9076, 0.6480, 0.9467である。

p値というのは、「理論が正しいとして、この比率よりも外れた比率が観測される確率」を意味する。つまり、 たとえばメンデルのデータから観測された0.9076という値は、「理論が正しい」ときに観測されたとしても問題のない範囲に存在している。

と、ぱっと見はそう思うだろう。

ただしこれを逆に考えると、「理論が正しいとして、この比率よりも『理論に適合する』 比率が観測される確率は1-pである」ということが言える。つまり、「理論」、3:1という比率、 これにメンデルの観測データはきわめて近い。そしてメンデルのデータよりもより「正確な」 データが得られる確率は1-0.9076=0.0924にすぎない。

コレンスはともかく、チェルマクのデータについてはメンデルよりも「正確」であり、 このデータが得られる確率は1-0.9467=0.0533、つまり5%程度である。

さらに悪いことに、メンデルが予想した「誤った理論」に「正確」なデータもメンデルは得てしまっている。 そこではある対立遺伝子についてF1個体の分離比をAA:Aa=201:399と報告していたのだが、 AAとAaの個体の識別方法に問題があった。メンデルは優勢形質が得られた個体からの種子10粒を育てて自家受精し、得られた種子にaa型 (劣勢)形質が見られれば親はAa型であり、全てAA型(優勢)ならば親はAA型であると判断を下した。 しかし確率的には(3/4)^10≒5%で親がAaにもかかわらず子が全てAAになってしまうという場合が生ずる。 よってメンデルの方法によって観測されるべきF1世代の遺伝子型の分離比は223:377でなければならないのである。

それに加えてなお悪いことに、メンデルは他にも数多くの形質について分離比を調査しているのである。そして、 フィッシャーの批判というのはこれらのデータ全てを総合して考察した結果にたいして向けられている。フィッシャーの計算するところによれば、 メンデルが報告した値よりも「正確な」データが得られる確率は、0.00005、 つまり10万回に5回という極めて小さな値となるというものだった(これについては「やさしい統計入門 」 から引用)。

以上から次のことが言える。

まず、メンデルのデータが当てはまりすぎているのは事実である。

つぎに、メンデルに続いた「三人の研究者」のデータが「当てはまりすぎている」かどうかというと、 これは微妙なところだ。ド・フリースのデータについては調べていないので分からないが、とりあえずチェルマクのデータはたしかに「きれい」 に当てはまっている。しかし、コレンスのデータはp=0.6480であり、これは実際の実験で「ありそうな」値だ。決して 「当てはまりすぎている」というものではない。

なので、続く研究者のデータも「きれい」なものだったからフィッシャーの批判はあてはまらない、 というようなことは言えない。

もちろん、植物が確率に従わないような何らかの制御機構を持っており、分離比が確率の示す以上に3: 1に近づくのだ、というような説は否定できない。ただ、「工学のためのデータサイエンス入門」には 「三人の研究者」 以外の研究データもいくつか載っており、 それらのp値だけを挙げると0.7408, 0.3779, 0.1742, 0.4488, 0.0950である。 理論に対して適当なずれを持っており、「きれい」というわけではない。

どうも私の調べた範囲ではフィッシャーに分があるように思えるがどうだろうか。

ちなみに、フィッシャーの意見は 「どのような結果が得られることが望ましいかを熟知している数名の研究補佐員によって、メンデルは欺かれていた可能性がある」というもので、 メンデルを直接に攻撃するものではなかったようだ。

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

デルタ関数の定義とグラフのプロット

ディラックのデルタ関数の基本的な定義と性質は前回見たけど、 数式としての定義は特にしていなかった。

定義の仕方はいろいろある。けども、大体以下のポイントを抑えたものが使われる(っぽい)。

  • -∞〜+∞の区間で積分すると1になる
  • x=0において無限大の値をとる

前回の(1)の定義によるとδ(x)はx=0以外の点で0にならないといけないんだけど、これはどうも無視されることがあるらしい。 というのも、デルタ関数を実際に使うにあたって重要なのは前回でいう性質の2つ目、式(4)のところで説明した「値を抜き出す」 ということであって、これが守られれば少々のことは無視してもおkということらしい。

とりあえず教科書に載ってるのを列挙する。過程とかあまり書いてないので説明は無理っぽ。 ただ少なくとも最初に書いた2つの点は抑えてる。

1個目。

20090224eq01

2個目。

20090224eq02

これは実質的に1個目と同じらしい。Imというのは複素数の虚部を取り出す演算子。iは虚数単位で√-1。

3個目。

20090224eq03

おなじみ正規分布の密度関数を平均=0のもとで分散を0に限りなく近づけた場合と同じ。 これはもともと確率の密度を与える式なんだから全区間積分すると1になるってのは分かりやすい。 分散が0になるとx=0の部分が無限大に近づくというのもイメージしやすい。個人的に一番「それっぽい」気がする。

4個目。

20090224eq04

これが最初に言った定義を「少々無視」する形式。sin axはaをでかくすると0付近でどんどん激しく振動するだけで収束はしない。 ただ、この形をデルタ関数としておくとフーリエ変換の際などに便利らしく、この定義が結構頻繁に使われるぽい。

教科書的には式(4)は三角級数からいくつかのステップを踏んで誘導する感じで紹介されてるけど、 そもそもその三角級数が出てきた経緯がちょっと分からなかったので完成したブツだけ。

この4つについて適当なaの値を与えてグラフをプロットし挙動を観察してみる。

まず式(1)。そんな無茶苦茶早く収束するっていうわけでもない感じ。

Delta1

 

式(2)。プロットすると式(1)と同じものだってことがよく分かる。ちなみにRで複素数の虚部を取り出す関数はIm()。 複素数のオブジェクトを作って虚部にだけ値を与えるにはcomplex(imaginary=1)などとする。 一応ここだけプロットに使った関数のソース。

Delta.2 <- function(a, x){
 a <- complex(imaginary=1/a)
 Im((pi*(x-a))^-1)
 }

Delta2

 

式(3)。収束が早いっぽい。

Delta3

 

そして式(4)。aの値を大きくするとブレがむしろ拡大していくのが分かる。

※「ブレが拡大」という部分ですが、軽い指摘を頂いたのでよくよく考えてみたら確かに微妙な表現だと思えてきました。 x=0の付近で値は激しく振動し、ピークは鋭くなっていくわけですが、周辺での振幅は拡大されるわけではありません。 試しにa=200とした場合のグラフを描いてみました(↓に追加してあります。 もっと大きくてもいいのですがこれ以上だと塗りつぶされてしまいます。)。周辺部の振幅で言えばこれは同じと言うべきでしょう。重要なのは 「収束はしないものの、激しく振動しつつx=0でのピークが高く鋭くなっていく」という点でしょう。ちなみに教科書↓の表現を要約しますと、 「sin axはa→∞としても0となるわけではない。しかしa→∞で激しく振動するため、『値を抜き出す』性質は成立する」。 (2009/02/24 追記)

 

Delta4

Delta5

で今日はおしまい。ええ、手抜きです。

…なんだけど、Rで数式書くのが案外面倒だった。なので数式の部分だけメモっておく。 使うときはプロットのmain=とかの後にそのまま突っ込むだけ。

#式(1)
 expression(lim(over(a, pi(x^2+a^2)), a %->% 0))
#式(2)
 expression(lim(paste(Im, over(1, pi(x-ia))), a %->% 0))
#式(3)
 expression(lim(paste(over(1, paste(a, sqrt(pi))),exp,
 bgroup("(", paste("-",over(x^2, a^2)), ")")), a %->% 0))
#式(4)
 expression(lim(over(paste(sin," ",ax), paste(pi,x)), a%->%infinity))
#ついでにプロット時の余白の設定
 par(mar=c(3,3,3,1), mgp=c(2,1,0))

試行錯誤しながら書いた部分があるので正しくないかもしれない。式(3)とかなんだかpaste()たくさん使ってるし。 expression()の使い方あまり把握しきれてない気がする。

やっぱ数式書くならLaTeXのが楽だ。いちいち書き方調べなくても勘で大体書けるし。 LaTeX使って式を挿入する方法を調べないといけない。 どうもPSfragとかいうのでポストスクリプトの中身を書き換えられるらしい。 それでいけるのかも。また今度。

教科書

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

ディラックのデルタ関数

ディラックのデルタ関数(δ関数)の定義の仕方はいくつかあるらしいけど、たとえば次のように定義できる。

20090222eq01

基本的に0だけど、x=0の点を含む積分では1になる。積分範囲は0から任意の実数ε (0とほとんど同じくらい小さくても、無限大でもいい)だけ離れた範囲。x=0の点で無限大の値を持つ、 というような表現がされることもある。

本来積分を含めてしか性質が現れてこなくて、そういう意味で「超関数」とか呼ばれてる(多分)のだけれど、 無理やりグラフにするとたとえばこんな感じ 。

Delta

εはすっごい小さい数で、とりあえず形は正規分布になってる。正規分布において平均=0、 分散→0の極限をデルタ関数として定義することもできるらしい。

まあとにかく、デルタ関数は普通0を中心に対称な、偶関数として定義される。よって、 偶関数として定義されたなら、という条件付で

20090222eq02

とりあえずここまで定義。以下特徴の説明。デルタ関数が導入されるときによく説明される「性質」をいくつか。

それで、今の図において−ε〜εの範囲で線の下の面積が1なんだから、−ε〜0、 0〜εの範囲に含まれる面積は1/2じゃないといけないだろう。よって、定義から導かれる性質の1つ目。

20090222eq03

これ以上特に説明も不要だろうから次へ。

性質の2つ目。

20090222eq04

任意の関数f(x)とδ(x-a)(この場合デルタ関数はx=aにおいて値を持つ)の積を積分する(区間はx=aを含む)と、 x=aにおけるf(x)の値が得られる、ということを言っている。 これは次のように説明できるだろうと思う(表現として正確じゃないかも)。

20090222eq05

積分範囲をx=aから±εの範囲とする。εは任意の数なので、非常に狭い区間を考えるとf(x)→f(a)。 右辺δ(0)×dxという部分は、リーマン和の一切れだとイメージしてほしい。 縦×横で面積が出てくるということを言っている。縦がδ(0)、横がdx。それで、デルタ関数の下の部分の面積というのは1だったから、 結局f(a)が残る。

次、性質の3つ目。

20090222eq06

「中」をa倍するのは、「外」から1/|a|倍したのと同じこと、ということを言っている。 これは次のように証明する。

まず、δ(y)というものを考えると、定義より∫δ(y) dy = 1。ここで積分変数をyからaxへ変更する。置換積分だ。

dy/dx = aなのでdy=a×dx。よって次式が得られる。

20090222eq07

2段目の=の後は変換ではなくて単に等しいということを言っているだけ。定義から積分の結果はどっちも1になるので。それで、 式(7)の後半2つを使うと次が得られる。

20090222eq08

一段目は式(7)をaで割っただけ。2段目は両辺をdxで微分している。積分して微分したら元に戻るので単純に積分記号をはずすだけ。

しかしaを-aとおいて式(8)まで誘導すると1/aが-1/aとなってしまう。

ここで最初の定義(2)を思い出すと、偶関数なのでδ(ax) = δ(-ax)。よって、aを-aとして計算した場合では、 式(7)の時点で得られる -a∫δ(-ax)dx積分の中身が必ず正となるため、 a < 0である必要性が生じる。逆に言えば、 a < 0であるとき、 式(7)の時点でaの前に-をつけなければならない

よって-aからはじめて式(8)の時点で得られる-1/aという値は非負でなければならない

なので式(6)のように絶対値記号のついた表現となっている。

次。性質その4。

20090222eq09

性質3の拡張なのだけれど、axをいきなりf(x)としているから分かりにくい。aiという定数が右辺に含まれるが、 これは関数f(x)の実根である。つまり、f(ai) = 0となるような値。

ちなみにf(x)は任意の関数というわけには行かない。右辺でゼロ割りが生ずる危険性があるので、 f(x)はf'(a)≠0という制限がつく。これは結構厳しい制限だ。 x軸と平行に交わる関数はすべてだめなのだから。たとえばx^2という単純な関数を考えるとその根は0、 導関数は2xなので、f'(a) = f'(0) = 2*0 = 0となってしまうために式(9)は適用できない。

さて、式(9)は根が複数あるような関数を想定しているが、とりあえずひとつからはじめよう。

20090222eq10_13

まず式(10)は部分積分。yをf(x)とおいて積分変数をxへ変更している。

次、この積分が非0となるのはデルタ関数の中身が0となる点を含む区間においてのみであるから、 f(x)が0となるようなx、 つまり実根aから±εに積分範囲を限る。この操作を表しているのが式(11)。

εは任意の実数なのでいくらでも小さくとることができる。εを十分小さくとったとき、 f'(a±ε)をf'(a)という定数で近似できる。 定数は積分記号の外へ括りだせる、 ということを言っているのが式(12)。

式(13)は式(7)のときと同じで変換ではなく、単に等しいということを言っている。式(11)で積分範囲をa±εに限ったので、 xがそのような区間しか動けなかったとしてもちゃんと値が1になるようにしなければならない。 そのためにδ(x)ではなくδ(x-a)としている。式(13)の積分範囲はa±εより広ければ無限大でもいい。 ちょっとずれてるだけで基本のデルタ関数だ。

それで、式(12)と式(13)から次が得られる。

20090222eq14

そして、実根が複数ある場合を考えると、 どの実根を含む範囲においても式(14)が成立する必要がある。 これを表現するにはすべての場合を単純に足し合わせればよく、したがって次式が得られる。

20090222eq15

そして、式(6)で絶対値記号をつけたのと同様の議論をすることにより、やはり絶対値記号つきの式(9)が得られる。

----

教科書なんか見てると四つ目の性質、もしくはそれに類似する「性質」や「特徴」の説明って結構証明なしに与えられることが多い。 もしくは証明することが課題になってたりして。しかしΣと絶対値記号が突然出てくるのは結構ビビる。混乱する。

そんなのちょっと考えればすぐに解決できる程度のレベルの人間を想定してるのか、 あるいはいい先生に教えてもらうのが普通なのかは分からない。

でも式(9)の形の変換ってグリーン関数とか使う場合に出てくる気がする (伝熱関係の論文で見たことある気がするだけで使ったことはないんだけど)ので、ちゃんと把握している必要がある事項だと思う。

ここに躓いたのいつだっけか。確か半年くらい前。そのときは何が道具として必要なのかの見当がつけられずにあきらめたけど、 結局は置換積分と「εをすっごい小さくする!」という微積分によくある操作だった。

やっぱ基礎だな。基礎を徹底的に固めないといけない。誰か固めるべき基礎をリストアップしてくれ。

----

参考文献

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

最小作用の原理〜ニュートンの運動方程式

力学系は次の積分を最小にするように運動する。

20090218eq01

というのが最小作用の原理。もしくはハミルトンの原理。積分Sは作用と呼ばれる。Lは一般化座標q, その導関数q'(一般化速度) , 時間tからなる系を特徴付ける関数で、ラグランジアンという。(「特徴付ける」って何なのさ。という感じだけど、とりあえずここは 「ラグランジアンとして具体的な形を考えなくてもしばらく計算が進められる」 ということを強調するためにこんな微妙な言い方をしているだけなので、細かいことは考えず先へまいりましょう。)

一般化座標というのは系において位置を決めるのに必要な任意の個数の量の組で、たとえば平面上で垂直方向へy、 水平方向へxと座標を設定すれば、一般化座標の一例は(x, y)となる。「一例」 といったのはたとえば水平面からの仰角θと直線距離rを用いて(θ, r)としても位置は一意に決定できるからだ。 ちなみに位置を一意に決定するのに必要な量の個数をその系の自由度と呼ぶ。

さて、積分Sが最小となるということの意味を数式で表してみよう。変分記号δを使えば、それは次のようになる。

20090218eq02

ここから前に変分法を説明したときと同様の手順を踏んでオイラー方程式を得る。

20090218eq03

このオイラー方程式はラグランジュ方程式と呼ばれる。

次にラグランジアンの具体的な形を示す。重力のある質点系においてラグランジアンは次の形で定義される。

20090218eq04

Tは運動ポテンシャル。Uは位置ポテンシャルを表す。

いきなり出てきた感があるけど、これは相互作用のない系におけるラグランジアンに、相互作用(質点ともうひとつの質点、 地球との相互作用である)によって決まる座標の関数Uを足し合わせることによって得られている。

運動ポテンシャルは力(F = ma = m dv/dt)を位置について積分することで得られる。 物理の教科書とかによく載っているので多少飛ばし気味で(ただしランダウ=リフシッツ「力学・場の理論」では以下と異なり、 先に述べたように「相互作用のない系におけるラグランジアン」として運動ポテンシャルを導いており、議論は最小作用の原理だけを用いて進む。 が、ちょっとまだ理解しきれてない部分があるのでそれはパス。気になる人は「力学・場の理論」でどうぞ!)。

20090218eq05

一列目から二列目の変換は置換積分であることに注意。 ここで積分変数がqからvへ変わっている。

位置ポテンシャルも同じノリで重力(F = mg)を高さhについて積分することで得られる。

20090218eq06

こっちは置換積分とかしない。

式(5)、(6)から重力のある質点系での具体的なラグランジアンの形を求めれば次の形になる。

20090218eq07

これをラグランジュ方程式(3)へ代入すると、次の式が得られる。

20090218eq08

ところでここで左辺第一項はma(aは加速度)である。第二項、計算前の形で∂U/∂qであるこれは力である。

そう、要するにこれはma = Fというニュートンの運動方程式に他ならない (別にニュートンの運動方程式を求めるだけならUとして具体的に位置エネルギーというものを挙げず、 座標と相互作用により決まるポテンシャルエネルギーとしておけばよかったんだけど抽象的すぎて分かりにくいので)。

たとえばここで垂直方向へy、水平方向へxとして座標を設定し、式(8)へ座標の情報を取り込むと次の式が得られる(空気抵抗は無視) 。

20090218eq0910

これはよくある(と思う)空気抵抗を無視した状態で物体を投射したときの運動を求める問題。まあ放物線が出てくるわけですが、 結構面倒なんで今日は解かずにここで終わり。

とりあえず、最小作用の原理からニュートンの運動方程式が出てくる、というお話でした。

----

参考文献

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

置換積分

関数 f(x) において x=x(t) などとおき、tに関する積分に変形することで積分計算が単純化される場合がある。

この、変数変換ののちにtに関して積分する、という作業を置換積分と呼ぶ。

まず、tという新しい変数がいかにして導入されるかを確認するため置換積分の公式を誘導する。

F(x)をf(x)の原始関数とする。

20090212eq01

このとき、xを別の変数tの関数とみなし、x=x(t)とする。

20090212eq02

こうして導入された新変数tにて式の両辺を微分する。合成関数の微分なので、「外側の微分×内側の微分」となる。

20090212eq03

右辺は積分の微分であり、したがって結果は次のようになる。

20090212eq04

この両辺を先ほどとは逆に新変数tにて積分する。これにより左辺は微分の形から原始関数へ戻る。

20090212eq05

ここで得られた右辺が置換積分という作業となる。被積分関数中のxをtの関数x(t)に変形し、 積分変数をxよりtに変換するためには、被積分関数に関数x(t)をtで微分したものを掛け合わせればよい、ということを示している。

この説明だけだと具体的な計算手順をイメージすることが困難な場合がある(昨晩深夜の僕とか)ので、具体例で確認する。

まず、例として被積分関数f(x)を次のものとする。

20090212eq06

ここで、たとえば次のようにおいてtを積分する形としたら計算が簡単になりそうだ。

20090212eq07

ただ、これだとx(t)というかt(x)だし、f(x(t))というよりf(t(x))のような気がする。 式(5)の使い方がどうもイメージしにくい。

しかし、実際には式(7)を代入するわけではない。tはまだ被積分関数の中にないのだから、代入はxに対して行われなければならない。 当然だ。

つまり、代入されるのは式(7)の逆関数である次の式となる。

20090212eq08

これを式(6)へ代入すれば、(2x+1)がtに置き換わる。式(8)、 つまりx(t)を代入するのだからf(x)がf(x(t))となることにもはや疑問はない。

再び式(5)に戻る。

20090212eq05

あとはdx/dtを求めればいいだけ。式(8)の両辺をdtで微分する。

 20090212eq09

そしてx→x(t)への変換とともにこれを式(5)へ代入し、あとは計算して答えを得る。最後にtをxに戻すのを忘れずに(※ 積分定数は省略)。

20090212eq10

----

それで、とりあえずこれで正解は正解なんだけど、実際にはt(x)の逆関数のx(t)をわざわざ求めて、 それを微分したものを代入して…などという手順は踏まないし、そんな回りくどい説明はしない(だからわからなかった!)。

まず、次の式は自明だろう。言い方が正しくないかもしれないが、左辺でdtを約分したらdxになるっぽいのは間違いないだろう。

20090212eq11

この式の左辺はdxに等しい上、これで積分したら積分変数はtである。要するにこの式の左辺を求め、与式のdxへ代入したら、 もうxをx(t)にしてしまってかまわないということ。

実際に何をするのかを先の例を使って説明する。まず置換は式(7)と同じく次のようにしたいものとする。

20090212eq12

ここで両辺をxにて微分する。

20090212eq13

dxについて解く。

20090212eq14

これが要するに式(11)。なので、式(12)の変換をしたあと、式(14)を与式へ代入したらあら不思議、 もう積分変数はtになっていると、そういうわけです。

あとの計算は一緒。もちろん答えも一緒。なので省略。

----

ところで、置換積分にあたって何を仮定しないとだめなのか。

なんか逆関数とか面倒なこと言ってたから、x(t)は単調じゃないとだめなんじゃないかなとかそんな風な気もする。

けど別にそういう仮定でなくともいいらしい。

最後に仮定を持ってくるのもどうかと思うけど、「解析概論」から丸写しすると置換積分の公式を得るに当たっての仮定は次の二つ。

  1. 積分区間 a ≦ x ≦ b を含む区間 c ≦ x ≦ d においてf(x)は連続。
  2. x(t)およびx'(t)は[α, β]で連続で、tがαからβまで変動するとき c ≦ x(t) ≦ d, かつ x(α) = a, x(β) = b。

多少ハミ出ても連続で着地点が一緒ならいいと、そういうことですね。

…多分。

参考

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

変分法 -デカルト座標における2点の最短経路を求める練習問題

前回変分法という方法をつかってオイラー方程式を求めた.

これを使って実際に問題を解いてみる.前回の冒頭の問題はこれ.

  >原点O(0,0)を出発し、点P(x, y)に至る最短の経路は何だろう。

x, yだと後々わかりにくくなるので,点Pは今回(a, b)とする.

最短の経路というのは,次の積分Iを最小とするようなものだった.

とりあえずこのままだと解けないので積分変数をxにする.xにするにはまずdsをdy, dxの関係から次のように変形する.おなじみ三平方の定理.

そうしたらこれを(2)へ代入する.積分変数はxになり,積分範囲は0〜aになる.

ここでの被積分関数がオイラー方程式におけるfなのだけれど,よく見るとこの関数にはyの値が関与していない.xとy'さえ与えればいいらしい.

ということは,オイラー方程式における第二項が0であって,次のように変形できる.

xで微分して0なのだから,関数fをy'で偏微分した値は定数じゃないといけない.また実際に関数fをy'で偏微分してみると,次の結果が得られる.

つまり,y' = dy/dx, xの増加量とyの増加量の比は常に一定でなければならない.

で,まあ残りは解くというようなレベルでもない微分方程式を解くだけなのでぱぱっと.

というわけで,原点O(0, 0)と点P(a, b)を結ぶ最短の経路は,y = (b/a)x という一次方程式でした.デカルト平面で最短経路をたどりたかったらまっすぐ進めばいいわけですね!!

----

…しかしこれだけだとあまりに簡単すぎる上,わざわざ計算する有難みがない.なのでもう一問くらいは練習する.

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

変分法 ―変分問題からオイラー方程式まで

原点O(0,0)を出発し、点P(x, y)に至る最短の経路は何だろう。

…まあ普通に考えて直線なんだけど、直線かどうかはひとまずおいておく(たとえば座標の中に「動きにくい場所」 とかあったらそこをある程度避けるほうが効率的になる。このときは「最短経路」は時間にかかわったものになる)。

目的とするのはxとyを結びつけるような関数y(x)を求めることだ。

とりあえず経路を微小素片dsに分割する。そうすると、最短経路というのはこの微小素片dsの合計、 すなわち積分が最小となるようなものだ。dsの積分、つまり経路長をIとし、数式で表現すると次のようになる。

20090208eq01

ここでdsの進み方はわからないとして(もちろんy=axという一次式なんだけど)、たとえばx, y, y'(=dy/dx)という3つの値で決定されるf(x, y, y')という関数で表現できるものとしてみる。 なぜyについて1階の微分までしか含めないかというと、座標と速度さえ与えればその場に働く力などによってより高次の微分、 つまり加速度やその微分量などは自動的に決定してしまうからだ。座標だけでは不十分。なぜなら、 同じ座標でも質点はさまざまの速度を持つことができるから。

まあ細かいことはおいておいて、今のf(x, y, y')を利用して式(1)を変形する。式がx軸方向へdx進むとすれば、 その間に微小素片dsはf(x, y, y')dxだけ進む。これにより積分変数をsからxに置き換えて次式を得る。

20090208eq03

繰り返しになるけど、この積分を最小にするような関数y(x)を求めるのが目的だ。これは変分問題と呼ばれる。

ここで、そのような関数y0がわかってしまったとする。そしてy0に「小さな関数」 δyを足し、OとPを結ぶ新しい関数yを作る。この「小さな関数」は関数yの「変分」という。

20090208eq04

δyというのはyにδという微小量がかかっているわけではないことに注意。δyという2文字で1つの関数を表現し、 これは関数y0からのy軸方向の距離をあらわしており、その形は任意である。 (3)におけるy式はy0式と同様に点Oと点Pを通らなければならないので、 式δyは次の境界条件を満たしている必要がある。

20090208eq05

ここでyにδyを足したときの増加量をδIとし、次式を得る。

20090208eq06

このとき、2つの積分の差はy方向への微小な増分、y'方向への微小な増分それぞれの多項式だが、 もしも変分δyが十分に小さな関数ならば、それはyとy'についての一次式で近似される(cf. 全微分)。そして合計の増分は両者の和であるので次式を得る。

20090208eq07

ここで右辺の第二項のみ部分積分する。δy' = d(δy)/dxであることを利用する。

20090208eq08

このとき式(4)の条件から右辺第一項は0になる。

積分Iが最小になるということはδyをどのような方向への関数ととっても積分Iが必ず増加するという意味であり (通常の微分における極小値と同様)、これは積分Iの増加量δIに対する式yの増加量δyの比が0となるということに等しい。

20090208eq11

このとき右辺の分母は0でなければならないが、δyが任意であるため、それにかかる括弧の中身がいたるところで0となる必要がある。

20090208eq12

こうして微分方程式が得られた。Iを最小にする関数を求めることは、この微分方程式を解くということとなる。

変分法における式(9)はオイラー方程式と呼ばれる。

話を物理に移すと、力学系において運動は時刻と位置、加速度からなるある関数f(t, x, x')をt=0からt=tまで積分した際にそれを最小とするように起こる。「ある関数」というのは単に経路だったり、仕事量だったり、 屈折率だったり、場合によりけりだがこれの積分をまとめて「作用」とよび、作用が最小になるというこの原理を最小作用の原理 (またはハミルトンの原理)と呼ぶ。そして最小作用における関数f(t, x, x')は普通大文字のLで表し、ラグランジアンと呼ぶ。 またオイラー方程式はラグランジュ方程式と呼ばれる。

----

変分はいわば「関数の関数」の微分。理解のポイントは微分、偏微分、全微分、 部分積分といったモノをしっかり把握してるかどうかですね。

そのあたりがあいまいなまま「力学・場の理論」の文庫を買って最初の数ページから数年積読になってた。何事も基礎が大事。 独学だったんでどこが基礎なのか分からなかったというのもあるんだけど。

なんか思ったより疲れてしまったのでオイラー方程式を使って最初の問題を解くのは次の機会に。

参考

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

卒論執筆でお世話になったWebサイト・本・ソフトウェア

世間では卒論などが追い込みの時期だろう。昨年の今頃は私もそんな中に。

ここで昨年を思い出しつつ、論文執筆の上でお世話になったリソースをまとめてみようかと思う (ちなみに私は分類で行くと理系だと思います)。

以下、自分が世話になったもののまとめなので、万人に役立つものではないだろうし、 また現在進行形で追い込んでいる学生に役立つものとも限らない。どっちかというとB1〜3、M1などが役に立つんじゃないだろうか。

Webサイト・本・ソフトウェア、それぞれ3つずつ。

Webサイト

東大で学んだ卒論の書き方論文の書き方

もはや説明するまでもない。年に一度、 秋から冬にかけてほぼ確実にはてなブックマークのホットエントリーに浮上してくる風物詩的ページ。

まずはここを読もう。いつでもここに戻ってみよう。

GoogleGoogle Scholar

論文・資料の検索に。

…は当然として、私は英語の表現方法がどうしても分からないときの答え合わせとして使った。 もちろん辞書や例文集は第一に頼るべきだが、その語句、節で検索することで、本当にそのような表現が論文において一般的であるか、 また自分の分野で使われているのかを確認することができる。また、思わぬ論文との出会いのきっかけにもなり、視野を広げるのにも役立つ。

Webcat Plus

キーワードだけでなく、文章でも文献検索ができる。「自分の大学の図書館にあの文献はあるだろうか?ないとしたら、どの大学、 図書館にお願いしたらいいんだろう?」そんなときに。

理科系の作文技術

もはや古典の感は否めないが、それでも論文執筆に当たって第一に読むべき本はこれ。冒頭で紹介した「東大で学んだ〜」と同じく、 執筆の前に読むべきで、執筆に詰まったら読むべきで、見直すときにも読むべき本。

基本的な部分はこの一冊がほとんど教えてくれる。

推計学のすすめ―決定と計画の科学

データ処理にあたって統計学が理解できていないと困るが、統計学はそれほど容易なものでもない。基本的な道具ですら、 少し突っ込むと手に負えない複雑さと難解さを垣間見せる。

そういったあまり見たくない一面を避けつつ、かつ雰囲気として理解しているべき最低限の部分を抑えているのがこの一冊。 統計学が推計学と呼ばれていたころの本だが、基礎部分が対象なので今でも通用する内容。

もちろん、統計学に自信のある人は特に読む必要はない。

卒論執筆に取り掛かってからというより、実験計画程度以前に読んでおくべき本。

統計でウソをつく法―数式を使わない統計学入門

データとはどういう形をしているもので、統計学をデータにどのように適用するか、また、 他人のデータをどのように読み取るべきかということを教えてくれる。 どちらかというとデータの取り扱いに関する哲学みたいなものを扱っていて、具体的な方法が詳細に述べられているわけではない。 それでも一度は目を通しておくべき。長く読まれているものにはやはりそれなりの理由がある。

半世紀もの昔からのロングセラーで、いまだにこの本に書かれているような「ウソをつく法」が使われているのかと思いたくなるが、 残念なことにいまだに使われている場合が多い。

これもできたら実験計画以前に読んでおくべき本。

ソフトウェア

(基本的にWindowsユーザー向けです。 MacとかLinuxとか使う人は大概が変態なので私よりずっとパソコン関係には詳しいだろうし、そもそも勝手にやるでしょう。)

TeX・ LaTeX

(リンク先はTeXに関する有用な情報源であるTeX Wiki)組版ソフト。特に数式の出力が綺麗。 TeXではHTMLのようなマークアップ言語を処理することで、 DVIやPDFといった形式の文書を作成できる。仕上がりは非常に綺麗で、 出版物に使用される例も多い。

学会によっては論文提出用にTeXのスタイルファイルを配布している場合もあり、覚えておいて損はない。 特に数式を使用する可能性のある人は覚えておくべき。Wordをメインに使う場合でも、 数式の可読性を高めるためにTeXで処理した数式を貼り付けると良いだろう。

綺麗な数式が書けるということ以外の利点は、文章の本体がテキストファイルとなることで、 複数のバックアップを用意に取れるということ、突然のエラーでデータが消失する危険性が減るということだろう。また、 完全に無料であるという点も重要。

欠点は導入と初期設定のわずらわしさにある。 ただし最近は必要なものをまとめてインストールしてくれるインストーラもあるのでそれほどでは無いように思う。また、 HTMLなどに慣れていない場合は覚えたりなれたりするのに多少の時間が必要かもしれない。

それとWeb上にはそれ単体で学習できるようないいリソースがあまり無いような気がする。私は[改訂第4版] LaTeX2ε美文書作成入門 という本を買って勉強した。後に辞書的にも使えるのでいい本だと思うが、 他のLaTeX関係の本を持っていないので比較はできない。まあでも一冊で済んでいるということはいい本なんだろう。

xyzzyというテキストエディタにKaTeXというモノを組み合わせて使うのがお勧め(参考: xyzzy + KaTeX [物理のかぎしっぽ])。

JabRef

(リンク先はJabRefの使用方法などを解説した「JabRefによるBibTeX文献管理とJab2HTML」) LaTeXではBibTeXというツールで文献管理ができる。 BibTeX形式で記述された文献リストのテキストファイルを所定の方法で処理することで、 LaTeX文書の中に自動で引用番号を振ったりできる。

また、BibTeX形式の文献情報提供がサポートされている場合も多い(たとえばGoogle Scholarなんかがそう)ため、 LaTeXとあわせて覚えておいて損はない。

JabRefはこのBibTeX形式の文献リストを管理するためのソフトウェアだが、検索機能はもちろんのこと、 論文ウェブサイトへのリンクを記憶したり、論文PDFをダウンロード・保存したりと文献の管理に有効な機能が多く含まれている。

そのため、たとえばBibTeXを使う気がなくとも文献管理ソフトとしてかなり有効なものである。フリーなのも重要。

R

(リンク先はRjpwiki)統計処理とグラフィックスに特化した言語、環境。

CUIなのでエクセルよりはとっつきにくいが、インタプリタ言語で入力内容を逐一解釈、実行してくれるので、 コンパイルが必要なC言語などに比べたらとっつきやすいだろう。 WindowsならRguiというRの作業を行うのに便利な環境とエディタが標準で付いてきて、 完全なCUIというわけでもなくなるのでとっつきやすさはさらに高い。

Rは多くの統計手法が「関数」の形で最初から用意されているのが特徴。t検定、分散分析、多重比較みたいな簡単なものから、 ブートストラップ、クラスタ分析など本当に数え切れないほどの手法がカバーされている。 テキストマイニングやマイクロアレイの解析だってできるし、画像も扱える。汎用性は極めて高い。

Excelとの連携性を高めた関数も多いので、別にExcelを窓から投げ捨てる必要もない。

最初にベクトルやデータフレームを丁寧に勉強し、次にExcelデータからそれらを作り出す方法を理解し、 その後に統計処理を行う関数の使い方やプログラミングへ進むと、統計に関しては比較的短期間で実用的なレベルの能力が身に付くと思う。

グラフは頑張れば綺麗なものも描けるけど、結構な勉強と慣れが必要。 ExcelやCalcの方が手っ取り早いし少ない労力で綺麗になると思う。

情報源として本を頼るのもいいけど、基本的なことならR-tipsでほとんど間に合う。 ありがたいことです。あと青木先生のWebサイトも同時にチェックしておきましょう。

フリー。驚くべきことにフリー。

----

これらは特に印象に残っているし、ほとんどは今でもお世話になっている。

あとはカフェインとブドウ糖とBGM。中毒にならない程度にほどほどに。

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

周期関数と、それの微分の内積

1.フーリエ級数展開

周期関数f(x)がディリクレの条件というゆるーい条件(現実の世界にあらわれるほとんどの関数はこの条件を満たす) を満たしているなら、 その周期関数の整数倍の周波数を持ったsine波とcosine波を足し合わせることによってf(x)をあらわすことができる。

簡単のために周期を2πとすると、

20071224eq9

もしも周期が2πでない場合には、[2π/周期]を掛けることで周期を2πになおすので、基本的にはこれと変わらない。 この展開をフーリエ級数展開と呼ぶ。ちょうど一年くらい前に勉強したので詳しくはそっちの記事で。

2.関数の直交

周期関数は内積を積分によって表現する。2つの周期関数、f(x)とg(x)があったとして、 これらの内積は両者の積を一周期分積分することで求められる。

081228eq1

もしもこの積分の値が0だったとき、2つの関数は直交している、と表現する。直交する周期関数の組み合わせとしては、 たとえばsin(x)とcos(x)、sin(x)とsin(2x)などが挙げられる。

3.関数と、それを微分したものの内積

ここで問題を。

周期関数f(x)と、それを微分したf'(x)は直交している。

これは正しいだろうか?

直感的には正しいと思える。sineとcosineは符号程度の違いをのぞいて微分によって互いに入れ替わる、 つまり直交する状態になるんだから。

まず、f'(x)をキチンと書こう。f(x)としては最初に書いたこれ

20071224eq9

を使う。これの微分は別に難しいことはない。項別に微分すればいい。

081228eq2

そうして内積を考えるのだが、f(x)×f'(x)に含まれる項の組み合わせを真面目に考えるのは厳しい。

まず、f'(x)の積分について考えると、これは[a0/2]という項を含まないので、 平均値は0であり一周期分の積分は明らかに0となる。これを定数倍してもやはり0なので、 まずf(x)の[a0/2]は無視していいことがわかる。

次に、級数部分どうしの積の部分を考える。この積の結果無限個のsin×sin、cos×cos、sin×cosが生ずるわけだが、 その中で一周期分の積分をしても値が0にならない組み合わせは限られる。要するに互いに直交していないものだけが残るわけだが、 具体的にはsin×sin、cos×cosのうちで周波数の等しいものだけ、ということになる。

で、その様なものだけを残した形で内積の計算結果を書いてみる。

081228eq3

うまいこと打ち消しあって見事に0になった。

というわけで、周期関数f(x)と、それを微分したf'(x)は直交している。

3.部分積分

ところで、今のやりかたはとても分かり難い。ちゃんと書きながら読まないと自分でも分からなくなりかねない。 実はこれは部分積分を使うと簡単に解けてしまう。

まず、部分積分が何か、ということの詳細はググってもらうとして、結論だけ書くと次のような積分手法のこと。

081228eq4

積の微分公式をちょいちょいといじって積分しただけのものだが、先の問題にそのまんま当てはめられることが分かるだろう。 当てはめてみると、

081228eq5

なんとまあ簡単な。必要なのはf(0)=f(2π)という条件だけ。周期関数なので始点が不連続点とかじゃない限りこの条件は大丈夫。

なんかちょっと今まで使いどころがさっぱり分からなかったけど部分積分って便利!

---

実のところ考えていたのはf(x)f(x)f'(x)みたいな3つの関数の積の一周期分の積分で、 これが0になってくれないと研究の都合上大変困ることになるところだった。

これを2.でやったように力ずくで解こうとしていた。そうすると、三角関数3つの積がどういう条件で0になるのか、 ならないのかを全部書いて、それらがちゃんと打ち消しあうのかを確かめないといけないわけなんだけど、しかし肝心の 「0にならない組み合わせ」がざっと無限の無限乗個ほど出てきてもうこりゃだめだと半ばあきらめていた。

で、いいわけとごまかし方を半日くらいかけて考えていたときにふと部分積分を思い出したんだ。 まあ半日くらいつっても最近ひどくぼんやりしている時間が多いので実質2時間くらいなんだろうけど。まあとにかく、 今まで単なる呪文だった部分積分の公式の意味と使いどころがようやく分かってきた。

なんというか、練習問題をたくさんこなす、といったような作業が足りないんじゃないか。

知っていても、慣れていないものは使えない。

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

微分方程式とロジスティック式(2)

前回に引き続いてネズミの問題。 もう一回引用すると、

代数と三角関数は知っていたので、どんな数学の問題が出てきても解ける自信があった。 でもシュワブ先生は次のような問題を出してぼくをへこませようとしたんだ。あれから五〇年たった今でも、はっきりと思い出せる。 ネズミがネズミの人口に比例して死ぬとしたら、ネズミの人口はどれくらいかを求める問題だった。
ポール・ ホフマン「放浪の天才数学者エルデシュ」  p.164)

「ネズミがネズミの人口に比例して死ぬ」という一文の意味はよく考えないといけない。まず、

  • ネズミはネズミの人口に比例して増える

ということは書いてないがまあ当たり前だろうということで考慮する。

で、次に「人口に比例して死ぬ」というのをバカ正直に捕らえると前回解いた微分方程式と大差ない問題になる。 次のように解釈すれば問題が面白くなるし、より現実的でもある。

  • ネズミが死ぬ割合はネズミの人口に比例して増える

こう言ってもいい。「ネズミの増殖スピードはネズミの個体数によってブレーキがかけられる」

これを数式で表現すると次のようになる。

081117eq1

前回と同じくxは個体数、tは経過時間、そして増加速度にかかわる比例定数aに対し、xに比例したブレーキがかかる。

これを変数分離によって解こう。

まずは変数の分離。

081117eq2

左辺はそのままだと積分しにくいので、部分分数分解という作業をする。何をするのかというと、 分母が積の形になっていることを利用して、左辺の分数を和の形に変形する。 (この場合はちょっと試行錯誤したほうが楽な気もするけど手順もあるので気になる人は「部分分数」などのキーワードで検索を。)

081117eq3

もうちょい整理して積分。左辺は和の形なので項別の積分ということになる。

081117eq4

ちょっと難しいのは左辺第一項だろう。ほかは前回と同じなので。

∫f'(x)/f(x) dx = ln |f(x)| という式を知っていれば (1/xの積分と置換積分という方法を組み合わせれば簡単に証明できるので、気になって気持ち悪いという人は「1/x」「置換積分」 などのキーワードで検索を)、これは次のように計算できる。

081117eq5

ここで、変域が負の数に及んでしまうことに目をつぶれば絶対値記号を括弧に変えてしまえる。…あれ?いいんだろうか。

とにかく、そんな感じで積分を計算してしまうと、

081117eq6

対数の差は真数の商なので

081117eq7

で、対数を外す(x=a/bは定義できない…?結果から見ると定義できなくても問題ないけど)。Cの中身はちょいと変わる。

081117eq8

ここまで来たらあとは簡単な変形なのでちょっと省略して、

081117eq9

t=0でx=x0として求め たCを代入すると、

081117eq10

というわけで問題が解けた。

ここで解いた微分方程式はロジスティック式と呼ばれ、個体群増加のモデルなどに用いられる。

ロジスティック式は初期値がa/bよりも小さいときにロジスティック曲線と呼ばれる特徴的な曲線を描く。ためしにRで描いてみよう。

logistic <- function(t){
 a <- 1.5 #増加率
 b <- 0.3 #抑制力
 x0 <- 2 #初期個体数
 1/(b/a - (b/a - 1/x0)*exp(-a*t))
 }
curve(logistic, -5, 5, xlab="Time", ylab="logistic(Time)")
grid()

logistic

個体数は0よりも小さい値にはならず、またa/bよりも大きい値にはならない(a/bより大きな値を初期値に与えた場合は…)。 a/bは環境収容力といわれる。

うーん…

答えは間違ってないはずだけどどうも気になる点が。対数の変域が負に行くのはどう説明したらいいんだ?

それに気晴らしのために考えてたのに別に晴れてないしどうしたもんか。

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

微分方程式とロジスティック式(1)

代数と三角関数は知っていたので、どんな数学の問題が出てきても解ける自信があった。 でもシュワブ先生は次のような問題を出してぼくをへこませようとしたんだ。あれから五〇年たった今でも、はっきりと思い出せる。 ネズミがネズミの人口に比例して死ぬとしたら、ネズミの人口はどれくらいかを求める問題だった。
ポール・ ホフマン「放浪の天才数学者エルデシュ」  p.164)

ポール・エルデシュの友人であるロン・グラハムが数学に興味をもつきっかけとなったのがこの問題だ。ちょっと解いてみよう。

が、その前にもうちょっと単純な問題を考えよう。「ネズミがネズミの人口に比例して増える」としたら、 ある時間後にネズミは何匹になるだろうか。(「微分方程式」の意味が分かる人は以下飛ばして次へ行ったほうがいい気がします。 行くべきです。下は飛ばしてー見ないでー。)

まず、ネズミの増殖速度と個体数の比例関係を数式で表そう。

081116eq1

ここでxは個体数、時間はt。個体数の変化を時間で微分したものが増殖速度、 そしてそれが適当な比例定数aによって個体数xと結び付けられている。これは微分を含む方程式なので微分方程式と呼ぶ。 これは特に変数分離形と呼ばれるものの最も簡単な例になっていて、適当な作業によって、xについての式に簡単に書きなおせる。

まずは変数を分離する。変数はxとtだ。これを左右に分別する作業を変数分離という。

081116eq2

dxやらdtやらも気にせず移動してしまう。

そして両辺を(不定)積分する。変数を整理しているので、左辺はxについての積分、右辺はtについての積分を考えることになる。

081116eq2_1

積分という作業がどういうものかを分かっていないとこのインテグラルの唐突感が気持ち悪い(僕がそうだった)が、 まあとにかく両辺にうにょっとインテグラルを書いてしまっておk。で、普通に積分記号からdx、dtまでの間にある式を不定積分する。

081116eq3

ここでCは積分定数で、右の積分で付くやつも左の積分で付くやつもまとめている。もうちょい整理する。

081116eq4

右辺ちょっとおかしいじゃないかと思うかもしれないが(僕がそうだった)、「積分定数」というのが 「面倒なものをまとめただけの入れ物」だと知っていれば納得できる。定数に定数掛けても定数なので、 まあ要するにaCを改めてCと書き直したということ。

左辺にln(自然対数。logと書く人もいるけど、lnと書いたほうが自然対数であることがはっきりする。)が付いている。 lnの意味を考えれば分かることだが、こいつは右辺をeの指数部にしてやると外れる。

081116eq5

また右辺がちょっとおかしいじゃないかと思うかもしれないが(僕が(ry )、ここでのCというのはe^Cを書き直したものだ。 e^C*e^at=e^(at+C)というわけだ。定数の定数乗はやっぱり定数なのでこんな書き換えができる。

これで微分方程式は解けた。のだが、積分定数Cが残っているのはちょっと気持ちが悪い。よくこれは「初期値によって決まる」 と説明される。

たとえば、t=0のときにx=x0などとして代入する。これは容易に

081116eq6

だと分かる。これを代入すれば

081116eq7

となる。

これで問題は終了、解決した。初期のネズミ個体数と、経過時間、それに定数a(これは最初の式を見れば分かるが、 ある個体数のときの増殖スピードをそのときの個体数で割ったものだ)さえ分かれば、一定時間経過後のネズミ個体数が計算できることになった。

ちょっと長引いたので一旦切って次の記事へ引き続き。 以下は余談みたいなもの。

で、この説明は、

  1. 微分記号はバラしていい
  2. バラした記号を使って積分していい
  3. 積分定数に定数のみを含むどんな演算をしても積分定数として1文字で表現できる

ということが分かっていればすんなり(1/xの積分とかは突っ込むと多少難しいが)理解できる。しかし、 これらは明示的に教えてもらわないと結構分かりにくいポイントではないかと思う。 対数あたりから数学は独学なのでこのあたりでかなりつまづいた覚えがある。まあ僕が鈍いだけかもしれないけど。

それで、この問題の答えが次の一般的な事実につながる。

  1. 「指数関数」と単に言った時にeを底とする関数を指すことが多い
  2. 「複利的な増加」「指数関数的増加」を表現するとき、「初期値×指数関数」の形のモデルを用いる

一度でもこの問題を解いたことがあれば本当になんでもないことなんだけど、無いと特に後者が気持ち悪い。 何せ突然eを含む指数関数が出てくるんだから。

それでまー言いたいことは何かというと、自然科学を扱う学課を持つ大学では解析学を必須にすべきだし、 微分方程式について簡単な例から教えるべきだということ。 いま高校では微分方程式を教えていないということを知っている大学教員がどれだけいるのだろうか。 現状の問題がどこにあるのかきちんと把握できているんだろうか、改善に向けた努力はしているんだろうか。

…愚痴を言っていても気分は晴れないので次へ。

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

最小二乗法

なにか2つの数値が一組となったようなデータがあったとしよう(体重と身長、温度と湿度、国語の成績と数学の成績...etc)。このようなデータセットは散布図によって図に表すのが適当だろう。例えばこのように(以下Rを使って説明します)。

x <- c(1, 3, 4, 5, 7, 2, 8, 9, 10, 6)
y <- c(0, 4, 3, 7, 3, 4, 7, 9, 12, 8)
plot(x, y, pch=16)

ここでは、xとyと名づけられたそれぞれのデータセット(それぞれを変数と呼ぶ)の間には、片方が大きければもう片方が大きいといったような関係、正の相関があるように見る。

ここで、yをxの関数として予想できたら便利と思うかもしれない。

例えば、次のようなモデル(特にこれは一次式による線形モデルとか呼ぶ)を考えよう。

ここでaとbというのはまだ決まっていない定数で、xiというのが個々のデータとなっている(添え字のiには1, 2, 3...と言ったデータ番号が入る)。順番はともかく、いわゆる普通の一次式だ。yの頭に傘が付いているが、これはxiに各々対応するyの値を正確には表現できない(n個のyを完全にxの式で表示したければn-1次式が必要となる)ためで、yの推定値であることを主張している(ちなみに「ワイハット」と読もう)。

さて、まずは見た感じで決めてみよう。こんな感じでどうだろうか。

a=0, b=1としてみた。ためしにこの式の線を引いてみよう。

なかなかのものではありませんか。

まあ、場合によってはこれで十分かもしれない。しかしながら、例えばa=0.01, b=0.99と言ったような値と、今の値とどっちのが「より適切」であるかはどのように判断するべきだろうか。

ここで判断の指標には今引いた「線」から「点」までの距離、特に垂直方向の距離を使う。線をグラフに書き入れよう。

この垂直方向に引かれた線、その合計が最も少なくなるような線が、最も適切な線であるとするのはそれなりに合理的だろう。というわけで、この垂直方向の線が最も短くなるような線を探索することを目標としよう。

まず、「垂直方向の線」を数式で表現する。これはy軸方向におけるモデル式とデータポイントの間の距離であり、誤差、あるいは残差と呼ばれる。残差をdという文字で表すことにする。あるyの値yiに対するdは添え字をつけてdiとする。つまり、

そうして問題はこのdiの絶対値の総和を最小にしようということだった。最小化するためにいじれるのはaとbという値だけで、それらはy[ハット]を表現するモデル式の中に含まれている。ところで「絶対値の」というのは関心があるのが差ではなく「距離」だからである。しかし、絶対値を含む数式というのは非常に計算が難しくなるため、ここではdを二乗することにより符号を正に揃えるという戦法を採ろう。

よって、目的は次の式を最小化するaとbの値を求めることに置き換えられた。で、置き換えられたついでにちょいちょいと変形してaとbを表に引きずり出してやろう。

展開して多少見づらくなったかもしれないが、ここで自由に変更できる文字はaとbだけであることを思い出すと、この式はaとbについてそれぞれ2次式になっていることに気づく。そして2次式として眺めれば、それぞれ下に凸であるということにも気づくだろう。

下に凸な2次式の最小化といえば…「微分して0とおく」というやつ。つまり、変化率が最小になるような点が極値であり、最小値となる。

「微分して」と言ってみたがここでは変数がaとbの2つある。変数が複数ある場合、偏微分という方法を使う。偏微分と聞くと構える人が居るかもしれないが、注目する変数以外を定数とみなして普通に微分するだけなので心配することはない。要するに、bを定数とみなして微分(aについて偏微分)、aを定数とみなして微分(bについて偏微分)の2回の微分をする。すると、次のように2つの式が得られる。

(ここの∂というのは偏微分であることを示す微分記号で、デルだとかラウンド・ディーだとか呼ばれる)

aを変化させたときの変化率、bを変化させたときの変化率というわけなので、こいつらを両方0と置く。

そしてこの2式を連立させ、aとbについて解けば問題は終了。要するに「解の公式」を求めるのが目的なのだ。以下は面倒なだけで特にこれと言って難しい部分もないので読み飛ばし推奨。

まず、式(1)より次式を得る。

式(2)へ代入。

bについて解く。

見やすい感じにする。

次に式(2)より次式を得る。

式(1)へ代入。

aについて解く。

見やすい感じに。

さて、これでaとbに関する「解の公式」が得られた。並べてみよう。


特にこの「形」が意味するものも無いのでさーっと眺めてもらえば問題ないのだが、とりあえず生データのみで答えが出せるというのがポイント。そして実際に計算すべき量は5つしかないのもポイント(よく観察!)。

まあとにかくこの辺は統計ソフトがちょいちょいとやってくれるので、流れだけ把握しておけばよいかなとか思う。モデルは一次式とも限らないわけだし(といっても根拠がなければ一次式にすべきといわれる)。

それで、このように残差の二乗を最小化することにより最適な(とは限らないのだが)モデルを求める手法を最小二乗法と呼ぶ。

最後に最小二乗法により得られたデータ直線を最初の方のグラフに重ねてプロットしてみよう。

まあ最初の「見た感じモデル式」もボチボチではありませんか。手元に電卓とかないときはこういうのも重要ですよ。

ちなみに、Rだと最小二乗法による回帰はこんな感じでやります。

lm(y~x)

これだけで切片と傾きが!で、もう少し詳しい情報が欲しければ

summary(lm(y~x))

と。傾きや切片の有意性、決定係数などといった情報を与えてくれます。

ちなみに最小二乗法は1805年、ルジャンドル先生が発表したわけだけれど、その10年前にはガウス先生が発見していたといういわくつき(?)の方法。ガウス先生はどうもメモに書き付けて満足して発表しないといったことがしばしばあったそうで、先生がちゃんと論文書いてたら数学の発展は何十年も早かったろうとかなんとか。

まあ、要するに論文はちゃんと書きましょうという話でした。

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

部分と全体を比較する

全体の比率を期待値にし、部分部分の比率がどれだけ期待値から食い違っているのかを観察するのが独立性の検定。

部分内のバラつき(グループ内の変動)を、部分を一つのデータと見たときの全体のバラつき(グループ間の変動)と比較し、「部分の塊のバラつき」が「部分内のバラつき」よりどれだけ大きいかを観察するのが分散分析。

似ているようだけど考え方の順序は逆だ。片方は全体から部分へ、もう片方は部分から全体へ。

「全体」から外れているからそれぞれは異なる。

「部分」よりもばらついているからそれぞれはばらついている。

「全体」と「部分」を同一のものさしで測っていいのか、ということは少し気になる。

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

三角関数などの根を小さい方から順に求める[Maxima]

前回RでやったのをMaximaで。

root:[0,0,0,0,0,0,0,0,0,0];
f(x):=cos(x)*(x*3);
myfindroot(f):=block([j:0,a:0.1,k:0],
 loop,k:k+1,
 for i:1 while f(j)*f(j+a)>0 do(a:a+0.1),
 root[k]:find_root(f(x),x,j,j+a),
 j:root[k]+0.1,
 if k=10 then return(root),
 go(loop));

前回は関数のエラーを利用していてスマートじゃなかったのでちゃんと符号を見て関数を使うようにした。 まーでも範囲を0.1ずつ増やしているので、0.1の範囲にいくつも根があると見逃してしまうだろう。もっと上手い方法はないものかしら。

しかしMaximaは使い方がよくわからないなー。一冊くらい本を買って勉強しないとだめかな。 Rみたいに任意の長さのベクトルをパッと作る方法は無いんだろうか。あとエディタからパッと命令を送る方法は無いんだろうか。 エディタをいじらないと駄目な気がしているけどどうなんだろう。やっぱRが一番使いやすいなー。

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

ブートストラップとジャックナイフ

うーん。分かってきたような分からないような。

以前ブートストラップに関してちょっとだけ書いた。 あのときから別に大して詳しくなったわけではないが、講義の関係でクラスタ解析なんてものを少し勉強している。 DNAの変異だとかmRNAの発現の程度だとかいった情報から種や遺伝子のクラスタを作るときは、 その配列をブートストラップ法によって複製し、系統樹が支持される確率なんてのを計算するのが今では当たり前になっているそうだ。 この辺詳しく解説したらプレゼンの時間が稼げるかなとか企んでいる次第で。

まあそれはそれとして、ちょっとジャックナイフ法とブートストラップ法について勉強したのでメモ。 なんか前と同じような内容になりそうだ。まあいいや。

1.推定値

母集団の分布というのは、平均だとか分散だとかいった特性値によって特徴付けられる。もちろん普通はそれを直接に知ることはできない。 母集団の全てのデータを得ることは大抵の場合不可能だろう。そこで、 多くの場合は母集団からサンプリングして得られた標本より母集団の特性値を推定する。

母集団の平均値が知りたければ標本の平均値を計算するし、母集団の分散が知りたければ標本の分散を計算するだろう(もちろん、 分散の場合は標本分散ではなく普遍分散が推定量としては適切である)。

ここで問題となるのは、これらの推定値をどの程度信頼してもよいのか、つまり、これらの推定値はどれだけばらつきうるのか、 ということだろう。

平均値については前回ネタにしたので、今回は分散を扱うことにしよう。

2.カイ二乗分布を用いた分散の区間推定

詳細は面倒なので省くが、分散の区間推定はカイ二乗分布を用いて行われる。カイ二乗分布を利用した場合の95%信頼区間は、 データの平方和をS、自由度をfとして

20080615eq1

で表される。分母はそれぞれ自由度fのカイ二乗分布における2.5%点、97.5%点の値である。

しかし、この場合は母集団の正規分布を仮定している。母集団が正規分布しているとはいえない場合、 導かれる区間は適切なものとならない場合がある。

3.ブートストラップ法を用いた分散の区間推定

ブートストラップ法を用いる場合は母集団に正規分布のような分布を仮定しない。分散の区間は、 サンプリングからの複製によってその場で作り出される。詳細は前回解説したので省略。

4.ジャックナイフ法による分散の誤差推定

ジャックナイフ法とは1940年代にM. Quenouilleにより考案され、1950年代にJhon Tukeyにより命名された手法である。ブートストラップ法と異なるのは、 リサンプリングにおいていくつかのデータを抜いた状態のサンプルを生成するということ。狭義には1つのデータが除かれ、 重複サンプリングを行わない。そのため、ブートストラップ法よりも計算量は少ない。

4.1 ジャックナイフ法によるバイアス推定

推定値にはバイアスがかかる場合がある。たとえば、母分散の推定値として標本分散を採用することは妥当なことのように思えるが、 よく知られているように母分散の不偏推定量は不偏分散である。標本分散を母分散の推定値としてみた場合、 サンプル数に依存したバイアスがかかっているということになる。

n個のサンプル全てを用いて計算された推定量をθとする。n個のサンプルから、 i番目のデータを除いて作られた複製より計算された推定量をθiとする。その平均をθmeanとする。このとき、 ジャックナイフ法によりバイアスを修正された推定量は次式で計算される。

20080615eq2

たとえばこの式へ標本分散を代入して計算してみると、不偏分散となる。

4.2 ジャックナイフ法による誤差の推定

i番目のデータを除いて計算された推定値をθiとして、θの分散は

20080615eq3

によって計算される。

分散が分かるんだから信頼区間も分かるんだろうけどちょっと求め方がよくわからない。 θは中心極限定理から正規分布すると考えて良いんだろうか。

5. Rによる計算

ともかく、Rを用いて確認してみよう。

まず母集団の用意。

x <- 1:10

1〜10までの値が等しい確率で出現する離散分布と考えていい。分散は8.25、平均は5.5の一様分布となる。

ここから20個のデータからなる標本を抽出する。

#サンプルの生成
x.samp <- sample(x, 20, replace=T)

そして、カイ二乗分布を利用した場合とブートストラップ法による場合の区間推定。

#カイ二乗分布を利用した区間推定
19*var(x.samp)/qchisq(c(0.975,0.025),19)
#bootstrap法による区間推定
x.var.samp <- numeric(2000)
for(i in 1:2000){
 x.b <- sample(x.samp,20,replace=T)
 x.var.samp[i] <- var(x.b)
 }
quantile(x.var.samp,c(0.025,0.975))

次に、ジャックナイフ法によるバイアス修正と、「分散の分散」の推定。

#ジャックナイフ・バイアス修正済み推定量
x.j.var <- numeric(20)
for(i in 1:20){
 x.j.samp <- x.samp[-i]
 x.j.var[i] <- sum((x.j.samp-mean(x.j.samp))^2)/19
 }
sum((x.samp-mean(x.samp))^2)-19*mean(x.j.var)
#ジャックナイフ・分散推定
(19/20)*sum((x.j.var-mean(x.j.var))^2)

手元で一回実行してみたところ、

  • サンプルの不偏分散: 7.671053
  • ジャックナイフ・バイアス修正済み標本分散: 7.671053
  • カイ二乗分布による信頼区間(95%): 4.43652 - 16.36442
  • ブートストラップ法による信頼区間(95%): 4.365592 - 10.450526
  • ジャックナイフ法による分散の分散の推定: 2.550627
  • ブートストラップ法による分散の分散: 2.371628

となった。まず、サンプルの不偏分散とジャックナイフ・バイアス修正済み標本分散が同じになるのは当然。計算の中身が同じなので。

カイ二乗分布による区間はブートストラップによるものより広く出ている。カイ二乗分布による区間はほぼ100% の確率で母分散である8.25を含んだ。母集団が正規分布でない影響だろうか。一方、 ブートストラップ法による区間は9割前後の確率で母分散を含んだ。2000回程度の試行をしてみたが、 母分散が信頼区間に含まれる確率がどうも95%よりも低い気がする。なんでだろ。

ジャックナイフ法による分散の分散の推定とブートストラップ法による分散の分散は大体同じになった。うーん。 両者がそれほど離れることはないし、これだとジャックナイフ法で十分じゃないかということになるなぁ。分散だから? ブートストラップは適用範囲が広いとかいう話も見たなー。

まーとにかく今日は眠いので寝よう。参考文献はまたあとで追加。追加しました。

参考

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

円分方程式

20080609eq1

以前、Rを使って正n角形を描く関数を定義したとき、 上記の円分方程式というものを利用した。この方程式の根は、複素平面上の単位円の円周をn等分する、というものであった。

この方程式の根は、つまり1のn乗根ということで、代数学の基本定理よりそれはn個存在する。

たとえば、次の3次方程式の場合ならば、

20080609eq2

この式の根は次の3つである。

20080609eq3

そしてこれを複素平面上にプロットし、それぞれ直線で繋げば正三角形が描かれる。

tri

そして、この根はnがいくつだろうと次の式によってすべて求められる。

20080609eq10

以下、ちっとだけ詳しく解説。

まず、この複素平面というやつだが、「a+bi」という形で表される複素数の、aを横軸の「実軸」に、bを縦軸の「虚軸」 に対応させることで平面上に複素数を表現している。平面上のひとつの点が、ひとつの複素数に対応する。

そして、この平面上の単位円ということだから、単位円周上の複素数のaとbはすべて

20080609eq4

を満足するということになる。左辺の式で表現されるものは複素数の絶対値と呼ばれる。つまり、 絶対値が1の複素数の集合が複素平面の単位円を形成する。また、複素数の絶対値とは極座標表現における「長さ」 であることもわかる。

ここで、原点から複素数へ引かれた直線と実軸とがなす角、複素数の絶対値の2つの情報があれば、三角関数を用いてaとbを表現できる。 絶対値をr、なす角をθとし、aとbはそれぞれ次のように表現できる。

20080609eq5

ここで角θは偏角と呼ばれる。ある複素数zの偏角である、ということを意識したい場合は、「θ= arg z」などと書く。複素数zは、aとbの代わりに偏角と絶対値を用いて次のように表現できる。

20080609eq6

このような複素数の表現方法を極形式と呼ぶ。極形式を用いて、複素数同士の積を考察する。

20080609eq7

最後の変形には三角関数の加法定理を用いている。この変形から明らかなように、 複素数同士の積の絶対値はもとの複素数の絶対値の積であり、 偏角はもとの複素数の偏角の和となる。

単位円の場合、絶対値は1であるため、偏角θを持つ単位円上の複素数を乗ずるということは、 複素平面上におけるθ分の回転を意味する

話を1のn乗根へ戻そう。1のn乗根ということは、n回乗ずると1になるような数、ということである。 先の複素数の積に関する考察を利用して言い換えれば、 n回回転させると角度が0(=2π)となるような角θを偏角に持つ絶対値1の複素数、ということである。 単位円上の複素数をn回乗ずると、偏角はn倍される

さて、n倍すると2πになるような角はθ=0を含め、n個存在する。それらはすべて2π/nの整数倍である(よく確認してほしい)。 このことに注意すれば、

20080609eq1

の根は

20080609eq8

で表されるn個の複素数であるということがわかる。

ここに、オイラーの公式を適用して…ってブログ内に記事がないじゃないか。とっくに書いてるものだと。 仕方ないので今日のところは天下り式に次の式を。(オイラーの公式について簡潔にまとまっていて分かり易い解説:オイラーの公式[物理のかぎしっぽ])

20080609eq9

これがオイラーの公式と呼ばれるもので、θにπを代入したりすると大変なことになったりするが今は触れない。 オイラーの公式と上の式を見比べると、1のn乗根は次の式で求められることがわかる。

20080609eq10

まあ書き方がシンプルになっただけで中身は変わっていない。オイラーの公式を使うと、 極形式の複素数を指数の形でひとつにまとめられる、ということである。

それで、以前Rを使って正n角形を書く関数を定義したときはこの最後の式だけを使ったわけだが、 当時は今ひとつこの式の意味がわからないまま使っていた気がする。

さっき寝ようと思って布団の中で本を読んでいたら1のべき根の話が書いてあって、 1のn乗根を求める上の式の意味がなんとなく分かってきたので、いろいろなことを棚上げして夜更かしして記事にまとめた次第。

参考文献

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

比率の検定と角変換

中心極限定理によれば

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

とされる。だから、サンプリングされたデータの平均値について、「もとの集団の平均値は○○% の確率でサンプルの平均値±△△の間にある」とか言える。が、 平均値の分布が正規分布に従うために必要なサンプル数というのはかなり多いので、 多くの場合はサンプルが少なくても使えるt分布というものを使う。 t分布はサンプル数が少ないときに若干すそが広いという特徴があるが、サンプル数が多くなればほとんど正規分布と同じになる。

いずれにせよ、平均値の分布する幅を推定するには元となる集団の分散、母分散(t分布の場合は推定された分散である不偏分散) が必要となる。そして、平均値の差についての検定では、比較するサンプルの母分散が等しいという場合において、 2つのサンプルの差がこれだけになる確率はどれだけなのか?ということを問うのである。

要は、平均値の分布の仕方に基づいた「平均値の差の検定」を使うためには、 母集団の分散が等しい必要がある、ということである。

「母集団の分散」なんてものは実際には推定するしかないわけで、本当に等しいかどうかはわからない。それに、 平均値に依存してばらつきの仕方が変わるかもしれない。たとえば、小さいものの大きさよりも大きいものの大きさの方がばらつきが大きい、 というようなことは大いにありそうである。ありそうだが、差の検定をする場合はそんなことは無視する。無視しても大して問題はない、 と前向きに考える。

しかしながら、「平均値」というものが二項分布に従う「比率」である場合(たとえば、コインを10枚投げたときの表が出た枚数で、 このような試行はベルヌーイ試行と呼ぶ)は、そんな風に前向きに考えていてはいけない。なぜなら、「比率」の分散は「ありそう」ではなくて 「明らかに」比率に依存している。比率の分散の計算に試行が成功する確率が入ってくるのが原因だが、直感的には 「ほとんど成功する、ほとんど失敗するような試行ほどばらつきが少なく、成功と失敗が半々となるような試行でもっともばらつきやすい」 と考えると理解しやすい。上限と下限があるものだから、ばらつきが押さえ込まれるのだ。 天井効果とかフロアー効果とか言うらしい。

そのため、二項分布に従うようなデータについてt分布や正規分布を利用して平均値とそのばらつきを推定、差の検定をしてみても、 比率によってばらつきの仕方が異なるので検定の妥当性が問題になることがある。

実際に、成功と失敗の割合を0〜100%まで変化させて行ったベルヌーイ試行(試行回数20)において、 その平均値のばらつきを次に示す。

 Rplot001

50%付近でばらつきがもっとも大きくなる。平均のばらつき方そのものが実際にこのような変化をしているので、 t分布や正規分布を用いて平均値の分布を正しく予想できたとしても、それを差の検定に使うのが難しいのだ。

では、平均値から計算できるもので、比率によらず分散が同じになるようなものはないだろうか。そういうものがあるならば、 それの値を比較することで平均値の差の比較の代替とすることができるはずだ。

で、そういうものは実際にある。それは成功した回数の比率をpとして、次のような変換によって得られる値だ。

temp1

arcsinというのはsinの逆関数で、0〜1までの値を与えるとそれをsinとして持つ角度を返してくれる。そこで、 この変換は角変換、または逆正弦変換、アークサイン変換などと呼ばれる。

角変換後の値のばらつきを次に示す。

Rplot003

まあおおよそ同じになっていることがわかるだろう。変換前に比べれば大分マシだ。

ところでこの値のばらつきの大きさがわからないと検定に使えないが、 この値はサンプル数をnとして分散が1/4nとなることが知られている。よってこれらの情報を用いて統計処理をすればいいということになる。 なお、比率が0の場合と1の場合は当然分散が0となるのだが、0の場合は1/4nを、 1の場合は1-1/4nを分散の値として用いるという決まりになっている。

これにてめでたしめでたし。

なのだが、もう一度この図を見てほしい。

Rplot001

「別に中央付近だったら値はそれほど変わらないじゃないか。」

そのとーり。30〜70%くらいまでの範囲だったら別に角変換をしてもしなくても結果が変わることはまずない。

そもそも分散が等しくないことが問題なのだから、等分散という仮定をおかない(ノンパラメトリックな) 統計処理を行えばそれで事足りる。

しかしながら、等分散を前提条件として必要とする多重比較や分散分析を行いたいという場合もあるかもしれない。 そういうときには角変換をしたほうが良いような気がするが、別に無視して普通に統計処理した結果を参考にしても問題はないと思う。ただ、 誰かに発表するような場合には「比率の統計=角変換」 といった断片的な情報だけを持っている人間を黙らせるために効果を発揮するかもしれないので角変換をしておくといいかもしれない。

重要なのはその「差」にどんな意味があるのか、ということだろう。「成功率が1%違った。統計的に有意な差があった。」 と言ってみても、「1%の差」に統計的以外の意味がなければ、「そんな小さな差がなんだ」という話になってしまう。そして、「1%の差」 が統計的に有意となるか否かは、サンプル数次第なのである。 完全にidenticalな処理などというものは存在し得ないという立場に立てば、ほとんど同一に見える2集団から得られたデータでも、 サンプル数次第でほぼ確実に「有意な差」が検出される。大量のサンプルにより導かれた「統計的な差」というのは大した意味をもたない。

角変換は重要で有用な手法なのだが、統計的な結果に対する説得力をわずかに増す以上の効果はない、ということは忘れないようにしたい。

グラフも見ずに「比率の検定?角変換しないと!」みたいな条件反射はよくない。まずは標準誤差を付けたグラフの観察からすべき。 値次第では、手元の道具(Excelとか)を使って普通の平均値の差の検定をしたっていい。

最後に今回使ったRのスクリプトを。

Xvar <- numeric(101)
Xarcvar <- numeric(101)
Xmean <- numeric(101)
X <- numeric(200)
Xtheta <- numeric(200)
for(j in 0:101){
 x <- c(rep(1,j),rep(0,101-j))
 for(i in 1:200){
a <- sample(x,20,replace=T)
X[i] <- mean(a)
Xtheta[i] <- acos(sqrt(mean(a)))
}
 Xmean[j+1] <- mean(X)
 Xvar[j+1] <- var(X)
 Xarcvar[j+1] <- var(Xtheta)
 }
plot(Xvar)
plot(Xmean)
plot(Xarcvar)

説明は省略。まー動かしてみてください。

参考

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

広告


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

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

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


×

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