とーけい部

現在、このブログは工事中です。

統計学 統計検定 数学

統計検定1級合格したので感想でも書いとく

2021/11/21に行われた統計検定1級(数理・応用)を受けました。数理の方は合否怪しかった(計算ミスが判明したため)んですけど、まあ両方とも受かっててよかったなって感じです。

試験対策

前日に過去問をやりました。2017年度の数理を1周し、応用は過去の問題をぺらぺら眺めて終わりにしました。
試験当日は久保川『現代数理統計学の基礎』で忘れてるところをザーッと、特に検定周りは普段使わないので知識を確認しておきました。
応用は理工学を選択していましたが、過去問見た感じしっかり勉強しとかないと無理そうだな~と思ってたんでほぼ諦め状態。とりあえず数理が終わって応用までの空き時間に『点過程の時系列解析』を読んでポアソン点過程の復習をしておきました。この本おすすめです。

試験中

数理は1,2,4を選択。思ってたより簡単で1,2は完答、4は(4)まで進めて(5)は手を付けられずに終了。
この時点では「受かったな~」という感じでしたが、昼休みにお茶飲みながら問題眺めてたら大問1の(4)で計算ミスを発見。さらに家に帰ってから大問1の(3)もミスってることが判明したのでガチで合否が怪しくなってしまいました。残念。
応用は1,3,5を選択。1は記述中心で、知ってる内容ばかりだったので普通にできました。3は計算が多め。最後のギブスサンプリングの問題が何を答えさせたいのかがいまいち不明瞭だったのを除けば、計算ミスなければできたかなという印象。5はベイズの定理を使って、正規分布表とにらめっこしながらひたすら計算していくタイプの問題でした。過去問見ていた感じだと応用の方が難しい(というか応用とは名ばかりで、数理.ver2みたいな問題ばかり)かなって思ってたんで、傾向が今までと違い過ぎてビビりました。結果として合格できたんで良かったです。

結果

12/20に結果発表があったらしい。翌日に「そろそろ結果発表かな~」と思って確認したところ出てましたね。数理は落ちてるかなと思ってたんですけど受かってました。よかったよかった。
あまり統計学の勉強やったことないよって人なら、数理・応用の両合格には大体100~400時間程度の勉強時間が必要かなと思います(最短ルートで100時間程度かな?)。大学1年レベルの数学力がある前提でこのくらいの時間が必要だと思うんで、大学受験で数学使わなかった!って人がとるのはかなり難しい資格と言えると思います。

確率を意識してマキシマム召喚を狙おう!(遊戯王ラッシュデュエル)

※この記事はサイバースデッキを中心に考察していますが、基本的な考え方はどのマキシマムデッキでも同じです。最後にダブルマキシマムが揃う確率にも触れておきます。

マキシマム召喚

レベル10の真ん中のパーツと左右の[L] [R]、合計3つのパーツを手札に揃えて行うのがマキシマム召喚ですよね。

f:id:hasegawa-yuta289:20211211224942p:plain

僕はアニメ遊戯王SEVENSでネイル君が使う天帝龍樹ユグドラゴが大好きなので大会でも大体はユグドラゴ軸サイバースを使用しているのですが、かなりピーキーでキツイです。大体2回戦で負けます。すまねえユグドラゴ、俺が、俺が弱いばかりに……
まあ仕方ないですよね。普通のデッキと違ってマキシマム召喚に特化したデッキはマキシマム召喚できるまでクソ雑魚ですし、運が悪いとパーツがデッキの下の方に眠っていてもうどうしようもないとかザラです。

しかし! 最近サイバースデッキにも期待の新星(ニュービー)がやってきました。アトラシュート・ハイドロン君です。

f:id:hasegawa-yuta289:20211211225139p:plain

これにより、ユグドラゴが揃うまではハイドロン使って粘るとか、逆にユグドラゴ突破された後でハイドロン使って巻き返すみたいな動き方が可能になりました。幻竜の場合はガントリードラゴン→ビルドドラゴンみたいな動き方をすることでマキシマムが揃うまでの隙を補うことができると思います。

マキシマム召喚ができる確率

ユグドラゴとハイドロン、二つの勝ち筋があるのは良いことなのですが、例えば『天の加護』などで手札のハイドロンを捨てたい場合があったり、他には『シードクロトロン・ブラッセルン』で墓地に落ちたユグドラゴパーツを回収するかハイドロンを回収するかで迷うことなどがあると思います。
こういった二つの戦略が考えられる場面でマキシマムパーツが揃う確率を考慮して自分の行動を選択することは非常に有意義だと言えますから、マキシマムパーツ3種類が揃う確率を調べてみたいと思います。
統計解析ソフトRを用いてモンテカルロ計算によりマキシマムパーツの『3種類目が存在する位置(確率)』を調べると、以下のようになりました。

f:id:hasegawa-yuta289:20211211231717p:plain

累積確率は以下の通りです。

デッキの上からn枚目 累積確率
3 0.002764
4 0.010212
5 0.023053
6 0.042244
7 0.067663
8 0.099404
9 0.136772
10 0.179053
11 0.225332
12 0.274372
13 0.325929
14 0.379198
15 0.432719
16 0.485926
17 0.538515
18 0.589217
19 0.637713
20 0.683814
21 0.726323
22 0.766382
23 0.802874
24 0.835776
25 0.865415
26 0.891724
27 0.914408
28 0.933760
29 0.950086
30 0.963655
31 0.974593
32 0.983073
33 0.989410
34 0.993935
35 0.996937
36 0.998770
37 0.999697
38 1.000000

図示するとこんな感じです

f:id:hasegawa-yuta289:20211212002108p:plain

これを考えると、デッキの上から17枚目まででマキシマムパーツの3枚目を引き当てる確率が50%以上ですから、ここら辺を目標に掘り進めていけばよさそうですね。『七宝船』『天の加護』『強欲な壺』『天使の施し』などを積極的に採用したマキシマム召喚特化の構築なら、よほどひどい事故をしない限り2ターン目には15枚程度は掘り進められますから、何も考えずとにかく掘り進める方針でやっていけば2ターン目には40~50%程度の確率でマキシマム召喚できます。3ターン目にはさらに掘り進みますから、7,8割でマキシマム召喚できますね(そんなにドローだけやってて3ターン目まで生き残ってるのかは疑問ですが)。

更に、既にマキシマムパーツが2種類揃っている状況で次のドローが最後のパーツである確率は以下のようになっています(近似値です)。

現在デッキの上からn枚目 次のドローで揃う確率
2~4枚 約8%
5,6枚 約8.5%
7枚 約9%
8,9枚 約9.5%
10,11枚 約10%
12,13枚 約11%
14枚 約11.5%
15枚 12%
16枚 12.5%
17枚 約13%
18枚 約13.5%
19枚 約14.5%
20枚 15%
21枚 約15.5%
22枚 約16.5%
23枚 約17.5%
24枚 約18.5%
25枚 20%
26枚 約21.5%
27枚 約23%
28枚 25%
29枚 約27%
30枚 約30%
31枚 約33.5%
32枚 約37%
33枚 約42%
34枚 50%
35枚 60%
36枚 75%
37枚 100%

グラフでプロットするとこんな感じですね
(正確な値は3/iを計算すれば良いです・i = 38 , … , 3)

f:id:hasegawa-yuta289:20211212011307p:plain

10枚目以降は、毎ドローするたびに10%のガチャガチャを引き続ける感じでしょうか。
10%を"ワンチャンある"の分水嶺とすれば、10枚目程度まで掘り進んでいた状態でアトラシュートなどでどうにかするのが無理っぽい時に手札にあと最後のユグドラゴパーツが来ればなんとかなる……! という状況で奇跡のドローにかけるのは悪くない賭けかと思います。10%でなんとかなります。20枚まで掘り進んでおけば15%でなんとかなりますね。

逆に言えば90%程度で引けないので、「次のドローが最後のパーツだったら出せるんだけど、どうしよっかな~、アトラシュートにしようかな~」という状況では素直にアトラシュート出した方がいいですね。間違いないです。ハーディフェンス怖いなら棒立ちさせるとかの選択肢も考えられます。

期待値計算すると、上から17.3枚目にマキシマムの最後のパーツが眠っていることになり、これの標準偏差は6.9でしたので、ここから上下7枚程度の変動は割とよくあるということになりますね。サイバースと当たったときも、11枚目程度にはもうユグドラゴ飛んでくると思っておいた方が良いです。

パーツを一つ減らす

マキシマムを使っていると、事故率低減のために真ん中のパーツを2枚にしたくなる欲求に駆られます。この時、3枚が揃う確率がどの程度変動するかを確認してみましょう。
図だけ書いておきます。

f:id:hasegawa-yuta289:20211212004740p:plain

f:id:hasegawa-yuta289:20211212004752p:plain

平気は18.7枚目だったので、マキシマムパーツを1種類だけ2積みにすると、マキシマムが出るまでの枚数が平均して1.5枚程度遅くなることがわかります。平均だけ見ると大して変わりませんが、標準偏差を考えるとこれは7.3枚となり、安定性が低下していますね(そりゃそう)。アトラシュート軸で、ユグドラゴ出せそうなら出そうかな~程度のノリで構築を組むなら、ユグドラゴの真ん中は2枚だけで問題なさそうです。

ダブルマキシマム

マキシマムパーツを2種類、3枚ずつ突っ込んだダブルマキシマム構築に対して同様の計算をすると、以下のようになります。

f:id:hasegawa-yuta289:20211212005131p:plain

f:id:hasegawa-yuta289:20211212005145p:plain

平均13枚、標準偏差は5枚となり、遅くとも18枚目にはマキシマム召喚できることが分かります。
安定性も速度も格段に上がっているのがわかりますね。しかし、この値はあくまでもマキシマムパーツを全て手札で保存している場合に限りますので、現実にはダブルマキシマム構築の場合はどちらのマキシマムパーツを犠牲にするかの択が常に発生するので、この理論値ほどの速度ではマキシマム召喚できません。

待ち時間のパラドックス(数式無し)

とあるバス停にて

 今ここに、世にも奇妙なバス停があります。このバス停では、どの時間にバスが来るのかわかりません。ランダムにバスがやってくるのです。
 とんでもないバス停ですね、使い物になりません。
 さて、できるだけこんなバス停使いたくないところですが、割と高頻度でバスがやってくるおかげで、ある程度バス停で待っておけば、まあそこそこの時間待つだけで乗れることにしましょう。

 「バスが一度やってきてから、次のバスが来るまで」の時間間隔をストップウォッチを使って、「バスが合計で10回来るまで」測ってみたところ、以下のようなデータがとれました。
※とあるバスが到着した瞬間を0回目として、そこから計測開始したとしましょう。(当然ながら、0回目のバスの到着時間は0分00秒です)

set.seed(200)
tau <- rgamma(n = 10 , shape = 0.5 , scale = 8)
plot(cumsum(c(0,tau)) , 0:length(tau) , type = "s" , main = "「バスが来るまでの時間」と「バスが来た回数」",
     xlab = "時間(分)" , ylab = "バスが来た累積回数")
points(cumsum(c(0,tau)) , 0:length(tau),pch=16,cex=2)

f:id:hasegawa-yuta289:20210819123241p:plain  さて、この図によれば、バスが10回到着するまでにかかる時間は40分くらいなので、つまり4分に一回くらいは、バスが来るということになりますね。

 この「バスが一度到着してから次のバスが到着するまでの時間」のことをここでは「バス到着時間間隔」と呼ぶことにしましょう。あとで説明する「待ち時間」とは別の概念なので注意してください。

待ち時間のパラドックス

 あなたはこのバス停を利用したいと思っています。つまり、ふらっとこのバス停まで歩いてやってきて、次のバスがくるまで待つ……ということをしましょう。
 ここで、あなたは、先ほど確認したように「4分に一回くらいはバスがくる」ということを知っているとしましょう。

さて、問題です。

 あなたはどのくらいの時間このバス停で待たなければならないでしょうか?

フツーに素直に考えれば、この「歩いてバス停まで到着したときに次のバスが来るまで待たなきゃいけない時間(の平均)」は、「このバス停で降りてすぐに次のバスを待つまでの時間(の平均)」である4分よりは短いですよね。

答えは……実際に確かめてみましょうか。

 さて、先ほど確認したデータ上では、バスが10回到着するまでの累積時間(単位は「分」にしましょう)は以下のようになっています。

cumsum(tau)
> cumsum(tau)
 [1]  4.993815 10.672151 13.749107 14.014109 15.660872 19.262998 32.793644 35.170187 35.648799 41.419587

なので、あなたが歩いてやってきたとして、バス停に到着した時刻が0分~41.419分のどれかの値だということにしてみます。一様乱数を使って、あなたがバス停に到着した時刻を設定してみましょう。

set.seed(10)
TOUCHAKU <- runif(n = 1 , min = 0 , max = 41.419)
TOUCHAKU
> TOUCHAKU
[1] 21.01924

ということで、あなたがこのバス停に到着した時刻は21.0分ということになりました。次のバスが到着するまでの時刻まで待つと、その時間は32.8 - 21.0 = 11.8分でした。

あら、残念! バス到着時間間隔は平均して4分なのに、その倍以上の11.8分も待たなければなりませんでした!

運が悪かったですね……本当に、そうでしょうか?

 この、「バス停まで歩いてやってきて、到着してから次のバスがやってくるまでの時間」を「待ち時間」と呼ぶことにしましょう。あなたの今の待ち時間は約12分でした。バス到着時間間隔の平均は4分なので、その倍以上の時間を待ったわけですね。

 直感的に考えて、「バス到着時間間隔」は「待ち時間」よりも(平均的に)長くなるでしょう。しかし、実はそれが成立しない場合があります。今考えている例では、実はそのような状況を設定しています。

 これを「待ち時間のパラドックスと呼びます。

確かめてみる

 先ほどはバスが「10回」やってくるまで計測しましたが、もっとデータを増やすためにバスが「10000」回やってくるまでストップウォッチを使って気合で計測したとします。  

set.seed(1)
tau.10000 <- rgamma(n = 10000 , shape = 0.5 , scale = 8)
max(cumsum(tau.10000))
> max(cumsum(tau.10000))
[1] 39814.38

 バスが10000回到着するまでの合計時間は39184.38分でした。663時間、まる27日はかかってますね。頑張りました。

 ちなみに、バス到着時間間隔の平均は

mean(tau.10000)
> mean(tau.10000)
[1] 3.981438

4分だということがわかります。

 さて、あなたがふらっとこのバス停に到着した時刻を0分~39184.38分の間でランダムに取ることにしましょう。その「歩いてバス停にたどり着いた時刻」から「次のバスがやってくる」までの待ち時間を計算する、ということを仮想的に100000回行ってみます。

MACHIJIKAN <- c()
for(i in 1:100000){
  TOUCHAKU <- runif(n = 1 , min = 0 , max = max(cumsum(tau.10000)))
  MACHIJIKAN[i] <- cumsum(tau.10000)[(cumsum(tau.10000) - TOUCHAKU) > 0][1] - TOUCHAKU
}
mean(MACHIJIKAN)
> mean(MACHIJIKAN)
[1] 6.084652

こうすると、待ち時間の平均は6分で、明らかにバス到着時間間隔の平均である4分よりも大きいですね。待ち時間のパラドックスが起きていることがわかりました。

直感的な説明

 もう一度、先ほど計測したバスが10回来るまでのデータを見てみます。

plot(cumsum(c(0,tau)) , 0:length(tau) , type = "s" , main = "「バスが来るまでの時間」と「バスが来た回数」",
     xlab = "時間(分)" , ylab = "バスが来た累積回数")
points(cumsum(c(0,tau)) , 0:length(tau),pch=16,cex=2)

f:id:hasegawa-yuta289:20210819123241p:plain
 バスが一度到着してから、次のバスが到着するまでの時刻ですが、これに結構なばらつきがあるように見えますね。一度バスが到着したら即座に次のバスが来るか、一度バスが到着したら次のバスが来るまで結構な時間がかかるか、そのどちらかが起きる可能性が高そうです。つまり、バス到着時間間隔の分散が大きいわけですね。
 それら到着時間間隔の極端な長さ・極端な短さが互いに打ち消し合って、平均してみるとバス到着時間間隔は4分になっているわけです。

 ここで、ランダムな時間にバス停まで歩いてたどり着いたとして、たどり着いた時点がどのバス到着時間間隔に当てはまるか? を考えてみると、これは多分「バス到着時間間隔が長い間隔」内で到着する事の方が多いだろうと考えられますね。この図で言えば、20分~30分の間にバス停までたどり着いた場合は明らかに長時間待つ必要があります。

 このように、バス到着時間間隔の分散が(その平均に比べて)やたらと大きい場合は、待ち時間のパラドックスが発生してしまいます。

参考文献

 数式を使った具体的な説明は以下の本に載っているので、興味のある方は読んでみてください。全然話題を聞かないですけどめちゃくちゃ良い本だと思います。

カーネル密度推定

密度を推定する

 データが点の集まりとして手に入ったとしましょう(例えば身長・体重データとか喫煙歴とか)。それで、この点々の母集団の分布を考えるときは、例えばこの点々を基にしてQQプロット打ち込んでう~ん正規分布っぽいから平均と分散を推定するか! とかやっていくわけです。これがパラメトリックな密度の推定になります。
 さて、ノンパラメトリックな密度の推定を考えてみましょう。つまり、母集団の分布が属する"族"を(ほとんど)設定せず、ただ単に今手元にある点々から、それっぽい密度関数をつくりたい……そういう需要はあると思います。めっちゃあります。
 というわけで、今回はノンパラメトリックな密度推定法の、カーネル密度推定と呼ばれる手法について軽くまとめ&実装してみました。

カーネル密度推定

 1次元の場合を考えます。
 カーネル密度推定において、点x上での密度の推定値は、バンド幅をh、サンプルサイズをnとしたときに、

\displaystyle{
\hat f (x) = \frac{1}{nh} \sum _{i=1}^{n} K\left( \frac{x-X_i}{h}\right)
}

によって表されます。ここで、\displaystyle{
K()
}カーネル関数で、じゃあカーネル関数ってなんだよというと、

\displaystyle{
\int_{-\infty}^{\infty}K(t)dt = 1
}

なる関数のことです。
 バンド幅ってなんだよとなると思いますが、まあイメージとしてはヒストグラムのビン幅みたいなものです。このパラメータhのことをバンド幅って呼ぶんだなあぐらいの認識でも良いと思います。
 なお、通常、このカーネル関数は原点対称なものが使われ、例えば正規分布の密度関数なんかはガウシアンカーネルと呼ばれています。
 このカーネル関数を色々と設定することによって、ただのサンプル(点の集まり)から密度推定が可能になるのです。

 カーネル関数の定義より、この推定された密度関数がきちんとマイナス無限大から無限大で積分したときにその値が1になっているとわかりますね。

統計学カーネルというと文脈によって複数意味があるので注意。ここでいうカーネル関数というのはノンパラメトリック統計におけるカーネル関数のことです。ベイズ統計や機械学習を勉強したことある人なら実感したことあると思います。参考→カーネル(統計学) ちなみにマルコフ過程の推移核のことをマルコフカーネルとか確率カーネルとか言ったりもする(M(x,dy)みたいな表記)。 個人的に、カーネルって言われても何のこと言ってるのかよく分からないので正直この言葉使ってほしくないですね。

 カーネル関数の種類は多分上げてったらキリがないと思うんですけど、さっきのWikipediaの記事にいくつかまとまっているので気になって人は読んでみてください。

非心t分布の密度関数を図示したい

\displaystyle{
U \sim N(\lambda , 1) \\
V \sim \chi ^2(m) \\
UとVは独立
}

という状況を考えましょう。このとき、

\displaystyle{
\Large
\frac{U}{\sqrt{\frac{V}{m}}}
}

が従う分布のことを、非心度λ、自由度mの、非心t分布と呼びます。t(m,λ)で表すことにしましょうか。例えば仮説検定で母集団同士の平均を比べるときに平均に差があるような場合はt統計量がこの分布に従いますね。

 この定義から明らかに、非心度λが0である非心t分布、t(m,0)は、自由度mのt分布と一致することがわかります。
 さて、ではt(m,λ)の密度関数はどんな関数になっているかというと、以下のような関数になっています。(求め方は普通にt分布の密度求める時と一緒なので、見た目ほど難しくはないのですが、計算はめんどくさいです)

\displaystyle{
\LARGE
f(t\ ;m,\lambda)=\frac{exp(-\frac{\lambda ^2}{2})}{\sqrt{\pi m} \  \Gamma(\frac{m}{2}) } 
\sum _{j=0}^{\infty} \frac{\lambda ^j t^j \Gamma \left(\frac{m+j+1}{2}\right)2^\frac{j}{2}}{j\ ! \ m^{\frac{j}{2}} \ \left(1+\frac{t^2}{m}\right)^{\frac{m+j+1}{2}}}
}

なんだこれ!?

 さて、私はちょっと気になったんですけど、t(m,λ)って、λの値を動かしたときにどういった感じで形が変わっていくんでしょうか? 定義より明らかにλを動かすほど分布が横にスライドしていきそうなものですが……ほかには、expの係数から、λが0より離れていくほど、分布が潰されていきそうな感じがするでしょうか。まあ、ぱっと見じゃ、そのくらいのことしかわからないですね。

 ということで、非心t分布の密度関数を私は図示したいのです。お気楽に、お手軽に。
 この関数には無限和が含まれているのでどこかのタイミングで計算を打ち切って密度関数を書くしかないですが、もう一つ、密度関数を描く方法がありそうですね。

 この分布は式は複雑ですが、この分布に従う乱数の発生は比較的容易ですから、もう乱数を大量に発生させてしまって、そのサンプルから密度をノンパラメトリック推定してしまおう~! というのが、目標です。

 じゃあまず乱数を沢山作りましょう。自由度は10にします。非心度は今はとりあえず5にしておきます。後でいじります。

n <- 10000
m <- 10
lambda <- 5
U <- rnorm(n=n , mean = lambda ,sd = 1)
V <- rchisq(n=n , df=m)
RANSU <- U / sqrt(V/m)
hist(RANSU)

f:id:hasegawa-yuta289:20210814021612p:plain

 さて、これによってt(10,5)に従う乱数が10000個できました。
 ヒストグラムだけでも、まあ正直なんとなく密度関数の形が見えてきますけど……
 ここから、カーネル密度推定を、上でやったように行いましょう!

 バンド幅hは色々と設定できるので、とりあえずh=0.2にしておきます。
 カーネル関数にはガウシアンカーネルを利用しましょう。標準正規分布の密度関数なので、以下のような関数になっています。

\displaystyle{
K(t) = \frac{1}{\sqrt{2\pi}}exp \left(-\frac{t^2}{2}\right)
}

再掲:カーネル密度推定

\displaystyle{
\hat f (x) = \frac{1}{nh} \sum _{i=1}^{n} K\left( \frac{x-X_i}{h}\right)
}

これを参考にしてプログラミングしていきましょ~

h <- 0.2
f.hat.Gaussian <- function(x,h,X){
  n <- length(X)
  ATAI <- c()
  for(i in 1:length(x)){
  ATAI[i] <- sum(dnorm((x[i]-X)/h)) / (n*h)
  }
  return(ATAI)
}
hist(RANSU , freq=FALSE , main = paste("ガウシアンカーネルによる密度推定 h=" , h , sep="") , 
     xlab = paste("x (n=" , n , ")", sep=""))
curve(f.hat.Gaussian(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col="red" , lwd=2 , add=T)

f:id:hasegawa-yuta289:20210814021633p:plain

ということで、ガウシアンカーネルを使った非心t分布の密度関数の推定ができました。
 nを増やし、hを細かくしていけばより精緻な推定ができるでしょう。

 ちなみに、試しに三角カーネルを利用したときの密度推定は以下のようになります。Iは定義関数です。(その他の設定は全部同じ)

\displaystyle{
三角カーネル:K(t) = (1-|t|)I(t\in [-1 , 1])
}
h <- 0.2
f.hat.triangular <- function(x,h,X){
  n <- length(X)
  ATAI <- c()
  for(i in 1:length(x)){
    t <- (x[i] - X) / h
    ATAI[i] <- sum( (1 - abs(t)) * (abs(t)<1) ) / (n*h)
  }
  return(ATAI)
}
hist(RANSU , freq=FALSE , main = paste("三角カーネルによる密度推定 h=" , h , sep="") , 
     xlab = paste("x (n=" , n , ")", sep=""))
curve(f.hat.Gaussian(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col="red" , lwd=3 , add=T)
curve(f.hat.triangular(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col="green", lwd=2 , add=T)
legend("topright" , legend=c("ガウシアンカーネル","三角カーネル") , lty=1,col=c("red","green"),lwd=2)

f:id:hasegawa-yuta289:20210814021646p:plain

 まあ、大体同じ結果になりましたね。実用上はどっちを使っても問題なさそうですが、どうでしょう?

 さて、ここから本題です。t(m,λ)のλを色々動かしたときの、非心t分布の変化を見ていくことにします。
 バンド幅は0.05で、n=1000000にしましょう。自由度は10として、非心度を0,2,4,……,20と変化させます。

curve(dt(x,df=10),xlim=c(-2,24) , col="red" , lwd=2 , main="非心t分布(自由度10 , 非心度0~20)",ylab="密度")
n <- 1000000
m <- 10
h <- 0.05
for(lambda in seq(2,20,2)){
  U <- rnorm(n=n , mean = lambda ,sd = 1)
  V <- rchisq(n=n , df=m)
  RANSU <- U / sqrt(V/m)
  curve(f.hat.Gaussian(x = x , h = h , X = RANSU) , xlim = c(min(RANSU) , max(RANSU)) ,
      col=rgb(red=(20-lambda)/20 , blue=lambda/lambda , green=0) , lwd=2 , add=T)
}

f:id:hasegawa-yuta289:20210814021701p:plain

ということで、t(10,λ)の密度関数は、λを0から20まで動かすとこの図のようになることがわかりました。
 分布が右にスライドすることや、段々と分布が潰れていくことは、予想通りの結果になりましたね。

 ということで、今回はサンプルから密度を推定する、カーネル密度推定を行いました。ちなみに今回作ったRのコードをいちいち書かなくても、Rにはこれと同じことをやってくれるdensity()という関数があるのでそれを使えば簡単です。便利ですね。ちなみにバンド幅hの良い選び方とかもあるんで、気になった人は↓の本を読んだり調べたりしてみてください。(おわり)

参考

棄却サンプリング

棄却法

密度πをもつ確率分布からの乱数が発生したいとします。
πに対して、そのサポート(π(x)>0となるようなxの範囲のこと)を包むようなサポートを持っている密度qをもつ分布からの乱数が発生できるとします。
ここで、πのサポート上の任意のxについて、以下のような定数Cが存在するとします。

\displaystyle{
\frac{\pi (x)}{q(x)} \leq C
}

このような定数Cについて、以下のアルゴリズムにより発生される乱数xは、密度πをもつ分布に従っています。

[アルゴリズム:棄却法]
①qから乱数yを発生させる
②(0,1)の一様乱数uを発生させる
\displaystyle{
u \leq \frac{\pi(y)}{C q(y)}
} であれば、このyを乱数xとする。ステップ①に戻る。

ちなみに、棄却法に限らず、このような我々が乱数を発生させることができる(簡単な)分布qのことを、提案分布と呼びます。乱数が欲しい分布πのことを目的分布と呼ぶこともあります。

なお、棄却法による、πに従う乱数を1個発生させるために必要な、qから発生させる乱数の個数は平均C個ですので、このような値Cは可能な限り低くとると、効率が良いです。

ちなみに、πに対してqの関数の減衰が明らかに速いと、π(x)/q(x)を上から抑える数Cはいくらでも大きくなってしまうので、注意。

理論的な背景は、鎌谷『モンテカルロ統計計算』の58p~や、小澄『ベイズ計算統計学』の30p~を読むと詳しいので、読んでみてください。

本当かよ?

さて、今までの説明は本当でしょうか? 私が嘘をついている可能性がありますね。
ま、試しにちょっとやってみましょうか。

位置パラメータ0、尺度パラメータ1、形状パラメータ3のフレシェ分布を考えてみましょ~

\displaystyle{
\large
\pi (x) = 3 x^{(-3-1)}exp(-x^{-3})
}
curve(3*x^(-3-1)*exp(-x^(-3)) , xlim=c(0,7) , main = "フレシェ分布(形状パラメータ=3)")

f:id:hasegawa-yuta289:20210607190813p:plain

ちなみにこいつはタイプII極値分布と呼ばれる分布になっていて、極値統計学との関係性が深いのですが、まあ今はそこら辺の話は置いときましょう。

この分布からの乱数を考えてみます。

こいつは結構裾が厚いので、例えば同じく正の値のみをとる分布として指数分布を考えてみると、

\displaystyle{
\large
q(x) = \lambda e^{-\lambda x}
}

これについて、指数分布の減衰の方が圧倒的に早いので、π(x)/q(x)は、特に裾の方で以下のようになり、

curve(3*x^(-3-1)*exp(-x^(-3)) , xlim=c(5,15) , main = "フレシェ分布と指数分布" , col="red")
curve(exp(-x) , add=T , col="blue")
legend("topright" , legend = c("フレシェ分布" , "指数分布") , col=c("red" ,"blue") , lty=1,cex=1.5 , lwd=2)
curve(3*x^(-4)*exp(x-x^(-3)) , xlim=c(0,10) , main="比率")

f:id:hasegawa-yuta289:20210607190854p:plain

f:id:hasegawa-yuta289:20210607190930p:plain

これを上から抑えるような定数は存在しないので、棄却法は実行できませんね。

じゃあ、もう仕方ないんで、私たちが良く知っている最強に裾が厚いアイツを使いましょう。

(標準)コーシー分布を提案分布として利用します。

\displaystyle{
\large
q(x) = \frac{1}{\pi}\frac{1}{1+x^2}
}

これとフレシェ分布の密度との比率は以下のようになります。

curve(3*pi*(x^(-2)+x^(-4))*exp(-x^(-3)) , xlim=c(0,10) , main="比率")

f:id:hasegawa-yuta289:20210607190958p:plain

というわけで、コーシー分布を提案分布として、棄却法アルゴリズムを実行してみます。
図から見る限り、まあ大体7くらいで抑えられるみたいなんで、C=7.2として実行してみましょう。

res <- c(0)
N <- 10000
flag <- 0
for(i in 1:N){
  y <- rcauchy(n=1)
  u <- runif(1)
  if(u < (3*y^(-4)*exp(-y^(-3)))/(7.2*1/(pi*(y^2+1))) && y>0){
    res <- append(res , y)
  }else{
    flag <- flag+1
  }
}
res
flag

さて、手に入れた乱数resをプロットしてみましょう。

plot(density(res) , xlim = c())
curve(3*x^(-3-1)*exp(-x^(-3)) , col="red" , add=T)
legend("topright" , legend="フレシェ分布" , lty=1 , lwd=2 , col="red" , cex=1.7)

f:id:hasegawa-yuta289:20210607191028p:plain

このように、棄却法により手に入れた乱数が、フレシェ分布に従っていそうだと確認できました。

Quasiモンテカルロお試し:(超一様分布列による積分計算)

超一様分布列(LDS)と積分

超一様分布列というのは、凄い一様に分布してる列のことです(は?)。
これは英語で言うとLow-Descrepancy Sequence なので、略してLDSと呼ばれることもあるようです。

例えば、1/2 , 1/4 , 3/4 , 1/8 , 5/8 , 3/8 , ……
というような数列を第100項まで図示すると、以下のようになります。
f:id:hasegawa-yuta289:20210605231708p:plain

これに対して、ただ単に一様乱数を100個とっただけの数列を図示すると、以下のようになります。
f:id:hasegawa-yuta289:20210605231813p:plain

比べてみれば一目瞭然ですね、超一様分布列(LDS)の方が、ただの一様分布の列よりも、わざとらしくもっと"一様に"分布しているのがわかります。

こういった超一様分布列を利用することで、積分計算を通常のモンテカルロよりも早く終わらせることができる、可能性があります。便利ですねえ~
LDSを利用したモンテカルロ計算(いやまあ、もはやモンテカルロ計算ではないんですけど)のことを、Quasi -モンテカルロと呼ぶのですが、QMCと略すことにしましょう。
Quasiとは、"疑似"という意味らしいので、日本語訳すると疑似モンテカルロということになりますけど……超一様分布列は「決定論的なただの数列」なので、これを使った積分モンテカルロ計算と呼ぶとややこしいですね。モンテカルロ法は通常は(疑似)乱数からの積分計算を指しますからね。まあしょうがないですね、そういう用語なんで。
(追記:Quasi-モンテカルロの日本語訳は準モンテカルロらしいです。)

というわけで、今回はLDSを利用した積分をちょっとやってみましょう。

積分をやってみよう!

以下のような関数hを考えます。

\displaystyle{
\Large{
\begin{eqnarray}
f(x) &=& sin(10x)arctan(10x)^2 \\
g(y) &=& -cos(2\pi y)|log(\Gamma(3y))|^ {arcsin(y)} \\
h(x,y) &=& (f(x)+g(y))((x+y)/2+1/2)
\end{eqnarray}
}
}

関数hは、0<x<1 , 0<y<1 の範囲について、以下のような関数になっています。
(なんかうまくgif作れなくて、長時間眺めてると気持ち悪くなってきそうな図になりました。申し訳ない)

f:id:hasegawa-yuta289:20210605233104g:plain

これについて、0<x<1 , 0<y<1 での積分をやってみましょう。

超一様分布列には、Halton列を利用します。
計算結果は以下のようになりました。

f:id:hasegawa-yuta289:20210605233451p:plain

一様分布を利用した通常のモンテカルロよりも、Quasiモンテカルロの方が収束が速そうに見えますね。

今回の記事を再現するためのコードを、以下に置いておくので、試しにRで動かしてみてくださ~い。(おわり)

コード

## 設定
f <- function(x){
  a <- sin(10*x) * atan(10*x)^2
  return(a)
}
curve(f , xlim=c(-1.2,1.2),col="red",lwd=2)
abline(h=0,lwd=2)

g <- function(y){
  b <- -cos(2*pi*y) * abs(lgamma(3*y))^asin(y)
  return(b)
}
curve(g,col="blue",lwd=2)
abline(h=0,lwd=2)

h <- function(x,y){
  c <- (f(x)+g(y))*((x+y)/2+1/2)
  return(c)
}

x <- y <- seq(0,1,length=50)
z <- outer(x,y,FUN=h)
persp(x=x,y=y,z=z,
      theta = -40,
      phi   = 20)

## 超一様分布列の利用

# 2,3進数変換  
trans2 <- function(n){
  if(n <= 0){
    return(NULL)
  }else{
    return(append(Recall(n%/%2), n%%2))
  }
}
trans2(1)

trans3 <- function(n){
  if(n <= 0){
    return(NULL)
  }else{
    return(append(Recall(n%/%3), n%%3))
  }
}

## Halton列の作成  
H2 <- function(N){
  a <- numeric()
  for(i in 1:N){
    b <- rev(trans2(i))
    a[i] <- sum(2^(-1*(1:length(b))) *b) 
  }
  return(a)
}

plot(H2(100),main="超一様分布列")
plot(runif(100) , main="ただの一様分布列")

H3 <- function(N){
  a <- numeric()
  for(i in 1:N){
    b <- rev(trans3(i))
    a[i] <- sum(3^(-1*(1:length(b))) *b) 
  }
  return(a)
}


plot(runif(1000),runif(1000),main="(疑似)乱数列")
plot(H2(1000) , H3(1000) , main="Halton列")

N <- 10000
res1 <- numeric()
res2 <- numeric()
X <- runif(N)
Y <- runif(N)
for(i in 1:N){
  res1[i] <- mean(h(X[1:i] , Y[1:i]))
}
plot(res1,type="l")

X <- H2(N)
Y <- H3(N)
for(i in 1:N){
  res2[i] <- mean(h(X[1:i] , Y[1:i]))
}
plot(res2,type="l")

plot(res1,type="l",col="red",ylim=c(-0.4,0.4),main="通常のモンテカルロとQuasiモンテカルロの比較")
lines(res2,col=adjustcolor(alpha=0.7,col="blue"))
legend("topright",col=c("red","blue"),legend=c("MC","QMC"),lty=1,lwd=2,cex=1.5)

参考文献

デルタ関数のフーリエ変換は「1」←ほんとかよ?

雑談

今書いてる記事のシリーズに「ハミルトン本を読もう!」っていうのがあるんですけど(2021/4/14注記:いろいろ忙しくて更新できてないです。再開のめどがたつまで非公開にしておきます)これは時系列分析の本なのでスペクトル分析のところでどうしてもフーリエ変換が必要になるため、最近フーリエ変換を勉強してます。
理工系の人ならそんなん常識やろ~wってなるかもしれませんが、ぼく文系なんで今まで触れる機会が無くて、全く知りませんでした!w すいません許してください何でもしますから!

まあ、冗談はさておき、結構面白いですねこれ。全ての周期関数(は言い過ぎだけど、病気みたいな関数を除けば)はサイン波とコサイン波の重ね合わせで表現できるっていう。これを利用することで、時系列データについてその周波数領域での考察を行うことができます。まあ決定論的なアレから確率論的なアレへの拡張は必要ですが。
で、まあこのフーリエ変換なんですが、ディラックデルタ関数(x=0でボーン!と上に跳ね、それ以外の点では0になるやつ)をフーリエ変換すると「1」という定数になる らしいです。微分方程式の講義でまあラプラス変換はやってたんで結果自体は知ってるんですけど、フーリエ変換の意味を改めて考えると、「え、ほんとか?」て思ったんで、ちょっと調べてみます。当たり前ですが、結果自体はきちんと成立します。

フーリエ変換

以下を関数f(x)の角周波数へのフーリエ変換F(ω)とよぶ。

\displaystyle{
\large{
F(\omega) = \int_{- \infty}^{\infty} f(x) e^{-i\omega x}dx
}
}

そして以下を、F(ω)のフーリエ逆変換とよぶ。

\displaystyle{
\large{
f(x) = \frac{1}{2\pi}\int_{- \infty}^{\infty} F(\omega) e^{i\omega x}d\omega
}
}

このフーリエ逆変換の式の意味:exp(iωx)が周期(角周波数)の同じサイン波とコサイン波の和であることから、異なる周期のサイン波コサイン波の重ね合わせで関数f(x)が表現されており、その重ね合わせのうち、角周波数ωに対応するサイン波コサイン波の、f(x)に対する強さがF(ω)という関数(その波にかかっている係数)で表現されていることがわかる。
F(ω)が大きければ、f(x)のうち、このωに対応しているサイン波コサイン波のパワー!がでかいってことですね。

ディラックデルタ関数

これは超関数とか呼ばれてる(超関数で調べてもよくわからなかった)関数で、まあ非常にざ~っくりと書くと、以下のような関数になっている。

\displaystyle{
\large{
\delta (x) =
\begin{cases}
\infty & (x=0) \\
0 & (x \neq 0)
\end{cases}
}
}

図で書くと、以下のような関数ですね。これは近似ですが。

pe <- function(e,x){
  a <- ifelse(-e<x&x<e,1/(2*e),0)
  return(a)
}
curve(pe(0.00001,x),xlim=c(-0.001,0.001),main="デルタ関数っぽいもの")

f:id:hasegawa-yuta289:20201220215046p:plain

確かめてみる。

以下のように、exp(iωx)×dω を、ωをめ~~~~~っちゃ細かくとって、それを全部係数「1」で足し合わせてみます。
め~~~~~っちゃ細かくとるんですけど、その振る舞いを分かりやすくするために、だんだん細かくしていきましょう。

f <- function(x){
  len <- 100
  omega <- seq(-1000,1000,length=len) #-1000~1000までをlen分割する
  d_omega <- omega[2]-omega[1]
  a <- c()
  for(i in 1:length(x)){
    a[i] <- sum( 1 * exp( 1i *omega*x[i]) * d_omega )
  }
  return(Re(a))
}
curve(f(x),xlim=c(-10,10),main="フーリエ逆変換後の関数(ω=-1000から1000までを100分割)")

ω=-1000から1000までを100分割して、exp(iωx)を全部足します。これがf(x)(に2πをカケたもの)の近似になってます。
図にすると以下。
f:id:hasegawa-yuta289:20201220215503p:plain

おいおい全然デルタ関数じゃないじゃないっすか~、本当にデルタ関数になるんですかぁ~?って感じですが、もっと細かくしていきましょうか。100分割じゃ足りなかったか。

1000分割するとこんな感じ。
f:id:hasegawa-yuta289:20201220215627p:plain

おや、だいぶ近づきましたね。

まあ、もういいですかね。飽きたでしょ。1000000分割するとこんな感じ。

f <- function(x){
  len <- 1000000
  omega <- seq(-1000,1000,length=len) #-1000~1000までをlen分割する
  d_omega <- omega[2]-omega[1]
  a <- c()
  for(i in 1:length(x)){
    a[i] <- sum( 1 * exp( 1i *omega*x[i]) * d_omega )
  }
  return(Re(a))
}
curve(f(x),xlim=c(-10,10),main="フーリエ逆変換後の関数(ω=-1000から1000までを1000000分割)")

f:id:hasegawa-yuta289:20201220215757p:plain

これはもうさすがにデルタ関数に収束するとみてもよさそうですかね?

というわけで、デルタ関数フーリエ変換するとまじで定数「1」になるっぽいことが確認できました。
周期がほんの少しずつ異なるサイン波とコサイン波を足しまくるとx=0の点で爆発してそれ以外では0になる関数(つまりデルタ関数)になるっていうのは、ちょっと面白いですね。
まあよく考えたら、どの周期のコサイン波でも x = 0のとき cos 0 = 1だから必ず1になるし、あとまあ、それ以外の点では互いに打ち消しあって0に近づくってことなんでしょうね。当たり前な気がしてきた。(おわり)