Jan 28, 2010

統計学復習メモ18: 少数のベルヌーイ試行による発生確率の区間推定

問:ある機械が20個の部品を作ったら、不良品が2個出てきた。この機械が1個作る毎に不良品を発生させる確率は一定だとすると、不良品発生率は信頼度95%の信頼区間として何%〜何%の間と推測できるか。

観測誤差を考えなければ、すなわち点推定量としては、推定確率はもちろん2/20=0.1である。しかし、標本はたった20個である。これをもって、だから2000個作ったら200個の不良品が出ることが予想される、と言うのはあまりにも乱暴だし、直感的にナンセンスであろう。サンプルが少なすぎて信憑性が低く、説得力が無いのである。
信頼性を高めるにはサンプル数を増やせば良いが、例えば以前のエントリーに書いたような計算方法では、推定確率が0.1前後の場合に誤差を±0.05以内にするには、(1.96/0.05)2*0.1*0.9=約140個のサンプルが必要になる。サンプル数を先に決められない場合や、サンプルはそれしかないという場合は、推定発生率はこの範囲に入る、という信頼区間を求めるしか無い。

観測されたベルヌーイ試行の成功回数(発生回数)からの元の成功率(発生率)の区間推定については、以前のエントリーに書いたように信頼区間の公式がある。
I \in \hat{p}\pm Z\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} ...(1)
(p^は点推定量、Zは標準正規分布N(0,1)に従う値、nはベルヌーイ試行の標本数)

従って、p^=0.1、Z=1.96(右側2.5%点)、n=20とすると、95%信頼区間は0.1±1.96√0.0045 = 0.1±0.131、すなわち-3.1%〜23.1である。(完)

イエース、公式一発、ビバ統計学!
…ちょっと待て、確率なのにマイナスの値?何を間違えたのだろう?

発生率が一定で、1回の試行の結果は発生するかしないかの2通りしか無い試行(ベルヌーイ試行)を複数回行った時の発生回数は二項分布に従う。そして試行回数が十分大きければ二項分布は正規分布で近似できる。試行回数をn、発生率をpとすると、発生回数は平均np、分散np(1-p)の正規分布に(近似的に)従う。観測される発生率(推定確率)は発生回数を試行回数で割ったものだから、平均p、分散p(1-p)/nの正規分布に(近似的に)従う。母分散が既知の場合、正規分布に従う標本の推定値は毋分散÷標本数を分散とする正規分布に従うとするのが区間推定の定跡である(母分散が未知の時は標本分散を使うのでt分布に従うとする)から、信頼区間Iは
I \in \hat{\theta}\pm Z\sqrt{\frac{\sigma^2}{n}}
(θ^は点推定量、Zは標準正規分布のZ値、このnは正規分布に従う標本の数なので今回は1)
で与えられる。冒頭の問題については、θ^=p、σ2=p(1-p)/nなので、信頼区間はやはり(1)の公式の通りである。

という訳で、二項分布を正規分布に近似しているのが問題のようである。
次のグラフは、n=20として、p=0.5とp=0.1の二項分布とN(np, p(1-p)/n)の正規分布を比較したものである。

青線が二項分布、赤線が正規分布である。p=0.5だとかなり近いが、p=0.1だと差がある。比率が0.5から遠いほど、二項分布の正規分布への近似は悪くなるのである。
大体、正規分布に近似するということは0より小さい確率や1より大きい確率を認めるということであり、n=20でp=0.1の場合、正規分布の左側6.8%は負の値に入り、無効になってしまうのである。この6.8%を全てp=0の所に加算すれば(信頼区間の下限が負の値なら0に補正すれば)いいだけの話とも言えなくはないが、pの推定値が0に近づくと急激にp=0の確率が上昇することになるので、ちょっと強引であろう。

上記の比率の信頼区間の公式(1)は、np<5またはn(1-p)<5の場合は不適、すなわちベルヌーイ試行の成功回数(発生回数)も失敗回数(発生しなかった回数)も5以上でないと使えないと言われているらしい。視聴率調査や投票率予想なら最低5人が観て/投票していないと、この公式を使うのは適当でないのである。それに基づくと、冒頭の問題もnp=2なのでこの公式は使えない。

ではこの場合はどのようにして信頼区間を求めれば良いのだろうか。それを解決してくれるのがClopper-Pearson methodによる通称Exact Confidence Interval(Exact CI)である。
区間推定の原点に立ち返ると、信頼区間とは、ある推定値がある信頼度で収まる値の範囲であり、まともに計算するなら事前確率、事後確率を考慮して計算すべきものであるが、それに代えて、推定の元になった観測値及びそれより稀な値が発生する確率が(1-信頼度)以下であるような推定値の範囲を信頼区間として計算する、というのを二項分布のパラメーター推定について行うのがClopper-Pearson methodである。冒頭の問題なら、95%の確率でpが0.1の前後のどこまでの範囲に含まれるかを求める代わりに、pが0.1よりどこまで小さければ不良品が2個以上出る確率が2.5%になるか、pが0.1よりどこまで大きければ不良品が2個以下である確率が2.5%になるか、を求めるのである。
という訳で、ベルヌーイ試行の確率の推定値の信頼区間の下限Plbと上限Pubを束縛する式は次のようになる。(αは1-信頼度)
\sum_{x=0}^{k}\pmatrix{n \cr x}P_{ub}^x (1-P_{ub})^{n-x}=\frac{\alpha}{2}

\sum_{x=k}^{n}\pmatrix{n \cr x}P_{lb}^x (1-P_{lb})^{n-x}=\frac{\alpha}{2}
これを計算したものがExact CIで、BetaCDF(x,a,b)をベータ分布の累積密度関数として、その逆関数を用いて、
P_{ub}=1-BetaCDF^{-1}\left(\frac{a}{2},n-k,k+1\right)
P_{lb}=1-BetaCDF^{-1}\left(1-\frac{a}{2},n-k+1,k\right)
と書けるらしい。

Maximaにはベータ分布の累積密度関数の逆関数(β値を求める関数)がquantile_beta()という名前で存在するので、これを使って信頼区間の上限を計算できる。

(%i1) load(distrib);
(%i2) Pub(n,k,a) := 1-quantile_beta(a/2, n-k, k+1);
(%i3) Plb(n,k,a) := 1-quantile_beta(1-a/2, n-k+1, k);
(%i4) Pub(20,2,0.05);
(%o4) 0.317
(%i5) Plb(20,2,0.05);
(%o5) .012349

答:1.23%から31.7%の間

参考文献:Understanding Binomial Confidence Intervals

See more ...

Posted at 20:03 in 数学 | WriteBacks (0)
WriteBacks

Jan 09, 2010

PMXで楽譜作成

PMXとは、

こういう楽譜を作るツールである。フリーウェアであること、出力される楽譜が綺麗であることもさることながら、入力データをテキストファイルで書けるのが魅力である。

表示されている楽譜の上にマウス操作で音符を乗せていくようなソフトは多数あるが、入力し易いものに出会ったことが無い。筆者が細かいマウス操作が苦手なこともあるが、五線譜の上から2番目と3番目の線の間に音符を持っていって、次の音符は2番目の線の上に乗せて…という作業は苦痛である。ドラッグしているその音符をあと半音上に動かす必要がある時、今マウスが止まってるので、マウスが反応するぎりぎりの速度を狙ってゆっくりとマウスを前に動かし、マウスが反応したらすぐ手を止めないといけない。手を止めるのが少し遅くて目的地より半音上に行ってしまったら、今度は同様に慎重にマウスをバックさせてカーソルが動いた瞬間に手を止めないといけない。そうこうしていると、ボタンを押す指が疲れて一瞬震え、上か下かに半音ずれた所で音符が確定しまう。それを思い出して文章にしてるだけでイライラしてくる。
昔、MSXの「MUE」(ハル研究所)という、マウスが無かったので当然(パソコンの)キーボードで音符を入力する作曲ソフトを買ったことがあるが、そっちの方がずっと入力し易かった。音符はキーボードで入力できる方が便利なのである。

そこでMusiXTeXである。TeXのマクロであり、当然入力はテキストファイルなので音符はキーボード操作で入力できる。そして、フリーソフトであるにも関わらず、出力される楽譜は出版物並みに綺麗である。
MusiXTeXの書式は大変難しいが、そのプリプロセッサであるPMXM-Txを使うと、かなり容易である。M-Txは楽譜無いに歌詞も書けるようにPMXを拡張したものだが、若干PMXと書式が異なる。M-Txを使うにはPMXの知識も必要になるので、歌詞が要らないならPMXだけで充分と思う。

上の楽譜は、次のソースファイルから生成されたものである。

2 1 3 4 3 4 1 0
0 8 20 0.0

bt
./

It144ipipi
Ap

r4 rp ( a82 e+ a ) r r4 ( e82 e+ gs ) r r4 ( a82 e+ a ) r r4 /
( e85 ds e ds e b dn c a4 ) r8 ( c- e a b4 ) r8 ( e- gs b c4 ) r2 /


"./"までの部分は、PMXファイルとして必須の、楽譜全体に関する設定(ヘッダ)であり、"r4"より下の部分が各小節のデータ(ブロックの集まり)である。最低それらが含まれていれば、PMXのソースファイルとして成立する。
このファイルをforelise.pmxとして保存し、PMXがインストールされている環境で
pmx forelise.pmx

または、
pmxab forelise
musixflx forelise
musixtex forelise

とすると、TeXでお馴染みのDVIファイルができる。dvipsでPSファイルに、dvipdfでPDFファイルに変換できる。上の楽譜画像は、DVIファイルをdvipsでPSファイルにしてImageMagickでPNGに変換したものである。
1つ目のpmxコマンド一発でコンパイルする方法だと、ソースファイルに"I"コマンドがあると、ついでにMIDIファイルもできる。

上のソースファイルを解説する。
●ヘッダーの先頭から12個の数字(必須)


通称上での値意味
nv2パートの数
noinst1楽器の数
mtrnuml3(小節の境目の判定に使う)何分の何拍子の分子
mtrdenl4(小節の境目の判定に使う)何分の何拍子の分母
mtrnump3(実際の表示に使う)何分の何拍子の分子
mtrdenp4(実際の表示に使う)何分の何拍子の分母
xmtrnum01弱起(最初の小節が不完全)の場合の最初の小節の拍子数
isig0五線譜の左端の調を決める#や♭の数(負だと♭)
npages0出力ページ数、0は自動
nsyst8五線譜の段数(全ページ分)npages=0の時は1段当たりの小節数
musicsize20音符のサイズ
fracindent0.0五線譜の1段目のインデント量

●ヘッダーのそれ以降の行(必須)

項目上での値内容
楽器名 (noinst)行分の、楽器名(ここでは空行)
音符記号bt各パートの、ト音記号(t)かヘ音記号(b)かの設定(下のパートから)
出力先./出力先ディレクトリ("./"はカレントディレクトリ)

●楽譜全体の設定
"I"で始まる行はMIDIファイル出力の設定
 "t144"はテンポ
 "i...."は(nv)個分の楽器("pi"はピアノ)
"Ap"はMusiXTeXの「Type Kスラー」(またはPostscriptスラー)という綺麗なスラーを使う設定(要インストール)
●楽譜本体
・下のパートから書く
・パートの切り替えは"/"
・音符と音符の間はスペースで区切る
・音符について
 cdefgabはドレミファソラシ、rは休符
 音符内の数字は、音の長さ(1,2,4,8,16,32分音符はそれぞれ0,2,4,8,1,3)とオクターブ(4がト音記号の普通の所)
  長さを省略すると1つ前の音と同じ長さ
  オクターブを省略すると1つ前の音に一番近い高さ(c b(ド、シ)ならシは1つ下のオクターブのシ)
 +/-はオクターブの上下移動
 "s"は#(♭なら"f")
・()はタイ/スラー
・rpは全休符
・小節の区切りは自動的に計算される

もう1つサンプルを作ってみた。解説は省略する。
csikospost.pmx(PMXソースファイル)
csikospost.pdf(できた楽譜)
csikospost.mid(できた音楽ファイル)

ついでにもう1つ公開する。
Nukenin1.pmx(PMXソースファイル)
Nukenin1.pdf(できた楽譜)
Nukenin1.mid(できた音楽ファイル)

参考資料:マニュアル(英語、PDF)

viやEmacsで楽譜を入力する。TeXで楽譜を作る。何とエレガントな響きであることか。

See more ...

Posted at 22:57 in PC一般 | WriteBacks (0)
WriteBacks

Java3Dの光源と物体の表面属性を操作してみる

Java3Dで物体に色を付けるには、光源と物体の表面属性を設定する必要がある訳だが、光源色と表面属性の色成分の関係がよくわからない。光源色を変えても色が変わらなかったり、物体表面の反射成分を変えても色が変わらなかったりする。いろんなパラメーターをいじって試しても、物体が真っ白のままで何が悪いのかさっぱり見当がつかないと、全く理解が進まない。筆者はもう1年以上わけがわからないままである。

そこで、光源色と物体の表面属性をスライドバーで操作できるアプレットを作ってみた。

光源/表面属性のテストアプレット
ソースコード

点光源が上の方を旋回しており、平行光源は左手前にあって右奥に向けている。2つの小さな円錐は、光源のおおよその位置を表している。

操作できるパラメーターを説明する。

●光源について
・環境光源(Ambient light)
 空間全体に万遍なく存在する光。影を作らない。
・平行光源(Directional light)
 ある一定の方向から全空間を等しく照射する光。
・点光源(Point light)
 ある一点から放射状に放たれる光。光源からの距離に応じて減衰する。

●物体の表面属性(材質、Material)について
・拡散反射成分(Diffuse color)
 入射角に関係なく全方向に反射させる色成分。
・鏡面反射成分(Specular color)
 入射角と反射角の大きさが等しい反射をさせる色成分。
・環境光反射成分(Ambient color)
 環境光を反射させる色成分。
・発光色(Emissive color)
 光源に関係なく、自ら発する光の色成分。
・つや(Shininess)
 鏡面反射の鋭さ。

大体、自然界と同じように、光源色の一部の成分が反射光となる感じである。白色光に対しては、赤色光のみを反射する物体は赤くなるし、緑色光のみを反射する物体は緑色になるが、青色光に対しては、赤色光や緑色光のみを反射する物体は光を反射せず黒くなる。

バーをスライドさせても物体が真っ白のまま変化しなくても、根気よくいろんなバーをスライドさせていれば、どれかを触った時にきっと色に変化が起こる。
どうやら有効なパラメーターはTarget colorによって変わるようであるが、まだ関係がよくわからない。そのうち、ちゃんと理屈を知ってからいじって理解しようと思う。今回はこれを作っただけで満足する。

今回のアプレットでは、Java3Dのシーングラフ(オブジェクトの木構造)を途中で変更している("Change shape"ボタン押下時)。それに関して調べたことを記録する。
・J3Dスレッドが使用中のシーングラフから削除できるノードはBranchGroupノードのみである。従って、後で切り離したいノードは切り離し用のBranchGroupノードの下に配置する必要がある。
・シーングラフを操作するためには、一度Universeから切り離す必要がある。そのため、シーングラフにはALLOW_DETACHを設定しておく。
・SimpleUniverseの場合、シーングラフのdetachは

SimpleUniverse#getLocale()#removeBranchGraph()

で行う。その後のattachは通常のSimpleUniverse#addBranchGraph()で良い。
・BranchGroupノードに繋がる子ノードを追加したり削除したりするためには、そのBranchGroupノードにALLOW_CHILDREN_WRITEを設定する必要がある。
・AppearanceやMaterialのインスタンスは、シーングラフをdetachせずに変更してもJ3Dの描画に反映されるようである。

See more ...

Posted at 22:57 in Java | WriteBacks (0)
WriteBacks

Java3Dの光源と物体の表面属性を操作してみる

Java3Dで物体に色を付けるには、光源と物体の表面属性を設定する必要がある訳だが、光源色と表面属性の色成分の関係がよくわからない。光源色を変えても色が変わらなかったり、物体表面の反射成分を変えても色が変わらなかったりする。いろんなパラメーターをいじって試しても、物体が真っ白のままで何が悪いのかさっぱり見当がつかないと、全く理解が進まない。筆者はもう1年以上わけがわからないままである。

そこで、光源色と物体の表面属性をスライドバーで操作できるアプレットを作ってみた。

光源/表面属性のテストアプレット
ソースコード

点光源が上の方を旋回しており、平行光源は左手前にあって右奥に向けている。2つの小さな円錐は、光源のおおよその位置を表している。

操作できるパラメーターを説明する。

●光源について
・環境光源(Ambient light)
 空間全体に万遍なく存在する光。影を作らない。
・平行光源(Directional light)
 ある一定の方向から全空間を等しく照射する光。
・点光源(Point light)
 ある一点から放射状に放たれる光。光源からの距離に応じて減衰する。

●物体の表面属性(材質、Material)について
・拡散反射成分(Diffuse color)
 入射角に関係なく全方向に反射させる色成分。
・鏡面反射成分(Specular color)
 入射角と反射角の大きさが等しい反射をさせる色成分。
・環境光反射成分(Ambient color)
 環境光を反射させる色成分。
・発光色(Emissive color)
 光源に関係なく、自ら発する光の色成分。
・つや(Shininess)
 鏡面反射の鋭さ。

大体、自然界と同じように、光源色の一部の成分が反射光となる感じである。白色光に対しては、赤色光のみを反射する物体は赤くなるし、緑色光のみを反射する物体は緑色になるが、青色光に対しては、赤色光や緑色光のみを反射する物体は光を反射せず黒くなる。

バーをスライドさせても物体が真っ白のまま変化しなくても、根気よくいろんなバーをスライドさせていれば、どれかを触った時にきっと色に変化が起こる。
どうやら有効なパラメーターはTarget colorによって変わるようであるが、まだ関係がよくわからない。そのうち、ちゃんと理屈を知ってからいじって理解しようと思う。今回はこれを作っただけで満足する。

今回のアプレットでは、Java3Dのシーングラフ(オブジェクトの木構造)を途中で変更している("Change shape"ボタン押下時)。それに関して調べたことを記録する。
・J3Dスレッドが使用中のシーングラフから削除できるノードはBranchGroupノードのみである。従って、後で切り離したいノードは切り離し用のBranchGroupノードの下に配置する必要がある。
・シーングラフを操作するためには、一度Universeから切り離す必要がある。そのため、シーングラフにはALLOW_DETACHを設定しておく。
・SimpleUniverseの場合、シーングラフのdetachは

SimpleUniverse#getLocale()#removeBranchGraph()

で行う。その後のattachは通常のSimpleUniverse#addBranchGraph()で良い。
・BranchGroupノードに繋がる子ノードを追加したり削除したりするためには、そのBranchGroupノードにALLOW_CHILDREN_WRITEを設定する必要がある。
・AppearanceやMaterialのインスタンスは、シーングラフをdetachせずに変更してもJ3Dの描画に反映されるようである。

See more ...

Posted at 22:38 in Java | WriteBacks (0)
WriteBacks