Dec 31, 2009

ペンローズ・タイル

ペンローズのタイルというものがある。鋭角72°、鈍角108°の同じ大きさの菱形を

このように分けてできる2種類のタイルで、同じパターンを繰り返すことなく、平面を無限に隙間なく敷き詰められるというものである。繰り返し無しというのがポイントで、例えば2種類のタイルを上の図のように合わせると菱形であるから、その菱形を繰り返し並べると無限に敷き詰められるのは当たり前である。そうではなくて、どの一部分を取っても、全体がそれを繰り返し並べたものではないように、平面を敷き詰めることが可能なのである。

10年以上前に、それを何かで読んで、ボール紙をはさみで切って試したことがある。いつかプログラムを作って自動的に埋めさせてみたいと思って、ついこの間までそのボール紙を取ってあったが、実現しないまま、もういいかと思って捨ててしまった。そしたら、今月、某TV番組でペンローズ・タイルが紹介されたのを見て、過去の思い出が蘇り、日々そのことが気になるようになってしまった。
頂点に並ぶタイルは最大5枚だから、単純に探索する順をランダムにして5手先まで縦型探索するのを繰り返せば良いのではないか…と思って、意図的かつ不可避に1人で過ごすことを選択したクリスマスイブの前日にその選択が正解だったことにするためにトライしたら、うまくいかなかった。やっぱり何十手も先まで読まないとNGとわからない置き方があるのと、それ以前に重なり判定が難しすぎるのである。(実際にイメージバッファに描きながら次に描く所が未使用かどうかを調べる方法でも難しい)

年を越す前になんとか一悶着、いや一段落を…と思ってインターネットを読み漁ったら、プログラむのにうってつけの簡単な方法が見つかった。
Conway's concept of "inflation" and "deflation"

"inflation"は少し面倒なので、"deflation"と切り抜きで繰り返せるようにしてみた。
"deflation"によるペンローズ・タイリングのJavaアプレットの起動ページ
ソースコード
一応、"deflate"と"zoom"を繰り返すと、無限にパターンが作れるはずである。

See more ...

Posted at 13:43 in Java | WriteBacks (1)
WriteBacks

今年(2011)のノーベル化学賞ですね。ペンローズ・タイルそして準結晶構造、美しいです。ありがとうございました。

Posted by kyu at 10/14/2011 04:56:39 PM

Nov 29, 2009

余弦定理、正弦定理の導き方

今週、TVで「余弦定理」という言葉を耳にした。高校でそういう公式を習って名前だけは覚えていたが、中身は雰囲気すら思い出せなかった。
15年以上前にセンター試験を受けた時は完璧に暗記していたはずである。三角関数の加法定理は今でも身体が覚えていた。漢字と同じように書き順で覚えていて、右手が紙に書き始めると、何も考えなくても勝手に手が動いて、筆記体のような加法定理を書き終える。ちなみに我が筆記体では+とxは同じ形であった。
倍角の公式や和積・積和の公式や合成の公式(sinθ+cosθ=sin(θ+α)という変形)を加法定理から無事に導き出せて、正弦定理も何とか思い出すことができたのだが、余弦定理は影も形も思い出せなかった。
仕方なく、インターネットで調べた。

●余弦定理

三角形の3辺の長さをa,b,c、各辺の向かい側の角度をA,B,Cとすると、
a^2=b^2+c^2-2bc¥cos A
b^2=c^2+a^2-2ca¥cos B
c^2=a^2+b^2-2ab¥cos C
が成り立つ。

●正弦定理
同じく、
¥frac{a}{¥sin A}=¥frac{b}{¥sin B}=¥frac{c}{¥sin C}=2R
(Rは外接円の半径)が成り立つ。

そうそう、そういえばそんなのがあった。何でそういうのが成り立つのかがわからないまま、丸暗記した記憶がある。何でそうなるんだろう?と、何度も不思議に思った記憶もある。
これらを初めて教わってから20年、ついに導き方を理解した。

●余弦定理の導出

このようにaに垂線を引くと、
a=b¥cos C+c ¥cos B
の関係があることがわかる。これの両辺を2乗すると、
a^2=b^2¥cos^2 C+c^2¥cos^2 B+2bc¥cos B¥cos C
である。sin2+cos2=1やコサインの加法定理
¥cos(B+C)=¥cos B¥cos C-¥sin B¥sin C
を使って整理すると、

となる。A+B+C=π(180°)だから、
¥cos(B+C)=¥cos(¥pi-A)=-¥cos A
であり、b sinC = c sinB = aへの垂線の長さで末尾の2乗部分は0なので、余弦定理の式が求まる。

●正弦定理の導出

このように同じ外接円を持つ1辺がaの直角三角形を描くと、円周角の定理より、aの向かいの角の角度はやはりAである。また、同じく円周角の定理より、円周角が90°なら中心角が180°なので、斜辺は外接円の中心を通り、長さは直径=2Rに等しい。従って、2R sinA=aである。

See more ...

Posted at 11:23 in 数学 | WriteBacks (2)
WriteBacks

 面白い!定義ばかりで覚えていたのでは直ぐ忘れます。

Posted by 村山良平 at 02/17/2011 03:58:27 PM

ご同意ありがとうございます。

Posted by ynomura at 02/18/2011 12:36:31 AM

Nov 22, 2009

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 12:52 in PC一般 | WriteBacks (0)
WriteBacks

Oct 03, 2009

統計学復習メモ17: 重回帰分析の計算結果の検定

以前に回帰直線を求めることと主成分分析は似ていると書いたが、ちょっと無理があり過ぎたかも知れない。なぜか主成分分析が因子分析より知名度が低いのが気に入らなくて、メジャーなものと結びつけたくて書いたのだが、やっぱり回帰分析と主成分分析は全然違う。回帰分析は目的変数を説明変数で説明しようとするものであり、主成分分析はデータの傾向を見るものである。通常、回帰直線を求める際に最小にする誤差は目的変数軸上の誤差であるが、第一主成分を求める為に最小にする誤差は全変数が構成する空間での誤差である。
筆者は、主成分分析は多変量解析の基本かつ最も重要であると考えており、学生時代には趣味のプログラミングのネタにするほど惚れ込んだが、因子分析はあまり好みでないのである。暴走ついでにさらに暴走すると、筆者の中では、

A群B群C群
因子分析主成分分析回帰分析
LinuxFreeBSDOpenBSD
RubyPythonPerl
WindowsMacUNIX
SleipnirLunascapeFirefox
秀丸サクラエディタEmacs
のような関係があり、因子分析は軟派で軽薄でチャラいのである。


閑話休題。今回は重回帰分析を復習する。一次回帰直線が、xを説明変数、yを目的変数とし、直線y^=a x+bの係数aと切片bをデータから求めたものだとすると、重回帰分析はそのxを複数次元にしたもの、
¥hat{y}=a_1x_1+a_2x_2+¥cdots+a_mx_m+b
の係数aiとy切片bを求めるものである。
便宜上、上のbをa0と書き、
¥hat{y}=XA=¥pmatrix{1 & x_1 & ¥cdots & x_m}¥pmatrix{a_0 ¥cr a_1 ¥cr ¥vdots ¥cr a_m}
と行列の積で表すことにする。

n個のデータの組(yi, x1i, x2i, x3i, ...)(1<=i<=m)が得られた時、誤差をeiとすると、

である。以下、これをY=XA+Eと書く。
yの誤差eiの2乗和が最小になるようにAを求めるとすると、
¥sum_{i=1}^{n}e_i^2=E^TE=(Y-XA)^T(Y-XA)
であり、これはaiの2次関数なので、これを偏微分して0になるAが答であるから、
¥frac{¥partial}{¥partial A}(E^TE)=¥frac{¥partial}{¥partial A}(Y^TY-2(XA)^TY-(XA)^TXA)=-2X^TY+2X^TXA=0
なので
A=(X^TX)^{-1}X^TY
という計算で求められることがわかる。(XTX)が正則でないと逆行列が存在しないという問題はあるが、異なるデータの個数が説明変数の数より少ないとか、中身が0だらけのぜろりとしたデータでない限りは正則だろうから、とりあえず気にしないことにする。

試しに次のサンプルデータのバグ数について計算してみる。

部品バグ数開発規模関係者数疲労度
A70.935
B152.118
C101.249
D71.217
E60.359
F50.693
G70.922
H50.363
I60.358
J70.634
これを"bugs.dat"というタブ区切りのテキストファイルにして、Maximaで
dat:read_matrix("bugs.dat");
n:10$ /*標本数*/
p:3$ /*説明変数の数*/
Y:submatrix(dat,2,3,4); /*1列目を取り出す(2-4列目を除去)*/
X1:submatrix(dat,1); /*2列目以降を取り出す(1列目を除去)*/
X:map(lambda([x],append([1],x)),X1); /*X1の左端に要素が全て1の列を追加*/
A:invert(transpose(X).X).transpose(X).Y; /*回帰係数*/
fpprintprec:4$
float(A);
を実行すると、次の結果が得られる。
a02.1
a14.612
a20.02896
a30.2436
バグ数=2.1 + 4.612*開発規模 + 0.02896*関係者数 + 0.2436*疲労度、という意味である。
この結果がどれくらい元データに合ってるかを、Yの値で見てみる。
/*元データ、回帰式による推定値、誤差*/
[Y, X.A, Y-X.A];
元データYYの推定値Y^Y-Y^
77.556-.5558
1513.761.237
109.943.05743
79.369-2.369
65.8210.179
55.859-.8589
76.7960.204
54.388.6115
65.577.4226
75.9291.071

YとY^が十分近いので、これでも満足なのだが、統計学にはこの回帰分析結果をきちんと検定する方法が存在するようである。
その前に一応、よくある指標を計算する。以下、データ数をn、説明変数の数をpとする。

load(descriptive)
load(distrib)
n:10$
p:3$
/*重相関*/
R:cor(addcol(Y,X.A)),numer;
/*決定係数(寄与率)*/
R^2;
/*自由度補正済み決定係数*/
1-(1-R[1,2]^2)*(n-1)/(n-p-1);
重相関0.9361
決定係数0.8762
自由度修正決定係数0.8143
相関係数の類は、値の解釈が難しいので、相関係数同士を比較する以外には使いづらい。一応、自由度修正済み決定係数が0.8以上というのがよく使われる基準であるらしい。
これは置いといて、まず、回帰式全体の有効性を検定する。これには、回帰分散÷残差分散、
¥frac{V_{¥hat{y}}}{V_¥varepsilon}=¥frac{¥frac{1}{p}¥sum_{i=1}^{n}(¥hat{y}_i-¥bar{y})^2}{¥frac{1}{n-p-1}¥sum_{i=1}^{n}(¥hat{y}_i-y_i)^2}
という値が自由度p,n-p-1のF分布に従う、というのが使えるらしい。YとY^の誤差の分散が十分に小さいかどうかを調べようとするものと言えるだろう。
E: Y-X.A$
/*回帰分散*/
Y_hat: X.A-mean(Y)[1]$
Vy: transpose(Y_hat).Y_hat/p,numer;
/*残差分散*/
Ve: transpose(E).E/(n-p-1);
/*回帰方程式のF値*/
f:Vy/Ve;
/*そのF値以上の発生確率*/
1-cdf_f(f,p,n-p-1);
結果
F値14.16
P(F(3,6)>14.16)0.03951
有意水準を0.05とすると、回帰分散=残差分散という仮説は棄却されるので、この回帰式は有効ということになる。

次に、回帰係数aiそれぞれの有効性を検定する。これには、回帰係数÷標準誤差=
¥frac{a_i}{¥sqrt{s^{ii}V_{¥varepsilon}}}
(sijは共変動行列(=n×共分散行列)の逆行列の要素)
という値が自由度n-p-1のt分布に従う、というのを使う。定数項a0については、siiの部分は次のような式になるが、同じようにt値が求められる。
s^{00}=¥frac{1}{n}+¥sum_{i=1}^{p}¥sum_{j=1}^{p}¥bar{X_i}¥bar{X_j}s^{ij}

/*共変動行列の逆行列*/
Sinv: invert(cov(X1)*n);
/*標準誤差(係数項)*/
Se[i]:=sqrt(Sinv[i,i]*Ve);
/*標準誤差(定数項)*/
Se[0]:sqrt((1/n + sum(sum(mean(X1)[i]*mean(X1)[j]*Sinv[i,j],j,1,p),i,1,p))*Ve);
for i:0 thru 3 do print(Se[i])$
iaiの標準誤差
01.863
11.023
20.2303
30.1661
/*回帰係数のt値、|t|値がそれより大きい確率*/
t[i]:=A[i+1,1]/Se[i]$
for i:0 thru 3 do print(t[i], 2*(1-cdf_student_t(t[i],n-p-1)))$
it値確率
01.1280.3025
14.5070.004074
20.12570.904
31.4670.1928
有意水準を0.05とすると、a1についてのみ、本当は0である(係数に意味が無い)という仮説が棄却され、有効ということになる。

See more ...

Posted at 20:46 in 数学 | WriteBacks (2)
WriteBacks

どうして定数項の標準誤差はサンプル数の逆数に、各説明変数の積を共分散でに近い形で除することで求まるのですか? どうしても導出方法が知りたいです。 ぜひ教えてください。

Posted by 理系学生 at 10/14/2010 09:04:54 PM

コメント欄に画像が貼れないので、上に追記しました。

Posted by ynomura at 10/24/2010 08:46:51 PM

Sep 02, 2009

統計学復習メモ16: 帰無仮説と棄却域

統計学の検定というのは、「帰無仮説」と「対立仮説」を設定して、「帰無仮説」を「棄却」できれば「対立仮説」が「有意」であるとする、というものである。
…というように言葉で書かれるとそんなに難しくなさそうなのだが、ここまでで既に「帰無仮説」だの「有意」だの独特の用語が出てくる上、帰無仮説が棄却できない時は「対立仮説が有意であるとは言えない」だの「差が無いという仮説は棄却されない」だの「差があるとは言えない」だの回りくどい表現が出てきて、混乱する。
棄却して「無」に「帰」するための仮説だから「帰無仮説」と言うんだよ、なんて説明は大して理解の助けにならない。むしろ、それで何か解った気にさせて思考を停止させる、本質の理解を邪魔する雑音にすら思える。
なぜ「差が無いとは言えない」ではないのか。そういう疑問を持ってもう一度理解しようとすると、「第1種の過誤」や「第2種の過誤」や「有意水準」や「危険率」や「検出力」が出てきて、阻止されてしまう。命題の逆が偽(矛盾する)なら命題が真となるような単純な代物ではないのである。

筆者は大学で「検定」を習って15年以上、統計学を復習して1年以上、まだ混乱している。そのため、自己流に整理して、悪あがきしてみる。

仮説として、ある統計量が、下の図の緑色の線で示されるような確率分布に従うと仮定されているとし、標本から得られた統計量がtであるとする。横軸がt、縦軸が発生確率である。

上記の仮説が正しいとすると、t=0.7(青い点の辺り)やt=1.5(紫色の点の辺り)はまだそれなりに発生確率が高いが、t=2.2(赤い点の辺り)は発生確率が低い。従って、t=2.2となるような標本が得られたら、それはかなり稀な事が起こったという事になる。その時、そうではなくて仮説が誤っているのだと考えるのが統計学における検定の考え方である。

そのために、どれくらい稀だと仮説が誤っているとする(仮説を棄却する)かを決める。その確率が有意水準である。何故だかわからないが、0.05や0.01がよく使われるようである。ここでは有意水準を0.05(5%)とする。
次に、統計量がどの範囲だと仮説を棄却するか(棄却域)を、その範囲の発生確率が有意水準に等しくなるように決める。例えば、上図の緑色で塗った部分は、0から遠い、発生確率が全体の5%の範囲である。すなわち、標本のtがこの範囲(大体t<2とt>2)に入れば、仮説を棄却する、という設定である。

注意すべきは、仮説を棄却できなかったからといって仮説が正しいとは言えないことである。仮説が正しいとすると稀にしか起こらない事が起こったかどうか、を調べているだけなので、わかるのは、仮説が誤っているか、何とも言えないか、のどちらかである。調べている仮説を帰無仮説、その反対の仮説を対立仮説として言い換えると、わかるのは、帰無仮説が誤っていて対立仮説が有意であるか、帰無仮説が誤っているとは言えず対立仮説については何とも言えないか、のどちらかである。

上記の設定では、もし標本から得られたtが青や紫の位置の値なら、仮説が棄却できず(検定失敗)、いわゆる玉虫色の決着となるが、もしtが赤い位置の値なら、仮説が間違っているだろうと言える。

以上で検定の計算は可能だが、設定した検定方法が妥当かどうか、また検定結果の意味を考えるためには、他に考える事がある。
帰無仮説にて統計量の確率分布を仮定するということは、標本の統計量にもばらつきがあることを仮定している。もし標本の統計量にばらつきが無いとするなら、棄却域に入るかどうかという以前に、ばらつきがあるとする帰無仮説が誤りになってしまう。
例えばもし、標本の統計量が帰無仮説の通りでないのに、平均的には棄却域に入らない場合は、次のような感じになる。

赤い点が統計量の平均、紫色の線が分布である(見やすさのため、縦に縮めている)。ピンクで塗られている部分が、帰無仮説(緑線)の棄却域に入る部分である。
もし帰無仮説が正しい場合は、紫の線は緑の線に一致し、ピンクの領域は有意水準と同じ5%となる。そのため、帰無仮説が正しくても、5%の確率で帰無仮説が棄却されてしまう(第一種の誤り)。帰無仮説は棄却するために設定する(誤りだろうと思うものを帰無仮説に採る)ものだが、棄却できたと喜んでも、誤って棄却した可能性が残るのである。有意水準は、この第一種の誤りを犯す確率という意味で、危険率とも呼ばれる。

上の図では、緑の線と紫の線はずれており、帰無仮説は誤りである。しかし、ピンクの領域、すなわち帰無仮説が棄却される確率は12.5%しかない。「検出力」が12.5%しか無い、と言う。
次の図は、統計量の平均が棄却域に入るくらいにずれている例である。

ピンクで塗っていないが、ここまでずれると検出力は70%に達する。しかし、逆に言うと、ここまで帰無仮説が誤っていても、30%の確率で棄却されないのである(第二種の誤り)。有意水準を下げて第一種の誤りを犯す確率(緑の領域)を小さくすると、第二種の誤りを犯す確率(藤色の領域)が大きくなる、という関係があることがわかる。

次に、帰無仮説と棄却域の取り方について考える。
通常、帰無仮説は「差が無い」方に取られる。これは、統計量が0から遠い方が「差がある」ことが多いというか、0から遠い方が「差がある」ように統計量を定める方が便利であることが多いからである。例えば標準正規分布やt分布やカイ2乗分布に従う統計量は0から遠いと平均と差がある。
そのため、もし帰無仮説を「差がある」方に取ると、統計量の棄却域が0に近い方になり、範囲が狭くなる。標準正規分布やt分布のような山形の分布だと、

このように棄却域が狭くなり、統計量が繊細で、計算上不便である。

それよりも、実際の検定統計量は和や積を使って計算されるため、「差がある」は単純ではなく、棄却するのは難しいのである。例えば、和で効いてくる平均値の場合、実際には標本に帰無仮説と全然違う大きなばらつきがあっても、相殺されて帰無仮説の平均値に近くなってしまうことがあり得るし、積で効いてくる場合は標本のどれか1つの確率が(帰無仮説の前提から外れて)0に近いために結果として統計量が0に近くなってしまうことがあり得る。「差がある」を棄却するには色々な可能性を考えないといけないが、統計量が0から遠ければそれだけで「差が無い」を棄却できるので、やりやすいのである。

そのため、上の例のように統計量が標準正規分布やt分布などの0中心の左右対称の確率分布に従うと仮定する場合は、棄却域を0から遠い左右両端に定めることが多い(両側検定という)。
カイ2乗分布に従う統計量のように、正の値しか取らない場合は、0から遠い方の(+∞までの)1つの領域を棄却域とする(片側検定という)場合が多いが、それでも、どこかに山がある場合は0に近い方の裾と0から遠い方の裾の両側に棄却域を設ける場合もあるようである。

次の図は、カイ2乗検定(片側検定)の棄却域の例である。

緑の線が、統計量が自由度2のカイ2乗分布に従うという帰無仮説であり、緑で塗られた領域が有意水準(αと書かれることが多い)5%の棄却域である。赤い点のように、標本から得られた統計量がここに入れば、緑の線を棄却する。
実際には自由度4のカイ2乗分布に従う場合は、帰無仮説と標本は次のような関係になる。

藤色の部分が第二種の誤りを犯す範囲(βと書かれることが多い)で、確率は約0.8であり、ピンク色の部分の大きさが検出力であり、(1-β)≒0.2である。

See more ...

Posted at 22:32 in 数学 | WriteBacks (0)
WriteBacks

Aug 29, 2009

(FreeBSD)かなキー入力をするための設定

FreeBSD 7.1を使っていて、どうしても(ローマ字入力でない)かな入力をしたくなったので、設定方法を調べてみた。
以下、キーボードは日本語106キーボードを使用することを前提にする。

●X Windowの設定
まず、デフォルトの設定では、バックスラッシュキー(「ろ」キー)と¥キー(「ー」キー)のkeysymが同じでアプリの設定ではどうしようも無いため、「ろ」キーが"_"になるよう、keysymを変更する。また、Shift+0キーがShift+^キーと同じチルダ記号だと「を」キーとしての設定がうまくいかないことがあるので、Shift+0が"|"になるように変更する。

~/.Xmodmapに以下の2行を加える。

keycode 211 = underscore underscore
keycode 19 = 0 bar kana_WA kana_WO

それをすぐに反映するには、次のコマンドを実行する。

xmodmap ~/.Xmodmap

おかしい時は、xmodmap -pkeとして、現在の状況を確認する。

●Emacs+Anthyの場合
まず、anthy-conf.el(筆者の環境では/usr/local/share/emacs/site-lisp/anthy/にある)の「ろ」「を」「ー」キーの割り当てを変更する。ついでに、Shift+"["、Shift+"]"が「」(かぎ括弧)でなく{}(波括弧)になっているが、「」の方が使う機会が多いので、Shift+"["、Shift+"]"、を「」に割り当てる。


("_" . "ろ")
("|" . "を")
;;("_" . "ー")は削除
("\\" . "ー")
("{" . "「") ("}" . "」")

ひらがな用とカタカナ用の設定があるので、カタカナ用の設定も同じように変更する。

次に、~/.emacs.elに次の3行を加えるか、Emacs上でこれを実行する。

(load-library "anthy")
(anthy-kana-map-mode)
(setq default-input-method "japanese-anthy")

これでC-\でかな入力ができるようになる。

日本語入力切り替えキーとしてC-\が使えない時は、Canna流にC-oを割り当てる。

(global-set-key "\C-o" 'anthy-mode)

●Emacs+Cannaの場合
(1) default.canna(筆者の環境では/usr/local/share/canna/)を~/.cannaにコピー
(2) ~/.cannaの"default.cbp"を"kana.cbp"に変更
(3) ~/.cannaのdefsymbolの行を全てコメントアウトし、以下を追加

(defsymbol
?_ "ろ" "ろ" )

(4) ~/.emacs.elに(canna)を加えて、emcwsコマンドでEmacsを起動
または、
(4') emcwsコマンドでEmacsを起動し、M-x canna
これでC-oでかな入力ができるようになる。

●jvim+Cannaの場合
jvimは、DIRECT_CANNAを有効にしたもの(ONEWを無効にしたもの)をインストールすること。
packagesなら ja-jvim-direct_canna-3.0.j2.1a_2 など、"direct_canna"が付いたものを選択する。portsならjvimのmake configでDIRECT_CANNAを選んでからコンパイルする。

上記のEmacs+Cannaの場合と同様に~/.cannaを設定してjvimを起動すると、インサートモードにてC-\でかな入力ができるようになる。

以上、FreeBSDには限らない設定だとは思うが、他の環境では試していない。

See more ...

Posted at 22:56 in UNIX | WriteBacks (0)
WriteBacks

Aug 08, 2009

8方向に掘り進めて迷路を作ってみる

昔、MSX-FANという雑誌に投稿されていたBASICのプログラムで、壁を8方向に掘り進めて迷路を作り、その一部分が画面に表示された状態でその中をさまようというゲームがあった。

ランダムな迷路を作る方法としては、格子状の壁を4方向に掘り進める方法はすぐに思いつく。既に掘り進めた通路のいずれかの壁をランダムに選び、その壁を取って通路が1マスだけ伸びるならその壁を取り払う(壁の向こうに既に通路があるとループができてしまうため)というのを繰り返せば良いのである。選んだ壁の向こうが未踏の地ならその壁を外す、というのを繰り返すとも言える。
他に、通路が4方向に進む迷路を作る方法としては、「棒倒し法」なる奇妙な方法があるようだが、いたずらに複雑な上にいびつな迷路になるので、取るに足りない。

しかし、雑誌に載っていたプログラムは、迷路を8方向に掘り進めた。そういうことができること自体が私にとっては不思議だった上、そのプログラムは今考えても戦慄が走る「1画面プログラム」すなわち40文字x23行に収まる激短プログラムだった。

プログラムのコードは見えていたのである。大概のプログラムは、手で打ち込んでいる内にどういう仕組みかは理解できた。理解できないのはマシン語が混ざっているものくらいだった。ALL BASICで手で打ち込んで実際に動かしてみてアルゴリズムが理解できなかったプログラムは、そのプログラムくらいしか記憶に無い。
加えて、その作者は当時の私と同じ奈良県民らしかったため、忘れられなかった。

それから約20年、そのプログラムのことを時々思い出しては、頭の中で考えて、その迷路を作る仕組みを不思議に思っていた。

昨日、またふとそのプログラムのことを思い出して、図を描いて考えてみたら、3分でわかった。

テストプログラム(Javaアプレット)の起動ページへ
ソースコード

ついでに自動的に解かせてみた。

要するに、隣り合う斜めの通路が接しないような幅で通路を掘れば良いのである。他に気を付けるのは、斜めのマスに通路を伸ばす時は直交する通路が無いことを確認することくらいで、基本的には4方向に掘り進める場合と同じように方眼のマス目を繋げていけば実現できる。

See more ...

Posted at 23:24 in Java | WriteBacks (0)
WriteBacks

Aug 02, 2009

正20面体の頂点座標の求め方

正20面体の頂点座標を簡単に求める方法があることがわかった。


同じ長方形3枚をこのように組み合わせて、近い頂点を結ぶと、何か20面体ができる。これの長方形の長辺と短辺の比を変えていくと、2つの長方形から3点取った三角形がどこかで正三角形になる。

その時、20面体の30本の辺の長さが全て同じになるので、その時の3つの長方形の頂点が正20面体の頂点である。

長方形の短辺の長さを2、長辺の長さを2xとし、点A〜Eを

このように取ると、AE=x、ED=(x-1)なので、AD2=x2+(x-1)2であり、DC=1なので、AC2=AD2+DC2=x2+(x-1)2+1である。ABCが正三角形の時、AC=2なので、
x2+(x-1)2+1=4
が成り立ち、これを解くとx=(1+√5)/2(=黄金比)が求まる。

従って、Gを黄金比とすると、
(±1,±G,0)
(0,±1,±G)
(±G,0,±1)
の12点が、1つの正20面体の全ての頂点となる。

See more ...

Posted at 19:27 in 数学 | WriteBacks (0)
WriteBacks

Jul 19, 2009

秘書問題の考察(4) 最小期待順位の極限値について

前回のエントリーで説明した、秘書問題(secretary problem)の期待順位を最小にする解に関して、もう少し追究する。

まず、前回は参考文献を鵜呑みにした、j人中x位の人がn人中何位であるかの期待値が
¥sum_{y=x}^{n-j+x}y¥frac{¥pmatrix{n-y¥cr j-x}¥,¥pmatrix{y-1¥cr x-1}}{¥pmatrix{n¥cr j}} = ¥frac{n+1}{j+1}¥,x
となること(左辺が右辺に等しくなること)を導出する。
これの左辺をy(j,x,n)と書く。j人中x位の人は、もう1人加えて(j+1)人とすると、x位か(x+1)位かのどちらかになる。(j+1)人目が上位x位に入ると(x+1)位、入らないとx位のままである。このことから、j人中x位の人の最終順位は(j+1)人中(x+1)位の人の最終順位と(j+1)人中x位の人の最終順位で表され、
y(j,x,n) = x/(j+1) y(j+1,x+1,n) + (j+1-x)/(j+1) y(j+1,x,n)
の関係が成り立つことがわかる。扱いやすいように、jを1つずらして、
y(j-1,x,n)=¥frac{x}{j}¥, y(j,x+1,n) + ¥frac{j-x}{j} y(j,x,n)
としておく。
y(n,x,n)=xは自明だから、そこからy(n-1,x,n)が求まり、再帰的に計算するとy(j,x,n)が求まる。
j=n-1については、
y(n-1,x,n) = ¥frac{x}{n}(x+1) + ¥frac{n-x}{n} x = ¥frac{n+1}{n}x
j=n-2については、
y(n-2,x,n) = ¥frac{x}{n-1}¥frac{n+1}{n}(x+1) + ¥frac{n-1-x}{n-1}¥frac{n+1}{n}x = ¥frac{n+1}{(n-1)n}¥left(x(x+1)+(n-1-x)x¥right)=¥frac{n+1}{(n-1)n}nx=¥frac{n+1}{n-1}x
なので、y(j,x,n) = (n+1)/(j+1) xではないかと予想できる。それを帰納法で確認する。
j=n-1, n-2については上のように成り立つ。あるjについて成り立つとすると、
y(j-1,x,n) = ¥frac{x}{j}¥frac{n+1}{j+1}(x+1) + ¥frac{j-x}{j}¥frac{n+1}{j+1}x = ¥frac{n+1}{j(j+1)}¥left(x(x+1)+(j-x)x¥right)=¥frac{n+1}{j(j+1)}(j+1)x=¥frac{n+1}{j}x
となり、(j-1)についても成り立つので、間違いなくy(j,x,n) = (n+1)/(j+1) xである。

次に、Chowという人の1964年の論文によるらしい、その筋では頻出の、筆者を秘書問題地獄に陥れたミラクルな定理を紹介する。
定理:期待順位最小の最適戦略による期待順位をV(n)とすると、
¥lim_{n¥rightarrow¥infty}V(n)=¥prod_{j=1}^{¥infty}¥left(¥frac{j+2}{j}¥right)^¥frac{1}{j+1}¥simeq3.8695
である。
つまり、応募者が何人に増えようが、千人だろうが1万人だろうが、平均的にその中の3.87位以上の人を選択できる採用戦略が存在するのである。その数字の小ささもさることながら、その数字がnを含まないある定数に収束するというのは驚きではなかろうか。

さらに、期待順位最小の最適戦略を、r(x,n)番目の応募者からは暫定x位以上の人を採用すると書くと、次の式が成り立つ。
¥lim_{n¥rightarrow¥infty}¥frac{r(x,n)}{n}=¥frac{1}{V(¥infty)} ¥prod_{j=1}^{x-1}¥left(¥frac{j+2}{j}¥right)^¥frac{1}{j+1}
r(x,n)/nは面接した応募者の割合だから、nが大きい時、全応募者のどれくらいの割合を面接した時点から何位以上の人を選べば良いかというのが、これで計算できる。x=1〜10について計算すると


xr(x,n)/n
10.2584
20.4476
30.564
40.6408
50.6949
60.735
70.7658
80.7903
90.8101
100.8265

となるので、最初の25.8%の人はパスし、25.8%の人が過ぎたら1位の人を、44.7%の人が過ぎたら2位以上の人を、56.4%の人が過ぎたら3位以上の人を(以下略)採用すれば良いということになる。(但し、これはnが小さい場合は正確な値との差が大きい。前回のエントリーでの計算結果と比較すると、n=20でもかなりずれている。)

参考文献
Optimal Stopping and Applications(最適停止理論の大家、Thomas S. Ferguson教授による解説)Chapter 2.

See more ...

Posted at 16:54 in 数学 | WriteBacks (0)
WriteBacks

Jul 12, 2009

ソフトウェアリリース問題

秘書問題」の参考文献として図書館から借りた本に、曲がりなりにはソフトウェアに関わっている筆者としては、読んどかないと気になって仕方なくなるに違いない、一際目を引く問題があったので、返却する前に一読した。
「ソフトウェアリリース問題」と呼ばれている、最適停止問題の1つだ。
単純化した概要をメモする。

問1:テストが未実施である、開発中のソフトウェアの全バグ数がBとし、n回目のテストで見つかるバグ数をX(n)とする。BとX(n)の確率分布は既知で、E(X(n))は単調減少とする。見つかったバグは、各回のテスト終了時にコスト0で直ちに直される(バグ修正のコストはテストのコストに含まれる)とする。テスト1回当たりの費用をCt、ソフトウェアリリース後に残されるバグ1件当たりのコストをCb、Ct<<Cbとする時、コストを最小にするには、何回テストしてリリースするのが最善か。

ソフトウェアの開発の仕事をしたことがある人なら1回は考えさせられたことがあるであろう、どこまでテストすればいいのかという問題である。「バグが0になるまでに決まってるだろ!」と思う人が居るかも知れないし、そのように上から言われてる人が居るかも知れない。筆者は過去にそんな理不尽な圧力を何度も目撃した。しかし、現実的にはバグが0であることを証明するのは不可能である。ごくごく一部の例外を除き、はるかに不可能であり、不可能オブ不可能ズであり、1秒でも考えるのが無駄なくらい十分に不可能なのである。バグが無いことが証明されてリリースされるソフトウェアなんて存在しないのであり、どれだけテストして修正しても、ソフトウェアにバグは必ず存在するのである。
現実社会できちんと真っ当なソフトウェア開発の仕事をする人にとっては、バグの数が何らかの確率モデルで仮定されることは至極当然であり、むしろ、そうしないと仕事にならない。

解1:n回テストした時のトータルコストC(n)はn*Ct + (B - ΣX(k))Cbであり、
C(n+1)-C(n) = Ct - X(n+1)Cb
であり、E(X(n+1))は単調減少なので、E(C(n+1)-C(n))は単調増加であり、どこかでテストのコストの平均が残バグのコストの平均を上回る。E(C(n+1)-C(n)) > 0を解くと、E(X(n+1)) < Ct/Cbなので、答えはE(X(n+1)) < Ct/Cbとなるnである。

問2:Bが平均λのポアソン分布に従い、X(n)がその時点のバグ数と確率pをパラメーターとする二項分布に従う時、問1の最適なnはいくらか。

ソフトウェアのバグ数は、ソフトウェアの規模に比例すると考えられることが多い。システムの規模が大きくなると単位コード量当たりの複雑さも増すので、比例するというのはおかしい、と考えるのも自然なことなのだが、現実にはソフトウェア開発の対価はコード量に比例して支払われることが多いため、ソフトウェアは冗長であってもなるべく単純に作られ、単位コード量当たりの複雑さは上がらないのである。コード量が減るように工夫して効率の良いコードを作ると逆に報酬が下がってしまうので、コードは最適化されずに単純で似たようなコードの積み重ねとなり、システムの複雑さがそのままコード量として表れるのである。
従って、ソフトウェアのバグ数が(規模×単位コード量当たりのバグ数)を平均とするポアソン分布に従うと仮定するのは、現実社会に居て違和感を感じない。
1回のシステムテストで見つかるバグの数は、1つ1つのバグについて確率pで見つかると考えると、バグ数×pを平均とする二項分布に従うと仮定するのは、至極妥当であろう。

解2:B(n)をn回目のテスト後のバグ数とすると、1回のテスト当たりに見つからずに残る割合の平均が(1-p)なので、B(n)は平均λ(1-p)nのポアソン分布に従う。X(n+1)の平均はB(n)*pなので、答えはE(X(n+1)) = λp(1-p)n < Ct/Cbとなるnである。

参考文献
「タイミングの数理 最適停止問題」穴太克則 著(朝倉書店)

See more ...

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

Jul 09, 2009

秘書問題の考察(3) 期待順位最小の後ろ向き帰納戦略

前回のエントリーの秘書問題をきちんと解く。

問:以下の条件で、結果としてなるべく順位の高い人が選べるようにするには、どのようにすれば良いか。
・n人の応募者から必ず1人を選ぶ。
・1人ずつ面接し、その場で採用するかどうかを伝える。採用したら、残りの人は面接しない。採用しなければ、その人は2度と選べない。
・応募者は、それまでに面接した人の中で、必ず順位付けされる。応募者の良し悪しは、その暫定順位でしか判断できない。

解:
まず、j番目の応募者の暫定順位(それまでのj人中の順位)がxの時、その応募者が全n人中の何位であるかの期待値を求める。全n人中の順位をyとすると、この期待値は、(順位×確率)の和なので、E(y)=ΣyP(y)のようになる。j人中x位の人がn人中y位である確率は、
1 ...... y ...... n
 ↓
1 ... x ... j
という対応を考えると、j人の組み合わせnCj通りの中で、yの左側(y-1)人の内(x-1)人がxの左側に入って、yの右側(n-y)人の内の(j-x)人がxの右側に入る確率なので、P(y)=(y-1)C(x-1) × (n-y)C(j-x) ÷ nCjである。これを使うと、yの期待値は
¥sum_{y=x}^{n-j+x}y¥frac{¥pmatrix{n-y¥cr j-x}¥,¥pmatrix{y-1¥cr x-1}}{¥pmatrix{n¥cr j}} = ¥frac{n+1}{j+1}¥,x
となるらしい。(参考文献(1)による。この変形は私にはできなかったが、このP(y)は「超幾何分布」の形をしてるので、その知識を使えば簡単な計算で求まるのだと思う。)このE(y)を、以後y(j,x,n)と書く。

次に、最後の応募者の場合から、具体的な手順を考える。j-1番目までの応募者を採用しない時の最終的な採用者の期待順位をR(j,n)とする。
n-1番目までの応募者を採用せず、n番目の応募者を面接する場合、R(n,n)=n番目の応募者の期待順位=1位〜n位の平均=(n+1)/2である。
n-1番目の応募者については、その人が暫定x位の時、その人の最終的な期待順位はy(n-1, x, n)であり、その人を選ばない時の期待順位はR(n,n)である。従って、y(n-1, x, n) ≦ R(n,n)ならその人を採用すれば良いということになる。n-1番目の人がある暫定x位である確率は1/(n-1)なので、min{a,b}をa,bの小さい方とすると、R(n-1,n) = 1/(n-1) * Σ min{ y(n-1, x, n), R(n,n) }と書ける。
これは、n-1番目の応募者に限らず、1≦j<nの全てのjについて同じである。すなわち、j番目の応募者が暫定x位の時、その人の最終的な期待順位はy(j, x, n)であり、その人を選ばない時の期待順位はR(j+1,n)である。従って、y(j, x, n) ≦ R(j+1,n)ならその人を採用すれば良いというのが答である。R(j,n) = 1/j Σ min{ y(j, x, n), R(n,n) }、yを展開すると
¥frac{1}{j}¥sum_{x=1}^{j}min¥left¥{¥frac{n+1}{j+1}x, R(j+1,n)¥right¥}
なので、R(j,n)は末端のR(n,n)から逆向きに帰納的に(backward inductionと言うらしい)計算できる。

n=10の場合、採用する人の期待順位を最小にするために、j番目の人が暫定何位以内ならその人を選べば良いかは、次のようになる。3列目は、その人まで達した時の最終的な期待順位である。


jthresholdR(j,10)
954.28
833.59
723.15
622.89
512.68
412.56
302.56

3人目まではパスし、5人目までは暫定1位なら、7人目までは暫定2位なら、8人目は暫定3位なら、9人目は暫定5位なら採用すれば良く、この戦略による期待順位は約2.56となる。

n=20の場合は次のようになる。


jthresholdR(j,20)
19108.01
1876.62
1755.70
1645.05
1534.56
1434.18
1323.89
1223.64
1123.46
1013.30
913.17
813.06
713.0021
613.0017
503.0017

5人目まではパス、10人目までは1位なら、13人目までは2位なら、15番目までは3位なら(中略)採用すれば良く、期待順位は約3.00となる。

参考文献
(1)「タイミングの数理 最適停止問題」穴太克則 著(朝倉書店)

See more ...

Posted at 22:18 in 数学 | WriteBacks (0)
WriteBacks

Jul 05, 2009

Java3DのInterpolatorを使ってみる

アニメーションの制御をしてくれるInterpolatorは、デフォルトの設定で使うなら簡単で便利この上ないのだが、ちょっと思い通りに動作を変えようとすると、途端に調べる事と考える事が増える。お仕着せのシンプルなI/Fは融通が利かないというのは世の常である。簡単に使えるようになっているものでも、結局どういう時にどれくらい便利なのか、不都合な副作用が無いのかどうかを綿密に調べる羽目になることはよくある。

過去に何度か、簡単なアニメーションのためにInterpolatorを使おうとしたが、結局はやりたいことをやる方法がわからなくて、諦めてWakeupOnElapsedTimeを使って時間毎の処理を自前で書くことにしてきた。

しかし、賢い人達が練り上げた、OpenGLより上位のJava3DのI/Fが、不便であるはずが無い。必ず理由があってこのようなI/Fになっているはずである。私はすっかりSUNの信者である。使いこなせば便利で安全で堅牢で可般で最適で信頼性が高いに違いないことは、宇宙の真理として決まっているのである。
それにしてもInterpolatorを使ったアプレットはよくクラッシュする。
という訳で、一度きちんとInterpolatorを使う練習をしてみることにした。

テストアプリの起動ページへ
ソースコード

今回は、
・PositionInterpolator
・RotationInterpolator
・ScaleInterpolator
・ColorInterpolator
・TransparencyInterpolator
を使ってみた。

Interpolatorを使う上で最も多く問題になるのは、おそらく、アニメーションの周期や時間毎の動き方を決めるAlphaが作れるかどうかであろう。一定の速度で一方向に動かすだけなら簡単だが、加速・減速させたり、端まで行ったら折り返したり、端まで行ったら直角に進路変更させたりしようとすると、AlphaのI/Fをよく理解しないといけなくなる。

Alphaのコンストラクタには、一番多いもので10個の引数がある。

Alpha(int loopCount, int mode, long triggerTime, long phaseDelayDuration, long increasingAlphaDuration, long increasingAlphaRampDuration, long alphaAtOneDuration, long decreasingAlphaDuration, long decreasingAlphaRampDuration, long alphaAtZeroDuration)

これの5つ目以降が、移動、加速、減速に関するものである。
Alphaというのは基本的には最小値(default=0.0)から最大値(default=1.0)までの間を動く値である。Alphaの動きは、次のようになる。



時刻速度時刻(パラメータ表記)
(1)最小値0(1)
加速(+)
(2)max(+)(1) + increasingAlphaRampDuration / 2
等速(+)
(3)max(+)(4) - increasingAlphaRampDuration / 2
減速(+)
(4)最大値0(1) + increasingAlphaDuration
停止
(5)最大値0(4) + alphaAtOneDuration
加速(-)
(6)max(-)(5) + decreasingAlphaRampDuration / 2
等速(-)
(7)max(-)(8) - decreasingAlphaRampDuration / 2
減速(-)
(8)最小値0(5) + decreasingAlphaDuration
停止
(9)最小値0(8) + alphaAtZeroDuration

AlphaのmodeがINCREASING_ENABLEの時(default)は(1)-(5)が繰り返され、(INCREASING_ENABLE | DECREASING_ENABLE)の時は(1)-(9)が繰り返され、DECREASING_ENABLEの時は(5)-(9)が繰り返される。
従って、端まで行ったら折り返す、というのは、Alphaのmodeを(INCREASING_ENABLE | DECREASING_ENABLE)にすればできる。

残りの引数については、
loopCount: 繰り返しの回数(-1だと無限)
triggerTime: Alphaが内部で動き始める時刻
phaseDelayDuration: Alphaが内部で動き始めてから値が動き始めるまでの時間
のようである。

See more ...

Posted at 16:24 in Java | WriteBacks (0)
WriteBacks

Jun 28, 2009

秘書問題の考察(2) 期待最高順位による停止戦略

「秘書問題」(secretary problem)というのは、非常に奥深い問題で、数理科学の「最適停止理論」における1つの重要な分野となっており、未解決の問題も少なくないらしい。
前のエントリーに書いた、37%の法則が成り立つ問題は、秘書問題の最も単純なものである。秘書問題に倣って秘書を面接で選ぶという表現に変えると、前のエントリーの問題は、応募者の中の最高の人以外は外れという条件の問題であるが、現実の問題としては、できるだけ最高の人を選ぶことだけを考えればいいとは限らない。誰か1人を選ばないといけない場合、最後の応募者まで採用を見送ると、最後の人はどんな人でも選ばないといけない。その1つ前の、残り1人(まだ面接していない応募者が1人)の時点で、前回の37%の法則の戦略に従うと、今面接している人がこれまでの中の2位であっても、残りの1人が1位である可能性に期待して見送るのであるが、なるべくいい人を選びたい場合は、残り1人で暫定2位であれば、残りの1人が2位以内に入る確率は(応募者が十分多いとすると)極めて低いので、暫定2位の人を採用するのが当然である。
このことからも直感的にわかる通り、1位の人を選ぶ確率を最大化するのとなるべく上位の人を選ぶのとは一般に両立しない。なるべく上位の人を選ぶ、というのを、順位の期待値を最小化する、と定義すると、1位の人を選ぶ確率を最大化するのと順位の期待値を最小化するのを同時に満たす最適停止戦略は存在しないことが数学的に示されているらしい。

という訳で、今回は次のような問題を考える。
・n人の応募者から1人を選ぶ。
・1人ずつ面接し、その場で採用するかどうかを伝える。採用したら、残りの人は面接しない。採用しなければ、その人は2度と選べない。
・応募者の良し悪しは、既に面接した応募者と比べてどうかでしか判断できない。絶対評価で何点というのは無く、あの人と比べてどうか、しかわからない。また、優劣は必ず付けられ、必ず順位付けが可能である。
・結果としてなるべく順位の高い人が選べるようにする。

必ず誰か1人は選ばないといけないということと、最後の1点以外は、前のエントリーの問題と同じである。
この問題は既に解かれていて、nを大きくすると順位の期待値の最小が約3.87に収束する(期待順位がそれより小さい戦略は存在しない)という摩訶不思議な定理になっているのだが、あまりにも難しいので、今回は単純な方法で考えてみる。

xを残りの応募者の数(まだ面接していない人数)とする。
x=0(今面接している人が最後の応募者)の時、その人を選ぶしか無い。
x=1の時、今面接している人を含むn-1人を1位〜n-1位に並べ、残りの人がその間または1位の上、n-1位の下のどこかに入ると考えると、入る場所はnヶ所なので、これを0.5位、1.5位、…、(n-0.5)位とすると、残りの人の順位の期待値はn/2なので、今面接している人の順位がn/2より小さければ選ぶ。(n/2以下なら、としても良い)
x=2の時、同様にこれまでの1位〜n-2位に対して残りの2人が0.5位、1.5位、…、(n-2+0.5)位の(n-1)ヶ所のどこかに入ると考えると、残りの2人のいい方の順位の期待値は、1人目の順位をkとすると
(2人目がkよりいい確率)*(2人目の順位の期待値)+(1 - 2人目がkよりいい確率)*k
であり、1人目の順位がkである確率は1/(n-1)なので、
¥frac{1}{n-1}¥sum_{k=1}^{n-1}¥left¥{¥frac{k-1}{n-1}¥,¥frac{k}{2}+¥left(1-¥frac{k-1}{n-1}¥right)k¥right¥}-¥frac{1}{2}
¥frac{n}{3}-¥frac{1}{2}+¥frac{1}{6-¥frac{6}{n}}
と計算できるので、今面接している人の順位がそれより小さければ選ぶ。

x>2の時、残りのx人の最高順位の期待値を同様に求めるのは難しいので、代わりに一様分布に従うx個の標本の最小値を求めてみる。
標本が0-1の一様分布に従う時、n個の標本の最大値がyより小さい確率は、全ての標本がyより小さい確率だから、ynである。これは最大値がyになる確率の累積密度だから、確率密度関数はそれを微分してn yn-1である。従って、最大値の期待値は、
¥int_{0}^{1}y¥,ny^{n-1}dy=n¥int_{0}^{1}y^ndy=¥frac{n}{n+1}
である。標本を全て(1-標本)にすると、最小値の期待値は1 - n/(n+1) = 1/(n+1)ということになる。
これを使うと、x>2の時、残りのx人の最高順位の期待値は、0.5+(0〜(n-x)の値を取るx個の標本の最小値の期待値)=0.5+(n-x)/(x+1)と概算できるので、今面接している人の順位がこれより小さければその人を選べば良い、と考えられる。

上の計算に沿って、n=10で残りがx人の時、今面接している人が何位以内ならその人を選べば良いかを表にしてみる。










xrank
15
23.17
32.25
41.7
51.33
61.07
70.88

x=1の時(9人目)は暫定5位以内なら採用、x=2の時(8人目)は暫定3位以内なら採用、x=3の時は暫定2位以内なら採用すれば良く、x=7の時(3人目)は暫定1位でもパスすべしということになる。

n=20についてもやってみる。

















xrank
110
26.5
34.75
43.7
53
62.5
72.13
81.83
91.6
101.41
111.25
121.12
131
140.9

x=1の時は10位以内なら採用、x=2の時は6位以内なら採用、(中略)x>13の時は1位でもパスすべしということになる。

See more ...

Posted at 22:46 in 数学 | WriteBacks (0)
WriteBacks

Jun 09, 2009

37%の法則

問:n個の箱にそれぞれ異なる額のお金が入っており、どれか1つの箱のお金をもらえるとする。最大の金額は不明である。1つ箱を開けたら、そのお金をもらって終了する(残りの箱は開けない)か、そのお金をパスするかを決めないといけない。1度パスしたお金はもらえない。最高額を得る確率を最大にするには、どのように箱を選択すれば良いか。

前項と同じく、「ミネルヴァの法則」の問題のアレンジである。番組の問題では、n=20だった。正解は、最初の37%(=7個)の箱は開けても選択せず、次に出てきた、それまでの最高額より大きい金額を選択するというものである。

前項の問題と同様、たかがテレビに出てくる確率の問題が解けないはずが無いと思って、その37%を(中略)降参した。
この度、やっと解けた。

答:最初のn/e個(eは自然対数、小数点以下切り捨て)の箱の中身をパスして、次に出てきた、それまでの最高額より大きいものを選択する。

解説する。
考えられる手順は、最初のr個の箱の中身をパスして、次に出てきた、それまでの最高額より大きいものを選択する、という手順に帰着する。これ以外を考えても無駄であることは容易にわかる(注:ここでの問題では最高額以外は全て外れである)。従って、その適切なrを選ぶのが問題である。

r個をパスして最高額が得られる確率をP(r)とする。最高額がx番目(x>r)にあり、それを選択できる確率は、
・最高額がx番目にある確率=1/n
・1〜r番目の最大より大きな額が(r+1)〜(x-1)番目には無い確率=1〜(x-1)番目の最大が1〜r番目に出てくる確率=r/(x-1)
の積である。P(r)はそれを全てのxについて足し合わせたものなので、
P(r)=¥sum_{x=r+1}^n ¥frac{1}{n}¥frac{r}{x-1}
となる。

P(r)が最大になるrを求める前に、P(r+1)-P(r)を計算して、P(r)の増減を調べてみる。

なので、上から下を引いて、
P(r+1)-P(r)=¥frac{1}{n}(-1+¥frac{1}{r+1}+¥frac{1}{r+2}+¥cdots+¥frac{1}{n-1})...(1)
である。これはよく見るとrが大きくなると項が減る一方で、どこかで0より小さくなる。従って、P(r)が最大になるのは、P(r)が減り始める直前のr、すなわちP(r+1)-P(r)<0となるrである。(1)の右辺の括弧内に注目すると、
¥frac{1}{r+1}+¥frac{1}{r+2}+¥cdots+¥frac{1}{n-1}<1...(2)
となる最小のrとも言える。

さらに、nが十分大きいとして、(2)の左辺を計算し易い積分に近似する。区分求積法より、
¥sum_{k=1}^{n}¥frac{1}{n}f(¥frac{k}{n})¥simeq¥int_0^1f(x)dx
なので、f(x)=1/xとして積分範囲を調節すると、
¥sum_{k=r+1}^{n-1}¥frac{1}{k}¥simeq¥int_¥frac{r}{n}^1¥frac{1}{x}dx
と近似できる。右辺<1を計算すると、-log(r/n)<1であり、r/n>1/eすなわちr>n/eが求まる。n/eは整数ではなく、r>n/eの整数だとP(r)がちょっと減るので、n/eを超えない整数がP(r)を最大にするrである。(nが小さい時は(2)で検算した方が良いかも知れない)

37%とは、1/e(≒0.368)のことであった。
なお、この問題は(classical) secretary problemと呼ばれ、問題の解は37%の法則とか37%ルールとかと呼ばれるらしい。

See more ...

Posted at 19:11 in 数学 | WriteBacks (0)
WriteBacks

Jun 06, 2009

Coupon collector's problem

問1:おまけとしてフィギュアか何かが入っているお菓子があるとする。おまけは全m種類あり、お菓子の購入時にはどれが入っているかわからないとし、m種類は同じ確率で入っているとする。このお菓子をn回買うと、おまけが全種類揃う確率はいくらか?

これはTV東京の「ミネルヴァの法則」という昔の番組で出てきた問題のアレンジである。番組での問題は、おまけは10種類だと、100回買うと99%以上の確率で全種類揃うか(○/×)、という感じだったと思う。これの正解発表の時に、全種類揃う確率が99%になるのは66回であると、計算方法抜きで出てきた。

大学を卒業し、数式の類から遠ざかって10年以上、数学力はすっかり鈍っていたとはいえ、大学は理系卒で、高校時代は確率が得意中の得意だった身である。たかがテレビに出てくる確率の問題が解けないはずが無いと思って、その66を計算で求めようとしたが、やり方がわからず、その番組を見た日からのべ10時間以上考えて降参した。

三十半ばで再び数学に目覚めて1年、紙の上で数式をぐりぐりする右手を取り戻して、やっと解けた。

答1:
P¥left( m,n¥right) =¥sum_{k=1}^{m}{¥left( -1¥right) }^{m-k}¥,¥pmatrix{m¥cr k}¥,{¥left( ¥frac{k}{m}¥right) }^{n}
または
P¥left( m,n¥right) =¥sum_{k=0}^{m-1}{¥left( -1¥right) }^{k}¥,¥pmatrix{m¥cr k}¥,{¥left( ¥frac{m-k}{m}¥right) }^{n}
(これらはΣを展開すると同じ式、(m k)が縦に並んでるのは組み合わせ(mCkと書くのと同じ))

解説を試みる。

n個買うと、結果のパターンはmn通りである。この内、m種類揃っているのは、(m-1)種類以下しか無い場合を引いたものである。
(m-1)種類しか無いパターンは、(無い1種類のパターンmC1)×((m-1)種をn回選ぶパターン(m-1)n)、(m-2)種類しか無いパターンは、(無い2種類のパターンmC2)×((m-2)種をn回選ぶパターン(m-2)n)、…のように計算できると簡単なのだが、そうはならない。例えば(m-1)nパターンの中にも(m-2)種類以下しか無いパターンが含まれてしまうからだ。

この問題には「包除原理」(Inclusion-exclusion principle)というのが使えて、正しくは、(m-1)種類以下しか無いパターンの数は
¥pmatrix{m¥cr 1}(m-1)^n - ¥pmatrix{m¥cr 2}(m-2)^n + ¥pmatrix{m¥cr 3}(m-3)^n - ¥ldots - (-1)^{m-1}¥pmatrix{m¥cr m-1}1^n
であり、従ってm種類あるパターンの数はm^n(=mC0 m^n)からこれを引いたもの、
¥pmatrix{m¥cr 0}m^n - ¥pmatrix{m¥cr 1}(m-1)^n + ¥pmatrix{m¥cr 2}(m-2)^n - ¥ldots + (-1)^{m-1}¥pmatrix{m¥cr m-1}1^n
という+と−が交互に出てくる和になる。これを全てのパターンの数m^nで割ってΣ記号で整理したのが上記の答である。

この問題に「包除原理」が使えることを理解するには、勘が冴えてれば何も要らないかも知れないが、勘が鈍い日は次のように考えてはどうだろうか。
包除原理とは、例えばn人の人がいて、A,B,Cの3つの特徴を調べ、1人の人が複数の特徴を持つことができる時、いずれの特徴も持たない人の数が、
全体の人数−(Aの人数+Bの人数+Cの人数)+(AかつBの人数+BかつCの人数+CかつAの人数)−(AかつBかつCの人数)
となるようなことを言うものである。なぜそのようになるかというと、

この図の白い部分が求める人数だが、全体から1つの特徴を満たす人数を引くと2つの特徴を持つ人を2回引くことになり、引き過ぎた分を足すために2つの特徴を満たす人数を足すと3つの特徴を持つ人を3回足すことになり(m個の特徴を持つ人は(m-1)個の特徴を持つ人としてm回数えられるため)、以下同様に全ての特徴を満たす人数まで相殺する必要があるからである。
話を戻して、m種類のおまけが揃うパターンの数を求めるには、おまけ1が欠けてる場合をA、おまけ2が欠けてる場合をB、おまけ3が欠けてる場合をC、…と考えると、AとBが同時に起こることもあるから、A,B,…は上の図のように重なるので、どれも欠けていない場合、すなわち上の図の白い部分の数を数えようとすると、必然的に包除原理のような計算になるのである。

さて、現実的には、何個買った時にコンプリートしてる確率がいくらになるかより、大体何個買ったらコンプリートできるか、すなわち、コンプリートするまでの回数の期待値の方が興味があるのではなかろうか。

問2:おまけが全種類揃うまでのnの期待値は?

この問題は、"Coupon collector's problem"と呼ばれるらしい。日本語訳は不明だが、「食玩問題」とか「クーポン収集問題」とか呼ばれることがあるようだ。
確率分布がP(X)の確率変数Xの期待値E(X)を求める公式はΣ(X P(X))であり、上のP(m,n)より、n回目にm個揃う確率がP(m,n)-P(m,n-1)なので、nの期待値E(N,m)(以降、E(m)と省略する)はΣ(n (P(m,n)-P(m,n-1)))である。コンピューターでそれなりに大きなnまで計算すると正しい値に近づくので、一応正しい式のようであるが、これを解こうとするとひどい目に遭う。
しかし、違う発想で解くと、かなり簡単に求められる。

答2:
E(m)=m¥left(¥frac{1}{1}+¥frac{1}{2}+¥frac{1}{3}+¥ldots+¥frac{1}{m}¥right)

k-1種類まで揃ってから、k種類目が得られるまでに、平均何回かかるかと考える。
まだ得られてないのは(m-k+1)種類なので、1回買ってまだ得られてないのが得られる確率pは、p=(m-k+1)/mである。当たりの確率がpの時、当たりを引くまでの回数は、幾何分布に従う。すなわち、x回目で当たる確率はP(x)=(1-p)^(x-1) pである。このxの平均はΣxP(x)であり、等比級数の和を求める要領で解くと、ΣxP(x)=1/pとなる。従って、k種類まで揃ってからk+1種類目が得られるまでの回数の期待値はm/(m-k+1)である。
E(m)はこれを1種類目からm種類目まで足したものだから、
E(m)=m/m + m/(m-1) + m/(m-2) + ... + m/(m-m+1) = m(1/1 + 1/2 + ... + 1/m)
となる。

See more ...

Posted at 22:47 in 数学 | WriteBacks (1)
WriteBacks

クーポン収集問題

昨日のキャットストーン問題ですが、実は久々にプログラム書いてみたいというのが先にあってまともに計算していませんでした。で、1日たって計算方法がわかった、と...[link]

Posted by at 08/04/2011 09:34:38 PM

May 16, 2009

統計学復習メモ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が妥当だと思う。

See more ...

Posted at 16:00 in 数学 | WriteBacks (0)
WriteBacks

May 14, 2009

統計学復習メモ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となる。

See more ...

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

May 10, 2009

統計学復習メモ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
であり、標本分散がθの有効推定量である。
標本分散は不偏ではないので、正規分布の分散の有効推定量は不偏推定量ではないことがわかる。

See more ...

Posted at 15:29 in 数学 | WriteBacks (0)
WriteBacks

May 06, 2009

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面体の頂点の座標をベクトルの成分とすることで得られる。

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

See more ...

Posted at 16:20 in Java | WriteBacks (0)
WriteBacks

Apr 29, 2009

統計学復習メモ12: 不偏推定量とは

以前にt分布を使うパラメーターの区間推定の方法を覚えたが、この前ちょっと推定の実践のネタを思いついてやってみようとしたら、そのパラメーターが正規分布に従うものではないことに気付き、1手目でつまずいた。これでは勉強した意味が無いと思い、推定というものの基礎を勉強し直すことにした。

区間推定をするにもまず、その区間の基準となる値(推定量)を点推定しないと始まらない。その点推定量が満たすのが好ましい性質として、不偏性と一致性がある。推定量θ^がパラメーターθの不偏推定量であるとは、
E(¥hat{¥theta})=¥theta ...(1)
が成り立つことをいう。

出たE(¥hat{¥theta})=¥theta。学生時代にこのθの上の屋根の意味が理解できなかったことが、私が推定・検定の1手目でつまずいた原因の1つであることは間違いない。そもそもなぜ統計学では推定するパラメーターをθと書くのか。高校の三角関数の授業で習って以来、私にとってはθといえば角度に違いないのであり、パラメーターをθと表すのは、円周率をxと表すとか、y=ax+bのxとyが定数でaとbがグラフの縦軸・横軸であるくらいの混乱の元である。

そもそも、θ^がθの推定値だとしたら、(1)が成り立たないとはどういう状態なのだろうか?それに、θが未知だからθを推定するのに、θ^の平均がθであることをどうやって確認するのか?
…という疑問を持ったのは筆者だけなのだろうか。

これは、θが確率変数のパラメーターであり、その推定量θ^を標本から求める方法を決めた時、仮にθによって決まる母集団からの標本を使った場合、θ^の平均がθに一致するなら、そのθ^は不偏推定量である、という意味らしい。

例題として、X1〜Xnが平均μ、分散σ2の母集団からの標本として、まず標本平均X~=E(X)=1/n ΣXiがμの不偏推定量かどうかを考える。
E(¥bar{X})=E¥left(¥frac{1}{n}¥sum_{i=1}^{n}X_i¥right)=¥frac{1}{n}¥sum_{i=1}^{n}E(X_i)=¥frac{1}{n}n¥mu=¥mu
従ってX~はμの不偏推定量である。

次に、標本分散S^2=¥frac{1}{n}¥sum_{i=1}^{n}(X_i-¥bar{X})^2、不偏分散U^2=¥frac{1}{n}¥sum_{i=1}^{n-1}(X_i-¥bar{X})^2がσ2の不偏推定量かどうかを考える。

ここで、Var(X)=E(X^2)-¥left(E(X)¥right)^2の公式とVar(¥bar{X})=¥frac{1}{n}Var(X)の公式を用いると、
E(X_i^2)=Var(X_i)+¥left(E(X_i)¥right)^2=¥sigma^2+¥mu^2
E(¥bar{X})=Var(¥bar{X})+¥left(E(¥bar{X})¥right)^2=¥frac{¥sigma^2}{n}+¥mu^2
が求まるので、
E(S^2)=¥frac{n-1}{n}¥sigma^2
E(U^2)=¥frac{n}{n-1}E(S^2)=¥sigma^2
となる。従って、S2はσ2の不偏推定量ではなく、U2はσ2の不偏推定量である。

See more ...

Posted at 23:34 in 数学 | WriteBacks (0)
WriteBacks

Apr 26, 2009

続・ポリゴンできれいな球を作る

先週は面積が均一なポリゴンで球を作れなかったので、Java3DのプリミティプであるSphereクラスの球を参考にすることにした。

テストアプリ(Java3D)の起動ページへ
(今回のはカーソルキーの左右とPageUp/PageDownキーで回転します)
ソースコード
画面サンプル

Sphereクラスのコンストラクタのdivisionsの値を変えることにより、ポリゴンの密度を変えることができる。右端から反時計回りに、divisionsの値が4,8,12,16,24,32の球である。
それぞれの三角形の面積が大体等しい、きれいな球面である。どういう仕組みなのだろうか。

まず、試しにdivisionsの値を変えながら重ねて表示してみると、5,6,7,8は同じ球、9,10,11,12も同じ球であることがわかった。どうやらdivisionsの値は4の倍数に切り上げられて使われるらしい。
また、divisions=4だと正8面体で、divisionsがいくらでも、その正8面体の辺を球面に投影した(球と同じ半径の)3つの円周で球面が8分割されていることがわかる。さらに、その3つの円周の1つと平行な円周で球面が輪切りにされていることがわかる。

より細かい分割方法の特徴を明確にするために、divisions=28の球面でさらに調べてみる。(ここまでで既に出てきた3や8の倍数は避け、divisions=20は正20面体と関係するかも知れないので避ける)

赤い線で囲まれた部分が、球面の1/8である。赤い線は内接正8面体の頂点を結ぶ円周の弦、水色の線は下の方の赤い線と接する円周(以下、赤道)と平行な円周(以下、緯線)の弦である。
赤道と残り2つの経線上に、それぞれ赤い線分が7つある。おそらく、divisionsの値は赤道や残り2つの経線を何等分するかという値であろう。この1/8球面において、赤道は7分割、その上の緯線は6分割、その上は5分割、以下同様に分割されている。1/8球面を三角形に見立てると、三角形を各辺に平行な線で等分割する要領である。

という訳で、その三角形を等分割する点(線形補間)を球面に投影した点でポリゴンを作ってみる。
テストアプリの起動ページへ
ソースコード
画面サンプル

大きな三角形の頂点の辺りは密で、辺の辺りは疎になった。先週のエントリーで書いたのと同じ問題である。やはり線形補間だとうまくいかない。

そこで、Sphereクラスに習って緯線に沿った円周補間に変える。
…と思ったが、クォータニオンを使った球面補間なら上のソースコードのVector3fをQuat4fに変えるだけだし、球面補間の方が1/8球面を均等に分けられるはずなので、球面補間にする。(Sphereクラスも球面補間か?)
テストアプリの起動ページへ
ソースコード
画面サンプル

うむ、丸い。

See more ...

Posted at 18:58 in Java | WriteBacks (0)
WriteBacks

Apr 20, 2009

ポリゴンできれいな球を作るのは難しい

ポリゴンで球を作ろうとすると、真っ先に考えるのは、地球儀の緯度・経度のように2つの角度で球面上の点を表す極座標モデルであろう。しかし、それだと、赤道の辺りのメッシュが荒く、極点の辺りが細かくなり、球面が均一にならない。何とか、なるべく均等な面積のポリゴン片で球を作れないものか…と考えてみた所、次のような方法を思いついた。

球面上の3点からなる三角形の各辺の中点を取り、それらの間を結んで、三角形を4等分する。そのまま、前述の中点を球面に張り付かせるように動かす。

上から、元の三角形、4等分して球面に張り付かせた4つの三角形、そのポリゴンに色をつけたものである。合同な三角形から成る正八面体から始めて、これを繰り返していけば、均一な球面ができるのではないか?と期待して、やってみた。

テストアプリ(Java3D)の起動ページへ
ソースコード
画面サンプル

右から反時計回りに、8, 32, 128, 512, 2048, 8192個の三角形ポリゴンからなる球であり、手前が線で表示したもの、奥が面で表示したものである。

うまくできているようにも見えるが、細かくしても正八面体の面影が残っている。よく見ると、元々頂点だった所が密で、面の中心だった所が疎であり、くっきりと境目がある。これは、上記の方法で三角形を4等分して各頂点を球面に張り付ける時に、真ん中の三角形の面積が他の3つより大きくなるのが原因である。このことは、三角形の辺の中点が球面から遠い、赤道に近い位置の大きな三角形で考えてみるとわかりやすい。
正八面体に代えて正四面体で始めてみると、その影響が顕著になる。

テストアプリ(Java3D)の起動ページへ
ソースコード
画面サンプル

こうなると、却って緯線・経線の四角形ポリゴンで作る球の方が自然に見えてしまう。

テストアプリ(Java3D)の起動ページへ
ソースコード
画面サンプル

何かいい方法がありそうな気がするのだが…

See more ...

Posted at 23:35 in Java | WriteBacks (0)
WriteBacks

Apr 14, 2009

JScrollPaneを使ったミニゲーム

どうしてもJScrollPaneを使ってスクロールする何かを作ってみたくなった。

Javaアプレットの起動ページへ
ソースコード

See more ...

Posted at 21:55 in Java | WriteBacks (0)
WriteBacks

Apr 05, 2009

mod_proxy+mod_perl+MySQLでProxy Error

2月末にDebianのaptitudeが言うまま全パッケージを更新して、Debian 5.xでサポート対象外となったApache 1.3+mod_perlを手作業でインストールした後、mod_perlで動かしているMovableTypeにmod_proxy経由でアクセスすると、時々

502 Proxy Error
Reason: Document contains no data

というエラーが発生するようになった。

しばらく放置した後に起こることから、mod_perl+MySQLでよく問題になる、時間が経つとMySQLへのコネクションが無効になってしまうために起こる、いわゆる"The morning bug"絡みだろうということはすぐに想像がついて、実際MySQLを再起動すると必ず起こったので、それだと確信した。
しかし、通常はmod_perlがDBD::mysqlを介してMySQLに接続すると、"auto_reconnect"の設定が自動的に有効になるので、この問題は起こらないはずである。
(参考:man DBD::mysqlの"mysql_auto_reconnect"の項)

if either the GATEWAY_INTERFACE or MOD_PERL envionment variable is set, DBD::mysql will turn mysql_auto_reconnect on

実際、DBD/mysql.pmでログ出力してみると、mysql_auto_reconnectの値は1になっていた。しかし、mysql.pmの先は.so形式のバイナリなので、実際に再接続が試されているかは調べられなかった。

ここで、"Proxy Error"は、たまたまこの環境で結果としてそうなるだけで、同じ問題に遭遇した人が皆それに行き当たる訳ではないことは落ち着いて考えると明白だったのだが、"mysql proxy error"で検索して情報を探したため、長期に渡って悩むことになってしまった。
それでも、mod_perl+MovableTypeの組み合わせにApache::DBIを加えてコネクション・プーリングすると何かが解決した、という情報を得て、これだ!と思ってApache::DBIをインストールして組み込んでみたが、解決しなかった。

Apache::DBIのドキュメントを読むと、DBDによってはtimeout時にpingメソッドがエラーを返すので、そういう時は自前のpingメソッドを書くと問題を回避できる、と書いてあった。

Most DBI drivers have a working ping() method, but if the driver you're using doesn't have one and the database handle is no longer valid, you will get an error when accessing the database. As a work-around you can try to add your own ping() method using any database command which is cheap and safe

これだ!と思ってDBD/mysql.pmにpingメソッドを加えてみたが、そもそもpingは呼ばれていなかった。

落ち着いてApacheのエラーログをよく見ると、

[notice] child pid 7987 exit signal Segmentation fault (11)

というのが出ていた。すぐさま、そういう時はPHPをコンパイルし直せという情報を見つけて、これだ!と思って作り直すべきPHPを探すと、どこにもインストールされていなかった。

その後の調べで、DBD::mysqlのv4.007で生じたバグであり、v4.006に戻すかv4.009以降に更新すると直るらしいことがわかった。

仕方なく、CPANを使ってDBD::mysqlをコンパイルしようとすると、mysql_configが無い、というエラーが出た。Debian 5.xではmysql_configはlibmysqlclient15-devというパッケージをインストールしないと使えないらしい。これはわかりづらい。

結局、CPANでDBD::mysqlのv4.010をインストールすることによって、うちの"Proxy Error"は解消した。

See more ...

WriteBacks

Hi, Excellent post

Pretty insightful post. Never thought that it was this simple after all. I had spent a good deal of my time looking for someone to explain this subject clearly and you’re the only one that ever did that. Kudos to you! Keep it up

Posted by Faye Ingram at 2010/7/20 19:25:46

Mar 29, 2009

簡易MMLプレーヤー

javax.sound.midiを使うと割と簡単に音階を鳴らせることがわかったので、簡単なMML(Music Macro Language)プレーヤーを作ってみた。BASICのPLAY文でお馴染みの、ドレミファソラシドがCDEFGABCで表されるあれだ。

Simple MML Playerの起動ページへ
ソースコード

とりあえず、CDEFGとでも入れてPlayを押すと、何か鳴ると思う。

以下の文字をサポートしている。(大文字と小文字の区別は無し、[]は省略可、xは数字を表す)
・C,D,E,F,G,A,B: ドレミファソラシド
 形式:(C,D,E,F,G,A,Bのどれか)[+/-][x][.][&]
 +/-は#/♭、xは長さ(全音符の何分の1か)、"."は符点(長さ1.5倍指定)、&はタイ/スラー
 例)C4 --- 四分音符のド、D+2. --- 符点二分音符の#レ
  B2&B8 --- 二分音符+八分音符の長さのシ
・R: 休符
 形式:R[x][.]
 xは長さ(全休符の何分の1か)、"."は符点
・O[x]: オクターブ(デフォルトは4)
・L[x]: 長さ省略時の音符/休符の長さ(1-64、デフォルトは4)
・V[x]: 音の大きさ(0-100、デフォルトはチャンネル0が100、チャンネル1が75)
・T[x]: テンポ、1分当たりの4分音符の数(デフォルトは120)
・@[x]: 音色(0-127、デフォルトは0(グランドピアノ))(参考資料
・<,>: オクターブを1つ下がる/上がる
・(): 和音
 形式:([C/D/E/F/G/A/B/+/-/>/<]の羅列)[x][.][&]
 例)(CEG)4 --- 四分音符の長さでドミソを同時に鳴らす
  (A>C+E) --- ラ、1つ上の#ド、ミを同時に鳴らす
・","(カンマ): チャンネル区切り(チャンネルは2つのみ)
 例)C8D8E8D8C4,E4F4G4 --- チャンネル0でドレミレド、チャンネル1でミファソ

もちろん独自仕様である。いわゆる作者 is Godである。

1行の中に","で区切って2つのMMLを書くと、それぞれが別のチャンネルで同時に鳴る。
同じ行の2つのMMLは必ず同時に始まる。一方が短い場合、長い方が終わるまで次の行が演奏されない。

エラー処理は極めて曖昧である。というか、ほとんど何もしていない。不正な文字列を入力や30和音などの無茶な指定をした場合にきちんと鳴るかどうかは無保証である。

ブラウザ上でなく、ローカルのファイルシステムから実行すると、演奏する前にMIDIデータ(SMF形式)がsimple_mml.midというファイルに保存される。

スライドバーの音量設定は、次の演奏から有効になる。

以下は、javax.sound.midiに関する覚え書きである。
・概略
SequenceにTrackを作って、Track#add()で時系列のMidiEventを登録して、Sequencerに渡してSequencer#start()を呼ぶとMIDIシーケンスが演奏される。
MidiEventは、ShortMessageのNOTE_ONがある音階の鳴り始め、NOTE_OFFが鳴り終わりである。ピアノの1つの鍵を押すのと離すのに相当する。音階は0〜127の整数で、60が基準のド、61がド#、62がレ…である。

・テンポについて
JDK 1.6のドキュメント(API spec, jdk-6-docs/technotes/guides/sound/)にもMIDI 1.0の仕様書(MIDI 1.0 Spec, SMF spec)にも見つけられなかったので正しいやり方かどうかわからないが、筆者のSun JRE1.6の環境では、javax.sound.midiでもmeta eventのFF 51で四分音符の長さを設定できるようである(デファクトスタンダードの類?)。MetaMessageで登録する場合は、type=0x51とし、 dataはbig endianの3バイトで、四分音符の長さをμ秒で指定する。
このメタイベントを使う場合は、Sequenceを作る時のdivisionTypeをPPQとする。
このメタイベントを使う以外の方法としては、divisionTypeをSMPTE_30、resolutionを100とかにしてtickの長さを絶対時間にして(この引数だと1 tickが1/3000秒になる)テンポに応じて音の長さをtickの整数倍で指定していく方法しか思いつかなかった。この方法だと誤差が大きいし、tickの値がすぐ大きくなるのでMIDIシーケンスのデータが長くなってしまうので、好ましくない。

・演奏停止について
Sequencer#stop()を呼んでも、既に鳴り始めた音は鳴り止まないことがある。確実に音を止めるために、Sequencer.close()を呼ぶようにした。
また、シーケンスが終了しても最後の音が鳴り続けることがあるので、シーケンスの最後にメタイベントのFF 2F(トラック終端)を付けて、それをMetaEventListenerで受けて、音を止めるようにした。

See more ...

Posted at 18:35 in Java | WriteBacks (0)
WriteBacks

Mar 20, 2009

統計学復習メモ11: なぜ特異値分解で主成分が求まるのか

前々項に書いたように、主成分分析は特異値分解を使っても計算できる。なぜだろうか。
それについても答えを得たので、ついでにメモしておく。

n個のサンプルデータを横に並べたものを

(m<n)とし、Pを単位主成分を横に並べたもの

とし、YをXの主成分とすると、
Y=P^TX
である。

特異値分解定理により、Xは何らかの正規直交行列U,Vと、m行m列部分が対角行列でそれ以外が0である行列Σによって、
X=U¥Sigma V^T
と分解される。それを応用すると、
XX^T=(U¥Sigma V^T)(U¥Sigma V^T)^T=U¥Sigma V^TV¥Sigma^TU^T=U¥Sigma¥Sigma^TU^T...(1)
が成り立つ。XXTはm×mの正方行列なので、ΣΣTもm×mの正方行列であり、
X=U¥Sigma V^T...(2)
と置くことができる。(σ1≧σ2≧...≧σmである)

前項に書いたように、Yの共分散行列は
Cov(Y)=P^TCov(X)P
であり、従って
Cov(X)=P\;Cov(Y)P^T
であり、Cov(Y)はCov(X)の固有値を対角成分に持つ行列である。Cov(Y)=Λを

1≧λ2≧...≧λm)と置くと、
...(3)
となる。一方、
Cov(X)=¥frac{1}{n}XX^T
なので、(1)と(2)を使うと、
...(4)
という形にすることができ、(3)と(4)を見比べると、PもUも正規直交行列なので、U=P、
すなわち ¥frac{¥sigma_k^2}{n}=¥lambda_k
の関係があることがわかる。さらに、
Y=P^TX=U^TU¥Sigma V^T=¥Sigma V^T
なので、Xの主成分YはXを特異値分解したものの一部として求められることがわかる。

See more ...

Posted at 19:42 in 数学 | WriteBacks (0)
WriteBacks

Mar 14, 2009

統計学復習メモ10: なぜ共分散行列の固有ベクトルが単位主成分なのか

前項に書いた通り、主成分分析における主成分の単位ベクトルは、共分散行列の固有ベクトルとして求まる。そのこと自体に昔から興味があったので、主成分分析の復習ついでに考察してみる。

まず、最小2乗法で考えてみる。簡単のために2次元で考える。n個のサンプルデータを

とし、第1主成分の単位ベクトルを

とすると、Xに対応する主成分軸上の第1主成分Yは
Y=P^TX
であり、そのYを元の座標系に戻したものX~は
¥tilde{X}=PY=PP^TX
である。このことは、高校で習った一次変換を思い出してやってみるとわかる。このX~が、Xを第1主成分の軸上に射影したものであり、これとXとの距離が、最小にしたい誤差ということになる。その誤差Eを、Xを直交座標とした場合の距離の2乗とすると、
E=(¥tilde{X}-X)^T(¥tilde{X}-X)
であり、p12+p22=1に注意すると、これは
=¥sum_{i=1}^n((p_1(p_1x_{1i}+p_2x_{2i})-x_{1i})^2+(p_2(p_1x_{1i}+p_2x_{2i})-x_{2i})^2)
=¥sum_{i=1}^n(p_2x_{1i}+p_1x_{2i})^2
と計算できる。

このp1とp2の2次関数であるEの最小値ををp12+p22=1の条件下で求めるので、「ラグランジュの未定乗数法」というやつを使う。すると、G=p12+p12-1と置くと、

である。これを整理すると、

となる。これの最初の行列は、よく見ると、Xの共分散行列

の逆行列のスカラー倍であるので、上の方程式は、両辺に左側からCov(X)をかけると

(aはスカラー)となり、a/λをλと置き直すと、

が得られる。従って、固有ベクトルの定義より、PはCov(X)の固有ベクトルであることがわかる。

さて、Xが2次元の場合については辛うじて求まったが、Xが3次元になると、この方法ではEが複雑になり、非常に辛い。筆者はMaximaに式の変形をさせながら3次元についても同じであることを調べようとしたが、E-λGを偏微分した所で変数が目測で200個以上ある状態になり、なんとなく雰囲気はあったが挫折した。この方法ではXをm次元に一般化するのはかなり困難か、何か特殊な技術が必要だと思う。

そこで、次に分散最大で考えてみる。こちらは幾分容易である。
n個のサンプルデータを

とし、第1主成分の単位ベクトルを

として、第1主成分

の分散
V=¥frac{1}{n}YY^T=¥frac{1}{n}(P^TX)(P^TX)^T=¥frac{1}{n}P^TXX^TP
を最大にするPを求めることを考える。Xのi行目とj行目の共分散をcov(Xi,Xj)と書くと、
V=¥sum_{i=1}^{m}¥sum_{j=1}^{m}p_ip_jcov(X_i,X_j)
とも書ける。このp1...pmの2次関数であるVの最大値をp12+p22+...+pm2=1の条件下で求めるので、やはりラグランジュの未定乗数法を使う。すると、G=p12+p22...+pm2-1と置くと、

である。これを整理すると、

すなわち

が得られ、PがCov(X)の固有ベクトルであることがわかる。

第2主成分以降についても、それより前の主成分と直交するという条件を加えて、同様の方法で解くと、やはりCov(X)の固有ベクトルであることがわかり、第1主成分から順番に分散最大の固有値が取られていくので、第k主成分はCov(X)のk番目に大きい固有値に対応する固有ベクトルだと求まる。

See more ...

Posted at 15:17 in 数学 | WriteBacks (0)
WriteBacks

Mar 04, 2009

Apache1.3+mod_perlのコンパイル

偶然にも、mod_perlという名前を久々にここに書いた日の夜に、このサーバーのmod_perlが動かなくなった。いつも通りに、Debianのaptitudeを使って"Upgradable"なパッケージを全てアップグレードしたら、mod_perlが消えてしまったのである。
先月半ばにDebianの5.0がリリースされ、Apache 1.3がサポート対象外になったのと関係しているものと思われるが、Apache 1.3自体は"Obsoleted package"に分類されながらも残っていたのに、なぜmod_perlだけ消えてしまったのだろう。Perlがバージョンアップされたからだろうか。アップグレードに伴い、依存関係の解消のために消されるパッケージとして表示されていたのだろうが、今回は"Upgradable"が多量だったのでろくに見なかったのが失敗だった。

仕方なく、面倒だがこの際Apache2+mod_perl2に移行するかと思って、mod_perl2をインストールして動かしてみると、mod_perlで動かしたいMovableTypeがmod_perl 1.xにしか対応していないことがわかった。

何としてもmod_perl 1.xをインストールしようと思って、Debianのapache1.3+mod_perlのパッケージを拾ってインストールしようとすると、それがおびただしい数のパッケージに依存していることがわかった。ついでに、依存関係にあるパッケージを1つ1つインストールすると、次々にそれらが"Obsoleted packages"に入れられてしまい、完全にやる気をそがれた。

という訳で、Apache 1.3とmod_perlのソースコードをコンパイルするしか無くなった。
このサーバーではDebianのaptitude頼りでほとんどオープンソースのコンパイルというものをしなかったので、簡単にコンパイルできるとは思っていなかったが、予想に反してあまりにも簡単にコンパイルできた。
感動したので、その手順を記録する。

1. 適当な所から以下の2つのソースを取得、展開
apache_1.3.41.tar.gz
mod_perl-1.0-current.tar.gz

2. mod_perl-1.XXディレクトリに移動して、以下を実行

perl Makefile.PL APACHE_SRC=../apache_1.3.41/src DO_HTTPD=1 USE_APACI=1 EVERYTHING=1 ADD_MODULE=proxy
make
make install
cd ../apache_1.3.41/
make install

3. Apache::Request, Apache::Cookieをインストール

perl -MCPAN -e shell
> install Apache::Request

以上である。
Apache::Cookieは、Apache::Requestをインストールすると自動的にインストールされる。
mod_perlをmakeすると、自動的にapacheもコンパイルされる。

Makefile.PLのオプションは、mod_perlのページの例を参考にした。例の通りではmod_proxyが有効にならなかったので、ADD_MODULE=proxyを追加した。

See more ...

WriteBacks

Mar 01, 2009

Java servletをPerl CGIに移植してみる

Javaのservletは、HTTPリクエストを受けてテキストを返す場合(HttpServlet)に限ると、役割はCGIと同じである。実行プロセスが常駐している以外は、CGIと大差が無い。CGIでも、例えばApache+mod_perlを使えばPerlのモジュールをapacheのプロセスに常駐させることができ、さらにJava servletと差が無くなる。従って、Java servletにすべきかCGIか、またCGIならどの言語を選ぶべきかは、動作速度と作り易さにかかっているのではないか…と思っていた。

先日、3択首都当てクイズのservletのコードを見直していて、文字列処理等でJava特有の面倒臭さがあり、PerlのCGIならかなり楽に書けるんじゃないか、と思った。
筆者はPerlは結構昔から使っているのだが、PerlのCGIを真面目に作ったことが無かったので、CGI.pmの勉強を兼ねて、3択首都当てクイズのサーバーサイドプログラムのCGI版を作ってみた。

CGIのソースコード
対応するクライアント(Java applet)側コード
対CGI版首都当て3択クイズの起動ページ

確かに幾分楽に書けたが、常駐型ではないことが前提になっている。当然、極端に遅い。速度を確保するには常駐型にする必要があるが、作ってみて思ったが、これを常駐型にするのは厄介である。
筆者はPerlのCGIを常駐型にする方法としてApache+mod_perlしか知らないので、これを前提に書くが、PerlのCGIはJavaと違って使用メモリ量を制限するのが容易ではない。このwebサーバーにもmod_perlを導入しているが、メモリに関しては結構苦労があった。適切なタイミングでガベッジコレクションさせる術が無いため、何回か毎にリセットするしか無いのである。1つのコンテキストで処理していると、リセット中のスループットが低下するため、2つのコンテキストで処理するようにしたいが、次にデータベースとの接続の管理方法が問題になる。もしTomcatのようなコネクションプールの仕組みが使えるとしても、そこまでしてマルチスレッドにしたいか?と悩んでしまう。
筆者はApache+mod_perlは重要な選択肢の1つだと考えているが、CGIをmod_perlでインターネットに公開するのは、怖くてできない。一般論として、CGIを常駐型で動かすのは、強引で裏技的な手段であり、好んで採用するものでは無いと思う。やはり、常駐型ならJavaのサーブレットで、非常駐ならCGIで作るのが良いのだろう、と考え直した。

See more ...

Posted at 23:04 in PC一般 | WriteBacks (0)
WriteBacks

Feb 21, 2009

ラグランジュの未定乗数法

統計学を復習していて、どうやらそれを使わないと非常に難しいらしい問題に突き当たったので、ラグランジュの未定乗数法を復習する。
その響きで一瞬怯んでしまうが、前提条件にさえ気を付ければ使う分には何も難しくない、ロピタルの定理のような便利技である。

どこにでも書いてあるので、わざわざここに書くまでもないが、ラグランジュの未定乗数法とは、x1...xnの関数fが束縛条件g=0の下で極値を持つ時、その極点のx1...xnは、f-λg(λは仮置きの定数)をx1...xnで偏微分したものが0になる値として求めることができるというものである。
数式で書くと、

である。ロピタルの定理同様、考えてみるとミラクルである。

例題として、2008年2月1日にフジテレビから放送された「たけしのコマネチ大学数学科」の第75回の問題を取り上げる。

問:

このような円柱と三角錐を組み合わせた鉛筆形の立体があり、体積が決まると、円柱の半径と円柱・円錐の高さは鉛筆形の表面積が最小になるように決まるとする。円錐部分の稜線の長さ(頂点から円錐面に沿って引いた線の長さ)が3の時、円柱部分の半径はいくらか。

解答:
円柱の半径をx、円柱の高さをy、円錐の稜線の長さをzとすると、表面積Sと体積Vは、中学校で習った公式を思い出すと
S=2¥pi x^2+2¥pi xy+¥pi xz
V=¥pi x^2y+¥frac{¥pi x^2}{3}¥sqrt{z^2-x^2}
と求まる。

まずは失敗の例から。普通に考えると、例えばVの式をyについて解いて、Sの式に代入してyを消して、
S=¥frac{2¥left( 3V-¥pi {x}^{2}¥sqrt{{z}^{2}-{x}^{2}}¥right) }{3x}+¥pi xz+¥pi {x}^{2}
このSを最小にするx,zを求めれば良いのだが、よく見ると1つの式に対して変数がx,z,Vと3つある上、例えばxについて偏微分すると
¥frac{¥partial S}{¥partial x}=-¥frac{¥sqrt{{z}^{2}-{x}^{2}}¥,¥left( 6¥,V-3¥,¥pi ¥,{x}^{2}¥,z-6¥,¥pi ¥,{x}^{3}¥right) +2¥,¥pi ¥,{x}^{2}¥,{z}^{2}-4¥,¥pi ¥,{x}^{4}}{3¥,{x}^{2}¥,¥sqrt{{z}^{2}-{x}^{2}}}
(右辺はMaximaの出力)となり、いわゆる強靭な計算力が必要な状況に陥る。ちなみに、Maximaでz,Vについても偏微分させてそれぞれ=0の連立方程式を解かせると、"system too complicated."というエラーを出してgive upする。

そこでラグランジュの未定乗数法である。
Vという条件の下でSの最小値を求める問題なので、

である。(正確にはVは(V-体積)だが、偏微分で体積が消えるので不要)
yについての式を実際に偏微分すると

となり、λ=2/xが得られる。次に、zについての式を偏微分してλを消去すると、

となり、これにz=3を代入すると、x=√5と求まる。
Vもyもわからないままで良い。なんということであろうか。

See more ...

Posted at 22:17 in 数学 | WriteBacks (5)
WriteBacks

Maximaの使い方を探していて、立ち寄りました。 ラクランジュ未定乗数法は確かに簡単ですね。 ただ、今回の場合は、最初の方法では解き方を間違っているために、おかしなことになっています。 単に、xを求めるだけであれば、Vが外生変数ですので、問題を解く時には定数で、 あるVに対するx、zの関係がでてきます。Vで偏微分して、解こうとしたら解けません。 1行目から4行目までは基本的に一緒です。 SYというのは、体積をYについて解いたものを代入した式です。 ここで、体積がzに対して最小化するように決まるのですから、zで偏微分してxについて解けば x=5^(1/2)/3zが求まりますので、z=3をいれれば、x=5^(1/2)となります。 S: %pi*x^2 + 2*%pi*x*y + %pi*x*z; solve(%pi*x^2*y + %pi*x^2*sqrt(z^2-x^2)/3=V,y); SY:%pi*x*z+2*%pi*x*(3*V-%pi*x^2*sqrt(z^2-x^2))/(3*%pi*x^2)+%pi*x^2; assume(z>0); solve(diff(SY,z)=0,x);

Posted by ぼんど at 01/31/2010 05:03:07 PM

ぼんど様ご指摘ありがとうございます。ほんとですね。Vはxやzと独立ではないので、Vで偏微分するのは意味が無いし、ここではVは定数扱いですね。紙の上でyを消したSの式(SYに相当)の形を見て計算不能と決めつけてましたが、これをzで偏微分するとVが消えるので、確かにdiff(SY,z)=0として解けますね。しかし、たまたまそうであっただけで、もしSYがxやzで偏微分してVが消えない形なら、やはりSYを最小にするx,zを求めるのは困難なのでは…

Posted by ynomura at 01/31/2010 10:18:39 PM

Vは与えられる定数で、その中で解きますので、Vが残った時は、V=の式を作って、それを連立させて消すことで、xとzとの式に帰着できます。 dS/dx=(d/dx) f(x,z,V)=0、dS/dz=(d/dx) f(x,z,V)=0となり、どちらかの式からV=を作り、これを使ってVを消去して、x=g(z)という式を導けばOKです。 ちなみに、未定乗数法というのは、代入法を簡単にしたと考えることができます。 未定乗数法が真価を発揮するのは、与えられた条件が陽に解けない場合です。 f(x,y,z)=0という式があったときに、これを、y=f(x、z)といった形に変形することを陽といいます。f(x,y,z)=0は、各変数の関係が式の中に隠れていますので陰(関数)です。 今回の問題の場合は、y=の形に式が変形できていますので、未定乗数法を使わなくても解くことが容易です。 ところが、式によっては、y=という形に式を変形して、代入できない場合があります。 こういうときは、未定乗数法を使うことで、y=という式を意識しないで問題を解くことができます。(ただし、陰関数定理によって、具体的には陽にできないけれども、理論的にはy=が存在することは保証されないといけません) 例えば、z=f(x,y)という式の極値を、g(x,y)=0という条件下で考える時、条件を書き換えて、y=j(x)という形にできるなら、これを代入することでz=f(x,j(x))という関数になりますので、これを微分して dz/dx=df/dx+(df/dy)(dj/dx)=0・・・(1)   とすれば、xが求まります。 ここで、面白いのは、y=j(x)としなくても、g(x,y)=0を用いて未定乗数法を行うと、dj/dxがj(x)という関数を使わずに表現することができます。 j(x)がどんな形かわわかりませんが、あるという前提で、g(x,j(x))=0とします。これを微分すると、dg/dx+(dg/dy)(dj/dx)=0となりますので、式変形すると dj/dx=-(dg/dx)/(dg/dy)です。 この結果を、(1)に代入すると dz/dx=df/dx+(df/dy)(dj/dx)=df/dx-(df/dy){(dg/dx)/(dg/dy)}=0 となります。 この式のすごいところは、g(x,y)=0という陰関数が与えられれば、それを各変数について陽の形にすることができなくても、変形して代入した時と同じ微分の条件が作れてしまうところです。 ラグランジュ未定乗数法は、未定乗数を使うことで、機械的に上記の計算ができるようになっています。また、未定乗数λにも意味があり、例えば今回の問題でλの値を求めると、それはVを変化させた時に、どれくらいSが変化するかを表しています。

Posted by ぼんど at 02/01/2010 05:39:38 PM

わかり易いご説明ありがとうございます。微分の公式を繙きつつ数式の変形を追って、感動致しました。 なるほど、(f-λg)'=0のλってのはdf/dgのことだったんですね。だからf+λgじゃなくてf-λgなんですか。 2つの陰関数の関係式が、1つの変数だけ陽でそれ以外は陰のままで表すことができるのがすごいですね。変数を1つだけ順に表に出してくることができる感じだと思いました。 dS/dx=0とdS/dz=0からのVの消去についても、理解しました。例えばそれらの2つの微分方程式の両方がVの2次関数だったりするとV=w(x,z)の形にするのは容易ではないのではないか、という疑問だったのですが、Vの式の%pi*x^2*yの部分を%pi*x^2*y^2や%pi*x^2*sqrt(y)にしたり、yの項を増やしたりしてMaximaで試した所、いずれも容易にV=w(x,z)の形にできました。y=f(x,z)と書けた上で本当に上の連立微分方程式からVが消去困難だと、きっとf(x,z)が陽関数とは言えない形をしているのだろうと思いました。

Posted by ynomura at 02/03/2010 09:32:31 PM

ラグランジュの未定乗数法

Lagrange multiplier EMANの物理学・解析力学・ラグランジュの未定乗数法 ラグランジュの未定乗数法 (Weblog on mebius...[link]

Posted by at 05/23/2011 08:06:45 PM

Feb 15, 2009

統計学復習メモ9: 主成分分析の計算方法

主成分分析の話になると大抵、固有ベクトルか特異値分解のどちらかが出てくる。主成分の単位ベクトルが固有ベクトルとして直ちに求められ、主成分の大きさが特異値分解を使って直ちに求められるからだ。

それにしても、行列には固有値や固有ベクトルがあるというのも、統計には分散や共分散があるというのも、理系の大卒なら誰でも知ってることだろうが、「共分散行列の固有ベクトルが主成分の単位ベクトルになる」と言われると、全く関係ない3つのものが突然シンプルに組み合わせられた感じで、衝撃的である。まるで
e^{-i¥pi}=-1
のような神秘的なものを感じるのは筆者だけだろうか。

そのこと自体は、至る所に証明というか説明があり、その方法も幾通りもある。基本的には、あるデータのど真ん中を通る直線として、最小2乗法のように2乗誤差を最小にする直線を求める代わりに、その直線方向の分散が最大になる(その直線上に射影した場合の分散が最大になる)直線を求める、という考え方である。そう言われればそういう考え方もありかな?とも思ったりするし、やっぱり最小2乗法で求める方がいいんじゃないの?とも思ったりするが、結局、主成分は分散最大で求めても最小2乗法で求めても同じ結果になるらしい。

とりあえず、固有ベクトル計算でやってみる。
m次元のデータn個を

(iは標本番号、pはパラメーター番号の意味)とし、共分散行列の固有値を

とし、対応する固有ベクトルを

とすると、W1〜Wmが第1主成分〜第m主成分、λ1〜λmがそれぞれの寄与率となる。
X,Wは行列である。プログラミング用語の「待ち行列」とかのqueueではなく、線形代数のmatrixである。フフン、行列くらい知っとるわい、という顔をしつつ、厄介だなあ、と思わざるを得ないが、主成分分析するのに行列は避けて通れない。
前回使用したデータについて、Maximaに計算させてみる。

Maximaへの入力

load(descriptive);
load(lapack);
fpprintprec:6;
X: read_matrix("data2");
[L,W,dummy]: dgeev(cov(X),true,false);
L;
W;

Maximaの出力(上のをbatch()で読み込んだ)
(%i7) L
(%o7) [.330129, .0531125, .00529463]
(%i8) W
[ .714134 .697955 - .0535787 ]
[ ]
(%o8) [ 0.4882 - .551443 - 0.67644 ]
[ ]
[ .501671 - .456912 .734546 ]

固有値(寄与率)は絶対値が一番大きいのが1つ目、二番目に大きいのが2つ目なので、第一主成分の単位ベクトルは(0.714, 0.488, 0.502)、第二主成分の単位ベクトルは(0.698, -0.551, -0.457)である。3つ目の固有値は小さいので、3つ目の固有ベクトルは無視することにする。
今回はたまたま大きい順に並んでいるが、固有値は絶対値の大きさ順に並ぶとは限らないので注意が必要である。
もしdgeev()が使えない場合は、eigens_by_jacobi()を使う。
(%i9) [L,W]:eigens_by_jacobi(cov(X)); L; W;
(%o10) [.330129, .0531125, .00529463]
[ .714134 - .697955 - .0535787 ]
[ ]
(%o11) [ 0.4882 .551443 - 0.67644 ]
[ ]
[ .501671 .456912 .734546 ]

2つ目の固有ベクトルの符号が逆だが、次に求める主成分の大きさの符号が逆になるだけなので、どちらでも良い。

次に、第一主成分、第二主成分を軸に取った時にXの各標本がどういう値になるか(主成分の大きさ)を求める。主成分の大きさYは
Y=XW
なので、

W: submatrix(W,3); /* 弟三主成分を削除 */
Y: X.W;

で求まる。これによって得られるYをグラフにすると、次のようになる。

前回のグラフと見比べると、視点が反対(前回のグラフ作成の時はengens_by_jacobi()を使ったため)だが、例えば楕円外の点の分布は大体似ていることがわかる。
これによって、今回使った3次元のデータが、情報の欠落を最小限にして、2次元で分析できる。

特異値分解による方法でもやってみる。
今度はXの中の標本は横に並べる。縦方向の1列が1つの標本である。

このXが特異値分解によって
X=U¥Sigma V^T
(Tは転置)と3つの行列(U=左特異行列、Σ=特異値行列、V=右特異行列)の積になる時、Uが主成分の単位ベクトルを縦に並べたもの、Σの対角成分(特異値)が[n×寄与率]の平方根となり、標本の主成分の大きさYは
Y=U^TX=¥Sigma V^T
となる。
Maximaで特異値分解するには、LAPACKのdgesvd()を使うと一発である。(LAPACKが使えないと、かなり厳しそうである)

Maximaへの入力

load(lapack);
fpprintprec:6;
X: transpose(read_matrix("data2"));
[s,U,V]: dgesvd(X,true,true)$
U;
s;

Maximaの出力
(%i7) U
[ - .713141 .698967 - 0.053619 ]
[ ]
(%o7) [ - .488765 - 0.55059 - .676727 ]
[ ]
[ - .502532 - .456395 .734279 ]
(%i8) s
(%o8) [5.77344, 2.30564, .727852]

(Vは100x100の行列なので省略)
行列Uに、先程と大体同じ主成分の単位ベクトルがあることと、3つ目の特異値が小さく、寄与率が低いことがわかる。

Maximaへの入力(続き)

S: zeromatrix(2,100); /* 行列Σ(3x100)の上2行分 */
S[1,1]: s[1]; /* 対角成分1 */
S[2,2]: s[2]; /* 対角成分2 */
Y: S.V;

これによって得られるYをグラフにすると、次のようになる。

さらに反対向きになったが、そういうものである。

さて、固有ベクトル計算による方法と特異値分解による方法をどう使い分けるかだが、一般に、特異値分解を使う方法の方が計算誤差が少ない、と言われる。固有ベクトル計算による方法だと、共分散行列にする時に一気に情報量が落ちるので、感覚的には納得できるが、使用するモジュールやライブラリの精度による所もあるので、一概には言えないと思う(例えばすごい桁数の共分散行列から固有ベクトルを求めれば精度は上がるはず)。

また、ライブラリの特性の問題かも知れないが、eigens_by_jacobi()やLAPACKのdgeev()では固有ベクトルが固有値の絶対値の順に並ぶとは限らないので、運が悪いと固有値を見て固有ベクトルを並び替える必要が生じる。それに対し、LAPACKのdgesvd()は特異値が必ず大きい方から並ぶので、LAPACKを使うなら特異値分解による方法の方が楽かも知れない。

See more ...

Posted at 19:15 in 数学 | WriteBacks (2)
WriteBacks

大変参考になりました

Posted by 岩淵 at 12/22/2010 05:55:32 AM

ありがとうございます。

Posted by ynomura at 12/22/2010 10:51:28 PM

Feb 14, 2009

EclipseにTomcatプラグイン導入

仕事で科学技術に恩返しできないならインターネットで、と思ってサンプルプログラムをここに書いてきたが、1年前に公開した「3択首都当てクイズ」のソースコードは掲載しなかった。なぜかというと、大体動くようになったが最後、そのコードを触るのが嫌になったからである。なぜ嫌になったかというと、Javaのサーブレットの開発が大変面倒だったからだ。

この24時間運転のwebサーバーは玄箱HGという小型の機械で動いている。いわゆるパソコンではなく、消費電力は低いが計算能力も低い。サーブレットのプログラムを書き換えると、それをこの機械で動かすためにTomcatを再起動するのに2分以上かかる。加えて、思い通りに動かないと、解析するにはTomcatのログを見るしか無い。ちょっといじって動かなくなると、調べるのも嫌になるのである。

そんな苦い思い出を、パソコンにインストールしたEclipseのTomcatブラグインが吹き飛ばした。どうせServlet APIの補完が効くのがいいくらいでしょ、と期待せずに一応インストールしてみたのだが、とんでもなく便利である。何が便利かというと、servletをデバッグモードで動かせることと、servletのソースコードを書き換えると勝手にコンパイルされてTomcatに組み込まれて動作開始することである。

Javaのサーブレットのソースコード上でブレークポイントを設定してTomcatを起動しておくと、クライアントがEclipseのデバッグモードで起動したJavaアプレットであろうが普通のwebブラウザであろうが、クライアントからの要求が来ると、サーブレットの動作がブレークポイントで一時停止する。デバッグモードのTomcatはバックグラウンドで動くので、Tomcatを起動した後も何の支障もなくEclipse上で他の作業ができる。クライアントのJavaアプレットにもブレークポイントを設定しておくと、サーブレットとアプレットとを交互に一時停止させることもできる。両方の動作を時間順に追跡できるので、大変わかりやすい。

さらに、一時停止したサーブレットのコードを、その場で書き換えて、そのメソッドの先頭からやり直させることができる。既に到着したクライアントの要求はそのままに、書き忘れた処理を追加して続行させたり、試しに処理をコメントアウトして続行させたりできるのである。ある1つの要求に対して、処理を書き換えながら色々な処理を試せるのが便利だ。
サーブレットのコードを変えたつもりなのに動作が思い通りにならない時、きちんとコンパイルし直したか、正しくTomcatに組み込んで再起動したかをいちいち気にする必要も無くなった。ブレークしてその場でコードを追加してそこを通る様子を見ればいいのである。

という訳で、EclipseのTomcatプラグインを使って見るとボロボロだった3択首都当てクイズのコードを、デバッグして整形したので、その準備過程の一部を書き残す。

●Tomcatプラグインのインストール
インストール方法は色々な所に書かれているので、省略する。注意点として、一次配布元のSysdeoのサイトは長期間アクセス不能になっているので、Eclipse Wikiのベージからダウンロードするのが良い。
今回はV3.2.1をEclipse V3.4にインストールした。

●プラグインの設定方法
EclipseのPreferencesの画面で、Tomcatの設定を埋める。
といっても、今回明示的に設定したのは、以下の項目だけ。
・Tomcatバージョン
・Tomcatホーム
・JVM設定のJREのバージョン
・Tomcatアプリケーションマネージャーのユーザー名とパスワード

●Tomcatプロジェクトの作成
取っ掛かりの部分はeclipse@Wkiのページに書いてある。
JavaのソースはWEB-INF/srcに、web.xmlはWEB-INFに、JSPのソースはプロジェクトの直下に配置する。

●web.xmlの作成
Tomcatのサンプルをコピーしてテキストエディターで書き換えても良いが、今回は、Eclipseで.xmlを開くと自動的に起動する、Eclipse組み込みのXMLエディターを試してみた。

1. EclipseでXMLファイルを新規作成(wizardはXML/XMLを使用)
2. 最初から入ってる1行("?=? xml")を消して空にする
 (これが残っているとサーブレットが起動できなかった)
3. 右クリックでAdd child→New elementを選択して"web-app"エレメント追加
4. 右クリックでEdit namespacesを選択
 Edit Schema Informationの画面でAdd→Specify New Namesp...を選択
 Prefix: (空)
 Namespace Name: http://java.sun.com/xml/ns/j2ee
 Location Hint: http://java.sun.com/xml/ns/j2ee/web-app_2_4.xsd
 を指定して追加
 これで定型部分が埋まってweb.xml用の補完が効くようになる
5. web-appを右クリック→Add attribute→versionで"version"属性追加
6. web-appを右クリック→Add Child→(略)→servletで"servlet"エレメント追加
 "servlet"エレメント以下に追加された項目を埋める
7. web-appを右クリック→Add Child→(略)→servlet-mappingで"servlet-mapping"エレメント追加
 "servlet-mapping"エレメント以下に追加された項目を埋める

●WARファイルのエクスポート
1. 作成したTomcatプロジェクトのフォルダアイコンを右クリック→Preferences選択
 →"Tomcat"の「WARエクスポート設定」タブにて.warのファイル名を設定
2. Tomcatプロジェクトのフォルダアイコンを右クリック→「Tomcatプロジェクト」→
 「プロジェクト設定に従いWARファイルを作成」を選択
3. できた.warをサーバーのTomcatのwebapps/の直下にコピー

●結果
・サーバー(Tomcat)側コード
TrichoiceServlet.java --- サーブレット本体
web.xml
QuizReader.java --- JSP用
listquiz.jsp
※ソースコード中のURL等は当方の開発用PCでのものです。

・クライアント側コード
Start.java --- メイン
MainPanel.java --- 3択画面
ServletAgent.java --- HTTP通信部分

・このサイトの起動ページ
3択アプレット
listquiz.jsp

・過去の関連エントリー
[1] [2] [3] [4]

See more ...

Posted at 16:44 in Java | WriteBacks (0)
WriteBacks

Feb 08, 2009

Maximaでの固有ベクトル計算について

筆者のMaxima 5.17+sbclの環境だと、標準のeigenパッケージを使って固有ベクトルを求めようとすると、エラーが出て求まらないことが多い。例えば、

load(eigen);
m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
eivects(m);

とやると、異様に時間がかかる上、
algsys failure: the eigenvector(s) for the1th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the2th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the3th eigenvalue will be missing.

というメッセージに続いて、虚数混じりのおぞましい値が固有値として出てきて、固有ベクトルは出力されない。そもそも固有値が間違っている。

固有ベクトルを英語でeigenvectorと言うことをたまたま覚えてたので、Maximaのヘルプでeigenvectorを検索して、出てきたのがeigenパッケージのeigenvectors()/eivects()だったので、それを使ったのだが、それが間違いだったようだ。Maximaのヘルプのeigenvectorsの項にも、固有値が正しく求められないことがあると、はっきり書いてある。対処方法も書いてあるが、とても面倒である。

インターネットで調べてみると、Maximaのメーリングリストのアーカイブに辿り着いた。どうやら2007年4月くらいからこの問題は原因も含めて知られているが、まだ解決していないらしい。

数式を含まない、純粋な数値行列なら、他の方法があることがわかった。
(1) eigens_by_jacobi()
標準関数として存在する。対称行列に限るという制約があるが、例えば共分散行列の固有ベクトルを求めるなら何の問題も無い。先に言ってくれという感じである。

(%i1) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i2) eigens_by_jacobi(m);
(%o2) [[- .7932213800442011, 1.180843735860609, .5243776441835921],
[ .6946378770918785 .5975383492309946 .4005323219261915 ]
[ ]
[ - .2001963037214474 .6953728293011738 - .6902014693159894 ]]
[ ]
[ - .6909411405362675 .3992549930407618 .6026572747810012 ]

eivects()と同様、固有値の行列と固有ベクトルの行列が一緒に返される。

(2) LAPACKパッケージのdgeev()
昔からある最強の線形演算ライプラリであるLAPACKのLisp版がMaximaに同梱されている。最初に使用する時にコンパイルがなされるが、これが環境依存であり、Lispの実行環境によってはコンパイルに失敗してしまうらしい。無事にコンパイルに成功すると、dgeev()で最強の計算結果が得られる。

(%i1) load(lapack);
(%i2) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i3) dgeev(m,true,true);
(%o3) [[- .7932213800442016, 1.180843735860611, .5243776441835922],
[ .6946378770918787 - .5975383492309948 .4005323219261915 ]
[ ]
[ - .2001963037214475 - .6953728293011737 - .6902014693159894 ],
[ ]
[ - .6909411405362672 - .3992549930407617 .6026572747810011 ]

この例は対称行列だが、もちろん非対称の行列でも問題ない。

筆者はgclとsbclで試したが、sbclではコンパイルが終了しなかったが、暴走したMaximaを強制終了すると、次からはコンパイル無しでロードできるようになった。(テスト段階で暴走していた?)

なお、root権限でインストールしたMaximaのLAPACKパッケージをコンパイルするには、root権限が必要である。(UNIX系のOSなら1回rootでMaximaを起動してload(lapack);とやれば良い)

See more ...

Posted at 23:19 in 数学 | WriteBacks (2)
WriteBacks

eigen_by_jacobi関数、役に立ちました。 ありがとうございます。 一つ質問があるんですが、数式を含んだ状態で固有値を求めた際、虚数を含まないようにする命令はありますか? A:matrix([-z,-1,0],[-1,-z/2,-0.25],[0,-0.25,-z/2]); B:eigenvalues(A); として固有値を求めると、対称行列なのに虚数が表われて、よく分からない数式になってしまいます。

Posted by yuji at 12/07/2010 12:51:00 PM

eigenvalues()はeigenパッケージにあり、まだeigenvectors()と共にバグがありますので、不正な答えが出力されることがあります。 また、3x3の大きさの行列の固有値は3次方程式の解であり、3次方程式の解の公式は(実数解でも一時的に)虚数が含まれますので、行列に変数が残っていると、実数解しかあり得ないとわかってても、虚数抜きで表現するのは困難だと思われます。(ちなみにMathematicaだと、こういう場合は「カルダノの公式」が適用されて、やはり実数解が虚数を使って表現されるようです) 固有値は固有方程式の解であることを使って、Maximaで solve(determinant(A-diagmatrix(3,x))=0,x); とすると、虚数は入りますが、正しい答えが得られます。 もし、zに何か実数が代入されてから計算されても良いのであれば、 A(z):=matrix([-z,-1,0],[-1,-z/2,-0.25],[0,-0.25,-z/2]); B(z):=eigens_by_jacobi(A(z)); としてから B(3) などとすると、実数表示の解が得られます。(ここでも残念ながら、eigenvalues()を使うと虚数を含む解の公式表示になってしまうようです)

Posted by ynomura at 12/07/2010 10:59:57 PM

Feb 06, 2009

統計学復習メモ8: 回帰直線と主成分分析

高校生の時に、必要も無いのに関数電卓を買った覚えがある。物理の先生に、理系ならずっと使うから早目に買っといて損は無いとか言われたからだと思うが、結局大学で物理実験をやるまでは遊びでしか使わなかった。それでも、元々電卓が好きだった私は、毎日毎日、関数電卓のプラスチックのボタンをカチャカチャと押して遊んでいた。
その電卓には、「一次回帰直線」の係数を求める機能があった。学校では習わなかったが、説明書にy=a+bxのaとbを求める機能と書いてあったので、意味は分かった。来る日も来る日も自室に逃げ込んで電卓と戯れてる内に、「一次回帰直線」がとても身近なものになった。

そんなに一次回帰直線に慣れ親しんだから、主成分分析には敏感に反応した。というのも、主成分分析は繰り返し回帰直線を求めることに似ていると思うからである。

一次回帰直線とは、パラメトリックな2変量の関係を直線で近似するものなので、2次元のデータを1次元で近似することに相当する。2変量をx,yとすると、標本(x1,y1)〜(xn,yn)の関係を直線y=a+bx(a,bは定数)で近似するので、その直線との誤差を無視して(x,y)=(0,a)+(1,b)tと置くと、2変量の標本(xi,yi)が1変量tiで表せるのである。
主成分分析の第一主成分だけを使うことを考えると、多変量の関係を1次元で近似することになるから、k次元の標本(x1i,x2i,…,xki)を(C1,C2,…,Ck)+(p1,p2,…,pk)tiに近似すると、k個の変量を持つ標本(x1i,x2i,…,xki)が1変量tiだけで表せる。(C1,C2,…,Ck)は定数ベクトルで、(p1,p2,…,pk)が第一主成分の単位ベクトルである。第二主成分の単位ベクトルを(q1,q2,…,qk)としてこれも使うとすると、(x1i,x2i,…,xki)≒(C1,C2,…,Ck)+(p1,p2,…,pk)si+(q1,q2,…,qk)tiという感じに、k変量の標本を2変量で近似することになる。

一次回帰直線の最小2乗法による求め方をおさらいする。標本(xi,yi)と直線y=a+bxとの誤差はyi-a-bxiであるが、この誤差の2乗和を最小にするようにa,bを決めたい。誤差の2乗和をEとすると、E=Σ(yi-a-bxi)2で、このEはaの2次関数であり、bの2次関数でもあるから、Eをaで偏微分したものが0になるaと、Eをbで偏微分したものが0になるbを求めれば良い。
\frac{\partial E}{\partial a}=0, \frac{\partial E}{\partial b}=0
これを解くと、
b\frac{\sum xy - \frac{1}{n}\sum x\sum y}{\sum x^2-\frac{1}{n}(\sum x)^2},a=\bar{y} - b\bar{x}
が得られる。bの分子はxとyの共分散、分母はxの分散なので、xの分散をVar(X)、xとyの共分散をCov(X,Y)と表すと、b=Cov(X,Y)/Var(X)とも書ける。実はシンプルである。

念のため、試しにやってみて確認する。例によってMaximaを使う。サンプルデータをdata1とし、

data1: read_nested_list("data1");
datam: read_matrix("data1");
load(descriptive);
b: cov(datam)[1][2] / var(datam)[1];
a: mean(datam)[2] - b * mean(datam)[1];

とすると、b=-.8338664607943058、a=5.848350705139143と求まる。そこで
plot2d([[discrete, data1], a+b*x],[x,0,2],[style,[points],[lines]]);

とすると、次のグラフが作られる。

主成分分析についても、詳細を省いて、とりあえずサンプルデータ(data2)を使って図示してみる。
下のように、標本が原点を中心に楕円盤状に散らばっているとする。(どれかの画像をクリックするとgifアニメが開きます)



少しずつ回転させてみたが、わかりづらいので、大体の分布の形を枠線で表示する。


第一主成分の単位ベクトルは最も誤差の和を小さくする直線の方向、すなわち楕円の長径方向の長さ1のベクトルになる。次の図の青線が第一主成分(をスカラー倍したもの)である。


標本から第一主成分を抜く(第一主成分を0にする)と、下の図のようになる。(回転方向は変更)


つまり、第一主成分と垂直な平面への射影である。元が楕円盤状だったのでこれも楕円状であり、この楕円を最も近似した直線、すなわち長径方向が第二主成分の向きになる。(下図緑線)


3次元なので、第二主成分も抜くと、直線状になる。その直線は第一主成分、第二主成分と直交し、その直線上に第三主成分がある。

空間が何次元でも同様である。最も誤差が小さくなる近似直線の方向の大きさ1のベクトルが第一主成分の単位ベクトルであり、その方向のベクトル成分を抜いた(超平面に射影した)ものの近似直線の方向の大きさ1のベクトルが第二主成分の単位ベクトルであり、その方向のベクトル成分も抜いて…を繰り返すと、全ての主成分が出てくる。

そんな反復計算をせずに主成分分析する方法は、次回まとめる。

See more ...

Posted at 22:30 in 数学 | WriteBacks (0)
WriteBacks

Jan 24, 2009

統計学復習メモ7: t検定で独立性の検定

前のエントリーではカイ2乗検定を用いてノンパラメトリック(定性的)な独立性の検定を試したが、次はパラメトリック(定量的)な場合をやってみる。独立かどうかを確かめたい2つの属性が定量的な値を持つパラメーター、C言語で言うとenum型やchar型でなくint型やfloat型の変数である場合に、その2変数間に相関があるかどうかを検定するものである。

なお、統計学では独立、相関という言葉はノンパラメトリックかパラメトリックかで使い分けるのが通常なのかも知れないが、ここでは区別せず、独立である/ないと相関がない/あるは同じ意味で使う。

まず、実際に例題でやってみる。

ある事業分野において、過去に優秀なプログラマーが書いた、読んで理解し易いしその後の仕様変更にも強いと評判の良い、模範的なC言語のソースコード一式があったとする。ソースファイル中のコード量(文の数や式の数、a=0;やb=foo(a);は1つ、while(式)はwhileと式とで2つ、for(式1;式2;式3)は4つ等)とコメントの量(/*~*/内の文字数)の関係が以下であった時、コード量とコメント量との間に相関がある無いかを検定したいとする。













ファイル名コード量コメント量
gentoo.h88234
emperor.h98413
macaroni.h485687
rockhopper.h285438
fairy.h184194
magellan.h487174
humboldt.h818344
king.h69139
adelie.h306244
cape.h219222












ファイル名コード量コメント量
gentoo.c5191876
emperor.c5152995
macaroni.c78517
rockhopper.c3131839
fairy.c3922142
magellan.c3691783
humboldt.c3452762
king.c73489
adelie.c7292202
cape.c4452309

統計学で相関と言えば、まず思い付くのは相関係数である。とりあえず、ヘッダファイルとソースファイルに分けて、コード量とコメント量の相関係数を求めてみる。
よく使われる相関係数は、共分散÷それぞれの分散の幾何平均、というもので、式で書くと、2つのパラメーターを持つ標本が(x1,y1)~(xn,yn)の時、
r_{xy}=¥frac {¥sum_{i=1}^{n}(x_i-¥bar{x})(y_i-¥bar{y})} { ¥sqrt{ ¥sum_{i=1}^{n}(x_i-¥bar{x})^2 ¥sum_{i=1}^{n}(y_i-¥bar{y})^2 }}
と表されるものである。(ピアソンの積率相関係数という名前が付いているらしい。)
例によってMaximaに計算させると、ヘッダファイルについてはr≒0.327、ソースファイルについてはr≒0.746であった。

さて、この値をどう考えれば良いのだろう。2つのパラメーターに十分に直線関係があるかどうかの基準として使われるものの1つに、相関係数が0.99以上(強い正の相関)または-0.99以下(強い負の相関)か、というのがある。検量線などの回帰直線の適切さの基準として使われることがあるようだが、そのように基準を決める自由がある場合ならそれでいいだろうが、一般的には(一般論としては)0.99という数字に根拠が無い。では相関係数がどういう値なら何なのだろうか。

そこで、検定統計量
T=¥frac{r_{xy}¥sqrt{n-2}}{¥sqrt{1-r_{xy}^2}}
が自由度n-2のt分布に従うという定理を使うのが、統計学の尊い教えである。

仮説を「相関がない」とし、有意水準を0.05(信頼度を95%)とすると、自由度n-2のt分布で右側(tが大きい方)の棄却域はt≧2.3くらいである。それに対し、ヘッダファイルについてはT≒0.978、ソースファイルについてはT≒3.17なので、ヘッダファイルについては仮説が棄却されず、コード量とコメント量の間に相関があるとは言えない、ソースファイルについては仮説が棄却され、コード量とコメント量の間に(正の)相関がある、ということになる。

結局、この方法では、相関係数はいくらぐらいだと相関があると言えるのだろうか。
まず、相関係数rと検定統計量Tの関係を見ておく。n=5,10,15,20の時のTのグラフは、次のようになる。
r-T
単調増加、点対称である。
Tがt分布に従うと、rの分布は次のグラフのようになる。
P(r)
当たり前だが、rが-1(負の相関)や+1(正の相関)に近いのは稀だということになっていることがわかる。
Tがt分布の両側5%の境目になるrの絶対値は、以下の値である。

nr
50.878
100.632
150.514
200.444

相関係数の絶対値がこれらより大きければ、相関があるということになる。
思ったより緩い基準である。

See more ...

Posted at 19:30 in 数学 | WriteBacks (1)
WriteBacks

I think you've just cpatuerd the answer perfectly[link]

Posted by Cah at 04/08/2012 03:35:04 AM

Jan 18, 2009

FMS(JAVA PRESS Vol.11より)

JAVA Press VOL.11のp.92から始まる「Java2で作るネットワークプログラム」という記事のコピーが手元にあった。いつかJavaの勉強のために読もうと思っていた、約9年前に発行された雑誌のコピーである。UDPを使って相手を探し、TCPを使ってメッセージやファイルをやり取りする、簡単なチャットアプリのTCP/UDP通信の部分のソースコードが掲載され、解説されている。画面イメージも載っている。"File&Message Send"の頭文字を取って"FMS"という名前のアプリだそうだ。当時はGUIの部分と合わせてフルセットコードがJAVA PRESSのWebサイトにあったらしいが、今は入手できないと思う。

掲載されているのは8つのクラスの内5つのクラスのソースコードで、合わせて300行程度である。画面イメージから読み取れるGUIはシンプルなものである。似たようなGUIさえ作れば動くんだろう、と思って作ってみたら、4倍近いコードに膨れ上がってしまった。記事に載っていたコードも、1台のマシンでもテストできるよう、動的にポートを切り替えるように改造したりしたので、結局大半のコードに手が入ってしまった。
折角なので、ここで公開することにする。

FMSアプリ(改)のソースコードとクラスファイル一式(JARファイル)

・画面イメージ

●起動方法
まず、FMSアプリを動かすそれぞれのマシンにfms.jarをダウンロードして、以下のコマンドで展開してください。(JREはそれぞれのマシンにインストール済だとします。)

jar xvf fms.jar

次に、接続先として検索するIPアドレスのリスト(seekips.txt)とユーザ名のリスト(uesrs.txt)を編集してください。ブロードキャストアドレスも指定可能です。ユーザ名は、名前とキーワード(2つ合わせて1つのユーザ名となる)をコロンで区切ってください。
・seekips.txtの例
192.168.0.255
172.1.1.1

・users.txtの例
John:America
Yamada:Hanako
Penguin:Antarctica
Whitebear:Arctica

次に、コマンドラインから下記のjavaコマンドで起動してください。
java fms.Run (自分の名前) (自分のキーワード)

1つのマシンで複数起動してもOKです。

●使用方法
アプリを起動すると、seekips.txtに書かれたIPアドレスとローカルマシンの別ポートに対してUDPのパケットが送信され、相手から返信があると、ウィンドウの左上の領域に相手が表示されます。但し、相手の名前が自分のusers.txtに無い場合は、表示されません(自分の名前が相手のusers.txtにある場合は、相手のウィンドウには自分の名前が表示されます)。
相手から検索があるとユーザ名の後ろが[COMING]となり、相手から検索の応答があると[ONLINE]となります。

[COMING]または[ONLINE]の相手が見つかったら、リスト上の接続相手を選択して、メニューバーの"Users"→"Connect"で接続します。接続に成功すると、相手の名前の後ろが[CONNECTED]になります。

接続に成功したら、メッセージまたはファイルを送る相手を選択し、ウィンドウ下段の入力フォームにメッセージまたはローカルマシン上のファイル名を入力して、"Msg Send"ボタンまたは"File Send"ボタンを押してください。送信に成功すると、相手側のウィンドウにメッセージが表示され、ファイルの場合は相手側のローカルマシンに同じファイル名で保存されます。

●ソースコードの説明
Seeker.java
相手を探したり自分の状態を伝えたりするためのUDPパケットを発信するクラス。
オリジナルのコードに動的ポート対応を加えている。
Watcher.java
UDPパケットを受信するクラス。
オリジナルのコードではユーザ情報の管理を別クラス(Users)にしていたらしいが、想像がつかなかったので、ここで管理するようにした。
SendPack.java
TCPで送信する1回分の情報を持つクラス。
オリジナルのコードから変更無し。
Connector.java
1つのTCP接続のソケットと入出力ストリームを保持するクラス。
オリジナルのコードに動的ポート対応を加えた。
Accepter.java
サーバーとしてTCPの接続要求を待つ(acceptする)クラス。接続要求が来る度にConnectorのインスタンスを作成する。
オリジナルのコードから変更無し。
Run.java
FMSアプリを起動するクラス。Seeker, Watcher, Accepter,MainPanelのインスタンスを作る。
GUIとしては、SwingのJFrameのインスタンスを用意して、メニューバー付きのウィンドウを作る。
全て自作。
MainPanel.java
ウィンドウ内の各部品を作り、それぞれに対する操作や要求を処理するクラス。実質のメイン処理。
全て自作。

より詳しく書くとJAVA PRESSの記事と被るので、説明はここまでとする。
なお、Run.JavaとMainPanel.JavaはEclipse+VisualEditor(GUIエディタ)で作ったので、他の環境のEclipse+VEでもGUIが読み込めるはず。

See more ...

Posted at 11:53 in Java | WriteBacks (0)
WriteBacks

Jan 04, 2009

数独解きアプレット

一時期、数独にはまったことがあり、腕には覚えがある身であり、一般紙に載る程度の問題は朝飯前であるはずだったのだが、年末年始に帰省して暇を持て余して、ある新聞の正月版に載っていた、懸賞問題の数独を暇潰しに解こうとしたら、解けなかった。
数字を埋め始めた直後から決して順調では無かったのだが、5分くらいかかって1/3程度埋めた後、10分以上1つも埋められなかった。問題が難しいというよりは、解無しな感触だった。まさかそれなりに発行部数が多い新聞で不完全な問題は出るまい、と思って粘っていたのだが、諦めた。

しかし、その問題のことが気になって気になって仕方が無かった。
数独ソルバーアプレットのページへのリンク
ソースコード

See more ...

Posted at 18:13 in Java | WriteBacks (0)
WriteBacks