« 2009年04月 | メイン | 2009年06月 »

2009年05月 アーカイブ

2009年05月06日

Java3Dで正多面体を作ってみた

Java3Dのシェーディングがよくわからないので、勉強しようとして、そのための実験材料として、ちょっと多面体でも作ってみようと思った。正多面体というのはn=4,6,8,12,20の正n面体しか無いが、正8面体以下は点光源や平行光源の光が当たる面の数が少ないので、正12面体と正20面体にすることにした。

…ら、えらく手こずって、GWの後半をこれだけで潰してしまった。光源とか反射率とか色々実験する計画だったのだが、時間切れ及び肩凝り等のため、一区切りする。

Javaアプレットの起動ページへ
・ソースコード
LightTest1.java
MyHedron.java
MyHedron12.java
MyHedron20.java

・サンプル画像
LightTest11.png

LightTest12.png

ついでに、各面の中心を凹ませたり凸ませたりしてみた。

多面体の作成は、正12面体と正20面体の各面の法線を求めた上で、各面の中心と頂点を結んで三角形ポリゴンを作る方法で行っている。
正12面体の1つの面をXY平面と平行にすると、その面の法線は(0,0,1)であり、その面と隣の面とがなす角度をθとすると、tanθ=2である(説明は省く。調べると所々に出てくる63.4349...というのがこのθである)ので、(0,0,1)を360÷5=72°回転させると別の1つの法線が求まり、それをZ軸周りに72°ずつ回転させていくと、正12面体の上半分の法線が全て求まる。下半分は、それを逆さにしてZ軸周りに36°回転させて上半分とくっつけると求まる。
正12面体の天面の頂点の1つをX軸上に取ると、そのX座標は、sqrt((3/4) G/sqrt(5)) となる(Gは黄金比=(sqrt(5)+1)/2)(いい説明が思いつかないので説明は省く)。正12面体の各面は五角形なので、これを正12面体の法線で72°ずつ回転させていくと、全ての頂点が求まる。
正12面体の面と頂点を入れ替えると正20面体になるので、正20面体の頂点は正12面体の法線ベクトルのX,Y,Z成分をX,Y,Z座標にすることで得られ、正20面体の法線は正12面体の頂点の座標をベクトルの成分とすることで得られる。

光源は、環境光(赤)と平行光(青)を設定し、なんとなく点光源(黄)を回転させてみた。

続きを読む "Java3Dで正多面体を作ってみた" »

2009年05月10日

統計学復習メモ13: 有効推定量と不偏性、一致性

我ながら物好きであるという仮説が棄却されなさそうであるが、大学の時に買う羽目になった統計学の教科書を引っ張り出して、推定量というものを勉強している。そこには、推定量の好ましい性質として不偏性と一致性が挙げられ、次に推定量同士を比べる話が出てきて、その後に有効推定量、最尤推定量が出てくる。この教科書のこの流れのために、有効推定量や最尤推定量は不偏性と一致性を満たすんだろうとずっと誤解していた。何せ「有効」「最*」なのだから。

特に「最尤推定量」の名前の意味は「最ももっともらしい推定量」であり、まるでキャッチフレーズのように「最尤推定」に付け加えられ、そのインパクトによって、これぞ最強の推定量というイメージが刷り込まれてしまう。大学の研究室にも、最尤推定量の意味を理解せず、「推定」という言葉が出てくる度に、馬鹿の一つ覚えのように「最尤推定を使えばいいじゃない」と言っていた教官がいた。実データが取りようが無い(統計量が無い)ものの話をしていても「最尤推定」という単語を持ち出されて、私を含め、学生たちは皆、そのセリフに降参するしか無かった。

余談はさておき、よく読むと、どうやら有効性(efficiency)と不偏性(unbiasedness)、一致性(consistency)とは直接は関係ないようである。例えば標準分散は一致推定量であるが不偏ではない(不偏分散は不偏かつ一致推定量)し、逆に不偏だが一致性が無い場合もある。有効推定量や最尤推定量にも不偏でない(biasedな)ものがあるし、最尤推定量にも一致性が無いものがあったりするらしい。


一致性というのは、数式で書くと、実際のパラメーターをθ、推定量をθ^、Pを確率、nを標本数として
¥exists¥varepsilon ¥lim_{n¥rightarrow¥infty}P(|¥hat{¥theta}-¥theta|>¥varepsilon)=0(または不等号を逆にして=1)
とちょっとアレであるが、要するに標本の数が多ければ推定量が実際の値に近づく性質のことである。意図的に外すかよっぽどチョンボしない限りはそれを満たすのが当たり前な気がするが、推定量が多項式でない場合は気にした方が良いのかも知れない。

不偏推定量であるが一致推定量でない極端な例としては、標本をX1...Xnとして、平均θの推定量θ^=X1やθ^=Xnがある。標本がいくら増えてもその中の1つしか使わなければ精度が上がらないので、当たり前である。標本を全て使っても一致推定量でない例として、こういうのを考えてみた。
¥hat{¥theta}=¥sum_{k=1}^{n}¥frac{X_k}{2^k}+¥frac{X_n}{2^n}
これは、実際の分散をσ2とすると
¥lim_{n¥rightarrow¥infty}Var(¥hat{¥theta})=¥lim_{n¥rightarrow¥infty}¥left(¥sum_{k=1}^{n}¥frac{¥sigma^2}{2^{2k}}+¥frac{¥sigma^2}{2^{2n}}¥right)=¥frac{¥sigma^2}{3}
であり、nを大きくしても分散が0にならないので、θ^はθに収束しない(はず)。(ちなみに、普通にθ^を標本平均(X1+X2+...+Xn)/nとすると、Var(θ^)=σ2/nなので、¥lim_{n¥rightarrow¥infty}Var(¥hat{¥theta})=0となり、θ^はθに限りなく近づく)


有効推定量というのは、クラメル=ラオの不等式
Var(¥hat{¥theta}) ¥geq ¥frac{1}{n E¥left¥{¥left[¥frac{¥partial}{¥partial¥theta}¥log P(X)¥right]^2¥right¥}} ¥equiv ¥sigma_0^2
の等号が成立する(Var(θ^)が「クラメル=ラオの下限」σ02になる)ようなθ^のことである。確率の対数をθで偏微分して2乗して平均して逆数にして訳がわからないが、対数が出てくるのは、AとBとCが同時に起こる確率がP(A)×P(B)×P(C)なのでその対数がlogP(A)+logP(B)+logP(C)となるように、標本それぞれが同時に起こる確率が線形和になるためのもので、偏微分は飛ばして、2乗の平均は誤差の2乗和の意味で出てくるのだと思う。

上記の不等式を解く代わりに、
¥hat{¥theta}=¥theta+K¥sum_{i=1}^{n}¥frac{¥partial}{¥partial¥theta}¥log P(X_i)
の右辺がθを含まなくできるK(但しKはXiは含まない)があれば、その時のθ^が有効推定量、という定理が使える。それを使って、X1...Xnが平均μ、分散θの正規分布に従う標本の場合のθの有効推定量を求めてみる。確率密度関数は
P(x)=¥frac{1}{¥sqrt{2¥pi¥theta}}e^{-¥frac{(x-¥mu)^2}{2¥theta}}
なので、
¥frac{¥partial}{¥partial¥theta}¥log P(x)=¥frac{¥partial}{¥partial¥theta}¥left(-¥frac{1}{2}(¥log(2¥pi)+¥log¥theta)-¥frac{(x-¥mu)^2}{2¥theta}¥right)=-¥frac{1}{2¥theta}+¥frac{(x-¥mu)^2}{2¥theta^2}
であり、
¥hat{¥theta}=¥theta+K¥sum_{i=1}^{n}¥left(-¥frac{1}{2¥theta}+¥frac{(x-¥mu)^2}{2¥theta^2}¥right)=¥theta-K¥frac{n}{2¥theta}+K¥frac{1}{2¥theta^2}¥sum_{i=1}^{n}(x-¥mu)^2
なので、K=2θ2/nとするとθが消えることがわかる。従って、
¥hat{¥theta}=¥theta-¥frac{2¥theta^2}{n}¥frac{n}{2¥theta}+¥frac{2¥theta^2}{n}¥frac{1}{2¥theta^2}¥sum_{i=1}^{n}(x-¥mu)^2=¥frac{1}{n}¥sum_{i=1}^{n}(x-¥mu)^2
であり、標本分散がθの有効推定量である。
標本分散は不偏ではないので、正規分布の分散の有効推定量は不偏推定量ではないことがわかる。

続きを読む "統計学復習メモ13: 有効推定量と不偏性、一致性" »

2009年05月14日

統計学復習メモ14: ポアソン分布のパラメーターの最尤推定量

筆者はある会社のソフトウェア関係の部署に在籍しているが、一般にはソフトウェア開発とは呼ばれない補助的な仕事をしている。調達担当や購買担当という呼び方が最尤推定量であろう。そんな筆者にもお呼びがかかるのが、不再現バグの再現試験である。
バグ発生の報告があったが、開発担当者が解析するために再現させようとしても、再現しない。発生した証拠はある。となると、総出で再現試験である。

そのバグを発生させるべく、100回も200回も同じ操作をしていると、いつまでやらないといけないのだろうと不安になってくる。1万回やっても出ない可能性もあるのである。仮に発生確率が1/1000だったとしても、1000回やったら必ず1回は出る訳ではない。確率が低い事象の発生回数は「小数の法則」ポアソン分布に従う。発生回数(発生確率×試行回数)の平均がλの場合、発生回数がkになる確率は、
P_k=e^{-¥lambda}¥frac{¥lambda^k}{k!}
である。
λ=1のポアソン分布の確率密度
従って、1000回やって平均1回しか起こらないものが1000回やって1回も起こらない確率は1/e≒36.8%もある。1回起こる確率が同じく1/e、従って2回以上起こる確率が(1-2/e)≒26.4%であるが、1回も起こらない確率が36.8%もあるのだ。既に何百回も起こらないでなぜ+思考になれようか。

自分が毎回その1/eに入るのを不思議に思いながらも、誰かが再現させて解放される。最初に再現させる人は1人しかいないのだから、自分がその1人になる確率は低いので、左脳で考えると当たり前である。複数人でのべ2000回やって1回起こったとすると、結果として発生率は1/2000である。
さて、本題である。この時、本当の発生率も1/2000と推定できるのだろうか?実は本当の発生率が1/5000の時が2000回に1回発生する確率が一番高かったりしないのだろうか?
これはまあ、P(k=1)=λe

こんな形をしてるので、P(k=1)をλで偏微分して=0とすると、(1-λ)e=0でλ=1ということになり、1回発生する確率が最も大きいのはλ=np=1すなわちp=1/2000と力技で求められる。

では、もし、ある人は200回で発生させ、ある人は500回で発生させ、ある人は1000回で発生させたら、その発生率は全て足し合わせて3/1700と推定するものなのだろうか?
問題を整理すると、ポアソン分布のパラメーターλの推定量λ^は何なのだろうか?


今回は、最尤推定法を使うことにする。最尤推定法とは、あるパラメーターθを持つ確率分布を仮定して、結果として標本X1...Xnが得られた時、n個の標本を取ると前記の標本X1...Xnが得られる確率が最も高いθを求める方法である。X1...Xnが得られる確率をP(X1)...P(Xn)とし、X1...Xnが独立(1つ前の標本によって発生率が変わったりしない)とすると、それらが同時に起こる確率はL(θ)=P(X1)×P(X2)×...×P(Xn)である。これが最大になるθ=θ^を求めるのである。さらにL(θ)が単峰(山頂=極大点が1つしかない)なら、そのθは
¥frac{¥partial}{¥partial¥theta}L(¥theta) = ¥frac{¥partial}{¥partial¥theta}¥prod_{i=1}^{n}P(X_i) = 0
の解である。積の形だと何なので、対数を取って
¥frac{¥partial}{¥partial¥theta}¥log L(¥theta) = ¥frac{¥partial}{¥partial¥theta}¥sum_{i=1}^{n}¥log P(X_i)=0
の解でもあるのであるのである。

ポアソン分布のλについては、
¥frac{¥partial}{¥partial¥lambda}¥log P_k = ¥frac{¥partial}{¥partial¥lambda}¥log e^{-¥lambda¥frac{¥lambda^k}{k!}} = ¥frac{¥partial}{¥partial¥lambda}(-¥lambda+k¥log¥lambda-¥log(k!)) = -1 + ¥frac{k}{¥lambda}
なので、L(λ)が最大になるλを求めると、
¥frac{¥partial}{¥partial¥lambda}¥log L(¥lambda)=¥sum_{i=1}^{n}(-1+¥frac{X_i}{¥lambda})=0
-n+¥frac{1}{¥lambda}¥sum_{i=1}^{n}X_i=0
¥lambda=¥frac{1}{n}¥sum_{i=1}^{n}X_i=¥bar{X}
となるので、標本平均が最尤推定量である。

元の問題に戻って、X1=1/200、X2=1/500、X3=1/1000と置くと、平均λ=(X1+X2+X3)/3=8/3000=1/375となる。

続きを読む "統計学復習メモ14: ポアソン分布のパラメーターの最尤推定量" »

2009年05月16日

統計学復習メモ15: 幾何分布のパラメーターの最尤推定量

前のメモで間抜けなことをしたついでに、同じ問題についてもう少し考察する。

再現率が低いバグを、(1) Aさんは200回試行して1回発生させ、Bさんは500回試行して1回発生させ、Cさんは1000回試行して1回発生させた場合と、(2) Aさんは200回目で1回発生させ、Bさんは500回目で1回発生させ、Cさんは1000回目で発生させた場合とでは当然数字の意味が異なる。(1)では発生させた後も試行してる可能性があるからである。

発生させたら試行を終了すると考えた場合、発生率をpとすると、k回目で発生する確率は、(k-1)回失敗して1回成功する確率なので、p(1-p)k-1である。簡単にするため、k+1回目で発生する確率をP(k)=p(1-p)kとする。このkの分布には「幾何分布」という名前がついているらしい。幾何級数(等比級数)だからということらしいが、この用語は一般に通じるものなのかどうか疑問だが、今回はこの用語を使う。発生するまでの試行回数の標本X1〜Xnが得られた時、それが幾何分布に従うとすると、発生率pの推定量は何になるだろうか?

前回に続き、最尤推定法を使う。幾何分布の確率をθ、尤度関数(結合確率)をL(θ)とすると、
¥frac{¥partial}{¥partial¥theta}¥log L(¥theta) = ¥frac{¥partial}{¥partial¥theta}¥sum_{i=1}^{n}¥log(¥theta(1-¥theta)^{X_i})
なので、これが0になるθを求める。
¥frac{¥partial}{¥partial¥theta}¥sum_{i=1}^{n}¥log(¥theta(1-¥theta)^{X_i})=¥sum_{i=1}^{n}(¥frac{1}{¥theta}-¥frac{X_i}{1-¥theta})=¥frac{n}{¥theta}-¥frac{1}{1-¥theta}¥sum_{i=1}^{n}X_i=0

¥theta=¥frac{n}{n+¥sum_{i=1}^{n}X_i}=¥frac{1}{1+¥bar{X}}
となる。

従って、上記(2)の場合だと、X1=199, X2=499, X3=999なので、p=3/1700と推定される。
(1)の場合は、前のメモの通り、(1/200+1/500+1/1000)/3=1/375が妥当だと思う。

About 2009年05月

2009年05月にブログ「Weblog on mebius.tokaichiba.jp」に投稿されたすべてのエントリーです。過去のものから新しいものへ順番に並んでいます。

前のアーカイブは2009年04月です。

次のアーカイブは2009年06月です。

他にも多くのエントリーがあります。メインページアーカイブページも見てください。

Powered by
Movable Type 3.35