スクフェスで学ぶ確率・統計（ベイズ推定編）
今回はスクフェスのイベント「メドレーフェスティバル（以下メドフェス）」を題材にして，統計の基礎を学ぼうという企画です．ベイズ統計の対象としては，離散値を扱うので非常に単純ですが勉強にはちょうどいいでしょう．ちなみに私は統計に関しては初歩的な知識しかなく，ベイズ統計は今回初めて触りました．以下の説明でも間違っているところがあるかもしれない．
今回のお題
メドフェスで金報酬，銀報酬，銅報酬が出る確率p_g, p_s, p_cをそれぞれ見積もりたい．n番目の報酬E_nを収集した時点でのそれぞれの確率を推定せよ．
理論
今回のケースでは確率分布の形，即ち母数に対する推定を行いたいため，次に示す母数分布に対するベイズの定理を利用する．
今回の計算に利用する尤度関数f(E_n|θ)について，報酬は金(g)，銀(s)，銅(c)の３種であるから，3次の多項分布に従う．
事前に確率に関する情報を得られていないため，理由不十分の原則に従い，事前分布p(θ)は一様分布
に従うものとしたいところだが，これだと最尤推定法と変わらない結果になってしまう．上の定理に代入するとわかるが，事後分布と尤度関数が比例してしまう．何が問題かというと，最尤推定法では「金が全く出なかった場合」などの特殊な場合に対応できない（確率0となってしまう．この現象をゼロ頻度問題という）．そこでもっと常識的な推定をするため，事前分布と事後分布が変わらないような分布を仮定してあげることで解決する．これを自然共役分布といい，K=3の多項分布の場合には3次のDirichlet分布
となる．事後分布は尤度関数×事前分布であるから，
である．
（2015/9/7補足：読み返しているうちに「なんでガンマ関数の中に何の説明もなく多項分布のn_k!がねじ込まれてんだよおかしいだろ」と上式に疑問が生じましたが、よくよく考えてみたら重要なのは確率変数p_kの部分だけなので、あとのガンマ関数部分はDirichlet分布に無理矢理合わせるための単なる規格化定数と考えられます（多項分布のn_k部分が１段目に含まれていないのも同じ理由）。つまり上式の２段目は=ではなく∝で結ばれるべきです。）
このように自然共役分布を使うと事前分布と事後分布がパラメータのみ異なる同じ関数で表される．則ち，報酬が1個増えるごとにDirichlet分布の対応するパラメータがα_k→α_k +1とベイズ更新されることになる．たとえばE_i=Cだった場合，α_c→α_c +1という風に銅報酬に対応する分布のパラメータのみが更新される．なお，Dirichlet分布の主要な統計量は次の通り．
ちなみに，以上の話の結果だけを抽出すると，ゼロ頻度問題を解決するためのLaplaceスムージング
が，実は全ての初期値α_kが一定であるという仮定を入れたDirichlet分布を事前分布としたベイズ推定で算出した期待値と一致する，という結果が得られる．
やってみた
さて，実際に自分のログを元にExcelでベイズ推定してみる．最終結果だけを言うと，全132試行のうち，金が4回，銀が38回，銅が90回出現した．ただし自分のログでは，１回のメドフェスで２個以上貰える報酬も別々の試行としてカウントしており，友人サポートで貰える報酬も1試行としている．また当然ながら，金銀報酬UPのサポートは一切つけていない．
単純に割り算で確認すると，p_g=0.030, p_s=0.288, p_c=0.682である．まあこれでも悪い推定ではないのだが，せっかくなのでnが小さい時の途中途中の推定確率も出してみよう．
1. 事前分布パラメータの決定
パラメータが1より小さいとそのパラメータに付随する確率の最頻値がマイナス値をとってしまい，1だと最尤推定法と同様ゼロ頻度問題が生じるので，パラメータは1より大きくすべきである．また，金報酬は確率が非常に小さいこと，金銀銅の順に確率が大きくなることを考慮して，
α_g=1.05, α_g=2.35, α_g=3.65
とした．小数点以下は（両方の意味で）適当に決めた．
2. 報酬によるパラメータのベイズ更新
小難しい理論を並べたが，結局は上述したように，金が出たらα_g → α_g + 1，銀が出たらα_s → α_s + 1，銅が出たらα_c → α_c + 1とカウントしていくだけである．つまり今回のログでは，最終段階で
α_g=5.05, α_g=40.35, α_g=93.65
となる．
3. 各時点での統計量の計算
今回は事前分布，事後分布にDirichlet分布を仮定しているので，上述したDirichlet分布の統計量の公式に当てはめるだけで計算が終わる．則ち，最終段階での金報酬の平均値，最頻値，分散はそれぞれ，
である．銀，銅の場合も計算は同様．
4. グラフ化
事後確率の最頻値および平均値のベイズ推定値は以下のように収束していく．どうやら，金報酬のようにごく小さい確率では平均値が最頻値を上回る傾向にあり，金をとった瞬間に一気に上がってその後指数関数的に減少する，という動きになるようだ．当然ながら，最初の方は情報が少ないため推定値が荒れている．なお分散は，回数を重ねるたびに0に近づいていくという至極当たり前の結果が得られた．本当はルート取って標準偏差にすべきなんだけどだいたいの傾向を見るにはこれで十分．
番外編
さて先程は確率分布の統計量で推定したが，統計量だけだと全体像が見えづらいのでDirichlet分布を直接見てしまいたいという欲求に駆られる．幸いにして今回の調査対象は金・銀・銅の３種類だったので，全体の分布を直接２次元プロット可能である．「えっ３次元＋分布の４次元じゃんどうすんの？」って一瞬思ったかもしれないが，よくよく見ると確率の総和が１であるという拘束条件を使って３次元に落とし込める．今回特に知りたいのは金・銀の確率なので，これらを先頭行・列に持ってきて0~1まで自分の知りたい分解能で並べ，中のマス目をDirichlet分布の式
に従って計算する．ただしExcelで計算する場合は，オーバーフロー・アンダーフローが怖いので対数をとってから計算するのが無難．なお，計算結果は当然p_g + p_s ≧ 1の箇所とどちらかが0になる箇所は未定義となる（確率0の地点は対数を経由しなければ0として出る）．計算結果は以下のグラフの通り．等高線が高い場所ほどそれらしき場所である．
結論
金報酬が出る確率...約3%
銀報酬が出る確率...約30%
銅報酬が出る確率...約67%
が尤もらしいことが今回の結果から分かりました．次回は報酬UPに本当に差があるのかを統計的に明らかにしたいと思います．
参考文献
主に以下のサイトでベイズ推定について学びました．この場を借りて感謝します．
説明が丁寧で例題も多く，非常にわかりやすかったです．
従来の推定法とベイズ推定法の違い | Sunny side up!
ベイズ主義の考え方がよくわかりました．
歪んだサイコロでベイズ - Negative/Positive Thinking
ほぼ同じ問題だったので具体的な解き方について非常に参考になりました．
Bernoulli distribution and multinomial distribution (ベルヌーイ分布と多項分布)
Bernoulli分布と多項分布の関係がよくわかりました．
Dirichlet分布で感じた疑問に的確な答えを見せてくれました．
http://hosho.ees.hokudai.ac.jp/~kubo/ce/2013/kubo2013esj.pdf
統計的な態度がどうあるべきかが詳しく書いてあり，非常に勉強になりました．
（2015/9/7追記）
LDAでは何故ディリクレ分布を仮定するのか - SELECT * FROM life;
多項分布×Dirichlet分布の疑問が晴れました。
XAudio2をC++/CLIで使いたい 実装編
前回はXAudio2(というかSlimDX)を使う上での注意事項を散々書いてきたので，いよいよ今回は実装に移ります．C++/CLIの自分用覚書も兼ねてるので当たり前のことも書いています．
あんまり一般的な書き方じゃない気もしますが，とりあえず動くものを作ってみました．実行するとボタンが1つだけあるウインドウが現れて，ボタンを押すと押した瞬間にwave形式の波形データをメモリ上に作り，作ったデータを即再生します．
重要な部分をピックアップ
名前空間を宣言してるのでSlimDX::などは省略されていることに注意してください．
TestForm.h
using namespace System; using namespace System::ComponentModel; using namespace System::Collections; using namespace System::Windows::Forms; using namespace System::Data; using namespace System::Drawing; using namespace SlimDX::XAudio2; using namespace SlimDX::Multimedia; //----------------中略----------------------- #pragma endregion public: XAudio2 ^device; MasteringVoice ^xMVoice; WaveFormat ^format; AudioBuffer ^xBuf; IO::MemoryStream ^mStream; SourceVoice ^xSrcVoice; private: System::Void button1_Click(System::Object^ sender, System::EventArgs^ e); System::Void TestForm_Load(System::Object^ sender, System::EventArgs^ e); System::Void TestForm_FormClosing(System::Object^ sender, System::Windows::Forms::FormClosingEventArgs^ e); }; } #include <Windows.h> #define PI 3.1415926535f using namespace RealtimeRendering; using namespace std;
ヘッダファイルのSystem名前空間のところにしれっとSlimDXが混ざっています．これはXAudio2のハンドルをクラス内で宣言するのにこのあたりに書くのが都合がいいからです．また，ここでハンドルを宣言だけしておくとメンバ関数間で同じハンドルを使えます（つまり初期化と解放と処理が別々の関数で行なえる）．ハンドルはpublicにしないとコンパイルエラーになります．最後のusingは別にTestForm.cppに書いてもいいんだけどなんとなく．
TestForm.cpp
//XAudio2の初期化
device = gcnew XAudio2();
xMVoice = gcnew MasteringVoice(device);
フォームをロードしたときにXAudio2オブジェクトを生成します．SlimDXだからHRESULTのtry-catch例外処理で囲わずに作れるのが楽でいいですね．
// 波形フォーマットの生成
format = gcnew WaveFormat();
format->FormatTag = static_cast<WaveFormatTag>(WAVE_FORMAT_PCM); // フォーマット形式:PCM
format->Channels = 1; // チャネル数:モノラル
format->BitsPerSample = 16; // 16 bit/sample
format->SamplesPerSecond = 44100; // サンプリングレート
format->BlockAlignment = format->BitsPerSample / 8 * format->Channels; // 1ブロックあたりのbyte数
format->AverageBytesPerSecond = format->SamplesPerSecond * format->BlockAlignment; // 1秒あたりのバッファ数
PCMで波形フォーマットを生成します．簡単なモノラル波形を生成できればいいのでこのようにしました．ただしXAudio2とは異なり，FormatTagに定数を入れるときにSlimDXでは型が違うと怒られるのでcastしてやらなくてはいけません．
// Sin波生成(double -> short) float a = 0.5; // 振幅(0, 1] float freq = 880.0; // 周波数[Hz] float time = 0.1; // 時間[sec] array<short>^ waveData = gcnew array<short>((int)(format->AverageBytesPerSecond * time)); // t秒分のバッファ for (int i = 0; i < waveData->Length / 2; i++) { waveData[i] = (short)(SHRT_MAX * a * Math::Sin(i * PI * freq / (double)format->SamplesPerSecond)); } // short -> byte(PCMフォーマットでメモリ上に波形データを展開) array<unsigned char>^ byteArray = gcnew array<unsigned char>(waveData->Length * format->BitsPerSample / 8); // 16bit音源なので2倍 Buffer::BlockCopy(waveData, 0, byteArray, 0, byteArray->Length);
波形データはarrayで確保しておきます．最初vector
PCMフォーマットでは8bitが基本ブロックなので，16bitモノラル音源では2byteのshortを1byteのunsigned charに分割しながら入れていかなきゃいけません．そこで大活躍するのがBuffer::BlockCopy()です．これを知る前は自力でビットシフトしながらforでコピーしていました．1行で書けるって素晴らしい…
// バッファの作成 xBuf = gcnew AudioBuffer(); mStream = gcnew IO::MemoryStream(byteArray); xBuf->AudioData = mStream; xBuf->AudioBytes = (int)mStream->Length; xBuf->Flags = BufferFlags::EndOfStream;
AudioDataにどうやってデータを渡すかで苦労したのを憶えています．普通wavファイルから読み込むときはWaveStreamで何の苦も無くできたのですが，自分で波形を作るとなるとMemoryStreamにarrayを渡してそれを渡す，という格好になります．AudioBufferのメンバがどういうものなのか後でちゃんと調べないと．
// Source Voice の作成
xSrcVoice = gcnew SourceVoice(device, format);
xSrcVoice->SubmitSourceBuffer(xBuf);
多分半永久的に音を生成し続けるようなコードになると最初にSourceVoiceオブジェクトを作った方がいいと思います．勘だけど．
// 再生 xSrcVoice->Start(); MessageBox::Show("終了します．"); // 終了 xSrcVoice->Stop();
ボタンでStart(),Stop()メソッドだけを呼び出して，他の部分で波形生成・バッファ転送するような仕組みにしたいですが，今のところはこんな感じで．
delete xBuf; delete xSrcVoice; delete mStream; delete xMVoice; delete device;
C++/CLIのハンドルはいつ解放されるかわからない（解放されないかもしれない）ので，明示的に解放します．
ちなみに今回のプログラムでは，こいつを怠ると実行時に
0xC0020001: その文字列結合は無効です。
なんていう見かけないエラーが出る時があります（出ないときもあるから厄介）．
以下全ソースコード
RealtimeRendering
TestForm.h
#pragma once namespace RealtimeRendering { using namespace System; using namespace System::ComponentModel; using namespace System::Collections; using namespace System::Windows::Forms; using namespace System::Data; using namespace System::Drawing; using namespace SlimDX::XAudio2; using namespace SlimDX::Multimedia; /// <summary> /// TestForm の概要 /// </summary> public ref class TestForm : public System::Windows::Forms::Form { public: TestForm(void) { InitializeComponent(); // //TODO: ここにコンストラクター コードを追加します // } protected: /// <summary> /// 使用中のリソースをすべてクリーンアップします。 /// </summary> ~TestForm() { if (components) { delete components; } } private: System::Windows::Forms::Button^ button1; protected: private: /// <summary> /// 必要なデザイナー変数です。 /// </summary> System::ComponentModel::Container ^components; #pragma region Windows Form Designer generated code /// <summary> /// デザイナー サポートに必要なメソッドです。このメソッドの内容を /// コード エディターで変更しないでください。 /// </summary> void InitializeComponent(void) { this->button1 = (gcnew System::Windows::Forms::Button()); this->SuspendLayout(); // // button1 // this->button1->Location = System::Drawing::Point(98, 121); this->button1->Name = L"button1"; this->button1->Size = System::Drawing::Size(75, 23); this->button1->TabIndex = 0; this->button1->Text = L"button1"; this->button1->UseVisualStyleBackColor = true; this->button1->Click += gcnew System::EventHandler(this, &TestForm::button1_Click); // // TestForm // this->AutoScaleDimensions = System::Drawing::SizeF(6, 12); this->AutoScaleMode = System::Windows::Forms::AutoScaleMode::Font; this->ClientSize = System::Drawing::Size(284, 261); this->Controls->Add(this->button1); this->Name = L"TestForm"; this->Text = L"TestForm"; this->FormClosing += gcnew System::Windows::Forms::FormClosingEventHandler(this, &TestForm::TestForm_FormClosing); this->Load += gcnew System::EventHandler(this, &TestForm::TestForm_Load); this->ResumeLayout(false); } #pragma endregion public: XAudio2 ^device; MasteringVoice ^xMVoice; WaveFormat ^format; AudioBuffer ^xBuf; IO::MemoryStream ^mStream; SourceVoice ^xSrcVoice; private: System::Void button1_Click(System::Object^ sender, System::EventArgs^ e); System::Void TestForm_Load(System::Object^ sender, System::EventArgs^ e); System::Void TestForm_FormClosing(System::Object^ sender, System::Windows::Forms::FormClosingEventArgs^ e); }; } #include <Windows.h> #define PI 3.1415926535f using namespace RealtimeRendering; using namespace std;
TestForm.cpp
#include "TestForm.h" //---------------------------------------- // XAudio2の初期化処理 //---------------------------------------- System::Void TestForm::TestForm_Load(System::Object^ sender, System::EventArgs^ e) { //XAudio2の初期化 device = gcnew XAudio2(); xMVoice = gcnew MasteringVoice(device); // 波形フォーマットの生成 format = gcnew WaveFormat(); format->FormatTag = static_cast<WaveFormatTag>(WAVE_FORMAT_PCM); // フォーマット形式:PCM format->Channels = 1; // チャネル数:モノラル format->BitsPerSample = 16; // 16 bit/sample format->SamplesPerSecond = 44100; // サンプリングレート format->BlockAlignment = format->BitsPerSample / 8 * format->Channels; // 1ブロックあたりのbyte数 format->AverageBytesPerSecond = format->SamplesPerSecond * format->BlockAlignment; // 1秒あたりのバッファ数 } //---------------------------------------- // メモリ解放処理 //---------------------------------------- System::Void TestForm::TestForm_FormClosing(System::Object^ sender, System::Windows::Forms::FormClosingEventArgs^ e) { delete xBuf; delete xSrcVoice; delete mStream; delete xMVoice; delete device; } //---------------------------------------- // Sin波生成直後に再生 //---------------------------------------- System::Void TestForm::button1_Click(System::Object^ sender, System::EventArgs^ e) { // Sin波生成(double -> short) float a = 0.5; // 振幅(0, 1] float freq = 880.0; // 周波数[Hz] float time = 0.1; // 時間[sec] array<short>^ waveData = gcnew array<short>((int)(format->AverageBytesPerSecond * time)); // t秒分のバッファ for (int i = 0; i < waveData->Length / 2; i++) { waveData[i] = (short)(SHRT_MAX * a * Math::Sin(i * PI * freq / (double)format->SamplesPerSecond)); } // short -> byte(PCMフォーマットでメモリ上に波形データを展開) array<unsigned char>^ byteArray = gcnew array<unsigned char>(waveData->Length * format->BitsPerSample / 8); // 16bit音源なので2倍 Buffer::BlockCopy(waveData, 0, byteArray, 0, byteArray->Length); // バッファの作成 xBuf = gcnew AudioBuffer(); mStream = gcnew IO::MemoryStream(byteArray); xBuf->AudioData = mStream; xBuf->AudioBytes = (int)mStream->Length; xBuf->Flags = BufferFlags::EndOfStream; // Source Voice の作成 xSrcVoice = gcnew SourceVoice(device, format); xSrcVoice->SubmitSourceBuffer(xBuf); // 再生 xSrcVoice->Start(); MessageBox::Show("終了します．"); // 終了 xSrcVoice->Stop(); }
main.cpp
//------------------------------------------------------- // XAudio2によるリアルタイムサウンド合成 //------------------------------------------------------- #include "TestForm.h" using namespace RealtimeRendering; [STAThreadAttribute] int main(array<System::String ^> ^args) { Application::EnableVisualStyles(); Application::SetCompatibleTextRenderingDefault(false); Application::Run(gcnew TestForm()); return 0; }
XAudio2をC++/CLIで使いたい
初めて書くブログなので見づらくてもご容赦を…
VS2012から標準でUIプリセットが無くなり，Microsoftから見捨てられつつあるC++/CLIですが，OpenCVのようにC/C++ベースのAPIをフル機能で使いながらWindowsのGUIを作りたい俺にはC#という選択肢が無かった．
今回は簡単なリアルタイムファンクションジェネレータ（波形生成器）をGUIで作成したかったため，GUIに使い慣れたC++/CLIを，音響にXAudio2を利用することにしました．最初はWindows APIのPlaySound()で行けるかと思ったけど，細かい指定ができない上に，リアルタイム生成＆再生する都合上一旦WAVファイルに書き出すなんて真似はしたくなかったので，初の試みですがDirectXを使いました．また，DirectSoundなんか使うな！XAudio2を使え！という動きがあるみたいなのでビッグウェーブに乗ってみました．
本当はOpenALみたいなWinとMac両方で使えるクロスプラットフォームAPIを利用したかったのですが，なんか3D音源定位が主目的のライブラリっぽかったので（違ったらごめんなさい）今回はパス．機会があったら使ってみたい．
さて，C++/CLIでXAudio2を使うにあたって注意するべきポイント（つまり私が引っかかったところ）を列挙しておきます．なお，開発環境はVisualStudio 2013で，OSはWindows 8.1です．
-
MSDNの説明に潜む罠
XAudio2はDirectXの１ライブラリなのでMSDNにリファレンスがあります
（参考：XAudio2 移行ガイド ，XAudio2 ライブラリのバージョン）．そこには，
どちらのバージョンの XAudio2 も、xmcore ライブラリに依存しています。どちらかの XAudio2 ライブラリにリンクしているアプリケーションは、xmcore.lib にもリンクする必要があります。
なーんて書いてあるのですが，結論としては（CUIアプリケーションの場合でも）XAudio2を使うにあたってこいつらを特別に #pragma comment(lib, "xmcore.lib") かなんかで追加する必要はないです（XAudio2.libも不要）．私はこの"xmcore.lib"が一体どこにあるのかを探すために2日以上を浪費してしまいました．C:直下で検索しても出てこないし，新たにDLする必要があるのかとネットの海を探し回った挙句ウイルスに引っかかりそうになったり...皆さんも注意してくださいね．
-
XAudio2はC++/CLIと競合する
…タイトルを否定する事項が出てきました．
XAudio2を使うときにはどうやらプロジェクトのプロパティで「マルチスレッド(/MT)」を指定しなきゃいけないようで，こいつがC++/CLIの「共通言語ランタイムサポート(/clr)」と競合するんですね．これはどうしようもなかったので，仕方なくDirectXを直接使うのは諦めて.NETラッパーのSlimDXを使うことにしました．ここまで来ると何でわざわざC++/CLIにしたのか疑問が湧きそうになったけど，流石ラッパーといったところ，使い方がちょっと楽になってるのでまあ良しとしよう．
-
SlimDXの最新版（January 2012）にはサンプルソースがない
これには困りました．なんせどこのサイトを見ても当たり前のようにsamples内のサンプルプロジェクトに沿ってどーだかって書いてあるのに，俺のSlimDX SDKのディレクトリにはそのsamples自体が無い．漸くある一つのブログで上記の内容に触れてあったからよかったものの，それすら無かったらどうなっていたことか...
解決方法としては，過去バージョンのSeptember 2011をインストールするとちゃんとsamplesディレクトリが存在するので，その中身の"SlimDX Samples.exe"を実行・展開して事なきを得ました．
次回，Sin波のリアルタイム生成コードを載せます．
