底辺大学の院生がプログラミングや機械学習を勉強するブログ

勉強していることを雑にまとめるブログ。当然、正しさの保証は一切しない。

スマートポインタ

C/C++

qiita.com

スマートポインタについての知識がないC++使いは嫌らしいので猛省して勉強することとする。

スマートポインタ概要

ポインタを用いてメモリを動的確保をする場合、解放のし忘れが無いようにコーディングしなければならないことは周知の事実である。
ただし人間がコードを書く以上、どうしても解放のし忘れが起こりうる。
そこで、自動的にメモリの解放をさせようという仕組みがスマートポインタである。

C++11で利用できるスマートポインタは何種類かあるが、中には使うことが奨励されていないスマートポインタ（auto_ptr）もあるので、調べたところ主に使われている3種類について、使い方とその違いをまとめる。
なお、スマートポインタを利用するためにはmemoryヘッダーをインクルードする必要がある。

unique_ptr

先ほど述べた、よく使われる3種の中では最も基本的なスマートポインタがunique_ptrだ。
使い方はこんな感じ。

std::unique_ptr<int> int_ptr(new int(3)); // 宣言
std::cout << *int_ptr << std::endl;       // 間接参照演算子でアクセス

宣言のところだけが生のポインタと異なる部分で、その他使い方は同じ。
関節参照演算子*でデータが参照できるし、クラスであればアロー演算子->でメンバ関数へのアクセスや関数の呼び出しが行える。

これを見てわかる通り、スマートポインタはただのテンプレートクラスだ。
イメージとしては、unique_ptrクラスがポインタとそれに対する処理をラップしていて、コンストラクタで確保したメモリのアドレスを渡して管理を委託し、オブジェクトがスコープから外れると自動的に呼びだされたunique_ptrのデストラクタでポインタがdeleteされる感じか。仕組みがわかれば意外と単純な作り。

ただ、unique_ptrにはこれに加えて、オブジェクトをコピーできないという特徴を持つ。
つまり、こういうことができない。

std::unique_ptr<int> int_ptr1(new int(3));
std::unique_ptr<int> int_ptr2(int_ptr1);   // error
std::unique_ptr<int> int_ptr3 = int_ptr1;  // error

スマートポインタには所有権という概念があって、この例ではnew int(3)で確保したメモリの所有権はint_ptr1が持っている。
unique_ptrはその名の通り、複数のunique_ptrが1つのメモリの所有権を持つことを禁止しているユニークなスマートポインタだ。
だから、上の例ではint_ptr1が所有権を持つメモリの所有権をint_ptr2は持つことが出来ず、コンパイルエラーとなる。

ただし所有権の移動（ムーブ）はできる。

std::unique_ptr<int> int_ptr1(new int(3));
std::unique_ptr<int> int_ptr2(std::move(int_ptr1));  // ムーブ
std::unique_ptr<int> int_ptr3 = std::move(int_ptr2); // ムーブ

ムーブは所有権をムーブ先に委譲し、ムーブ元は管理していたポインタについての関連の一切を断つ。
故に、所有権は1つのunique_ptrしか持たず、所有権のユニークさは保たれる。


メモリの確保のし直しや解放はreset関数で行える。

std::unique_ptr<int> int_ptr1(new int(3));
std::cout << *int_ptr1 << std::endl; // '3'と表示
int_ptr1.reset(new int(5));          // 確保のし直し
std::cout << *int_ptr1 << std::endl; // '5'と表示
int_ptr1.reset();                    // 解放
std::cout << *int_ptr1 << std::endl; // error

メモリを確保し直した場合でも古いメモリは解放されるので安心。

shared_ptr

基本的にはunique_ptrでいいと思う。ただ、設計上複数のスマートポインタで同一のメモリを管理したい場面もあると思う。
そこで登場するのがshared_ptrだ。

基本的な使い方はunique_ptrと同じだけど、shared_ptrはunique_ptrで禁止されていたコピーができるようになってる。

std::shered_ptr<int> int_ptr1(new int(3));
std::shered_ptr<int> int_ptr2(int_ptr1);   // ok
std::shered_ptr<int> int_ptr3 = int_ptr1;  // ok

// make_shared
std::shared_ptr<int> int_ptr4 = std::make_shared<int>(3);

内部的には同一のメモリの所有権を持っているスマートポインタの数がカウントされていて、カウントが0になると自動的に解放が行われるらしい。

また、実用の際はint_ptr4のようにmake_sharedを使ったほうが高速らしい。
ただあんまり使うなって言っている人もいるしこの辺は議論の余地がありそう。

cflat-inc.hatenablog.com

weak_ptr

shared_ptrには循環参照という弱点がある。
循環参照は、AがBのshared_ptrを保持し、BもまたAのshared_ptrを保持するような状況を指す。
この場合、Aが所有権を破棄してもBがAの所有権を保持しているためAのデストラクタは呼ばれず、
またBが所有権を破棄してもAがBの所有権を保持しているため、やはりBのデストラクタは呼ばれない。
結果として、最後までAとBは解放されず、メモリリークとなる。

この問題を解決できるスマートポインタがweak_ptrだ。

weak_ptrはそれ自身はメモリの所有権を持たないが、shared_ptrの指すメモリを参照することができる。
使用の際は、参照している間に参照元のshared_ptrが解放されないようにlock関数を使ってアクセスする。

auto shared_p = std::make_shared<int>(3);
std::weak_ptr<int> weak_p = shared_p;

if (auto ptr = weak_p.lock())
{
	std::cout << *ptr << std::endl;
}

使い分けは？

スマートポインタそのものについても含めてこのページが詳しい。
qiita.com


簡単に書くと、基本的にはunique_ptrを使い、
どうしても複数のスマートポインタで管理したい場合は、使用している間は絶対に解放されてほしくないポインタについてはshared_ptrを使い、
参照できなくても問題のない場合weak_ptrを使うといった感じか。



うーむスマートポインタ……存在は知ってたけどこんなに奥が深いとは。
ただ使い方は簡単だしコードの書き換えもの負担も大きくないから、これからはちゃんと使っていこう。

おまけ

昨日の乱数を生成する関数オブジェクトをスマートポインタで書く。

// TEMPLATE CLASS UniformRand
template <class Ty_ = int>
class UniformRand
{	// wrap random number generator
	// from uniform integer distribution as function object
private:
	std::unique_ptr<std::mt19937> mt_;
	std::unique_ptr<std::uniform_int_distribution<Ty_>> ud_;

public:
	explicit UniformRand(const Ty_ min = 0,
		const Ty_ max = std::numeric_limits<Ty_>::max(),
		const long seed = static_cast<long>(time(0))) :
		ud_(new std::uniform_int_distribution<Ty_>(min, max)),
		mt_(new std::mt19937(seed))
	{	// construct from parameters
	}

	void init(const Ty_ min, const Ty_ max,
		const long seed = static_cast<Ty_>(time(0)))
	{	// initialize
		ud_.reset(new std::uniform_int_distribution<Ty_>(min, max));
		mt_.reset(new std::mt19937(seed));
	}

	Ty_	operator()(void) const
	{	// return random value
		return ud_->operator()(*mt_);
	}
};

コンストラクタ初期化子を使う

C/C++

多分知ってる人からしたら今更な内容。



今までコンストラクタでメンバ変数を初期化する際はこんなコードを書いてた。

MyClass(const int value) { value_ = value; }

これと同じことが下のコードでもできる。

MyClass(const int value) : value_(value) {}

何が違うの？という感じだし読みにくいから前者で書いてたけど、調べたら微妙に違うらしい。
どうやら、前者が宣言してから代入しているのに対して、後者は宣言と同時に初期化をしているようだ。

つまり、

int value_;
value_ = 1;

と、

int value_ = 1;

の違い。全く知らんかった……。
見ての通りでわかるけど処理時間は下のように書いた方が若干早いらしい。


もちろんメンバ変数がポインタでも初期化できる。

MyClass(const int value) : value_(new int(value)) {}

うん、最初可読性云々書いたけどこっちのほうが格好いい気がしてきた。


おまけ

一様分布に従う乱数を作る関数オブジェクトを書き換える。

#include <iostream>
#include <random>
#include <ctime>
#include <limits>

template<class Ty_ = int>
class UniformRand
{
private:
	std::mt19937* mt_;
	std::uniform_int_distribution<Ty_>* ud_;
public:
	explicit UniformRand(const Ty_ min = 0, const Ty_ max = std::numeric_limits<Ty_>::max(), const long seed = static_cast<long>(time(0))) :
		ud_(new std::uniform_int_distribution<Ty_>(min, max)),
		mt_(new std::mt19937(seed)) {}
	~UniformRand(void) { delete mt_; delete ud_; }
	Ty_ operator()(void) const { return ud_->operator()(*mt_); }
};

int main()
{
	UniformRand<> uRand(0, 10);
	for (int i = 0; i < 10; i++)
	{
		std::cout << uRand() << std::endl;
	}
	return EXIT_SUCCESS;
}

円周率を求めるシミュレータをマルチスレッド化して高速化する

C/C++

C++11を使ったマルチスレッドプログラミング。

モンテカルロ法で円周率を求める。
モンテカルロ法 - Wikipedia


一番簡単に非同期処理ができるstd::threadで実行してもいいけど、こちらは戻り値を返せない。
一応、結果を入れたい変数を参照渡しするとか、結果をグローバル変数に代入するとか色々手はあるけど、変数へのアクセスが他スレッドと競合しないよう排他制御を考慮する必要があるので非常に面倒。

そこで戻り値を返すことができる非同期実行であるstd::asyncを使用する。

std::future<int> f = std::async(std::launch::async, func);

結果はgetメソッドで取得できる。

const int result = f.get();

また、threadもasyncもSTLコンテナに入れることができるので、同じ処理を複数のスレッドで実行したい、みたいなときは非常に便利。

std::vector<std::future<int>> v;
v.push_back(std::async(std::launch::async, func));

関数部はラムダ式でもいいし、もちろんキャプチャもできる。

int a = 1;
std::vector<std::future<int>> v;
v.push_back(std::async(std::launch::async, [a]{ return a; }));

というわけでこれらを駆使してコードを書く。

#include <iostream>
#include <future>
#include <random>

#define _USE_MATH_DEFINES
#include <math.h>


int main(int argc, char** argv)
{
	// 引数が不足していればエラー
	if (argc != 3)
	{
		printf("$N_thread $N_trial\n");
		return EXIT_FAILURE;
	}

	// スレッド数
	const int N_thread = std::atoi(argv[1]);
	printf("# N_thread\t= %d\n", N_thread);

	// 試行回数
	const long N_trial = std::atol(argv[2]);
	printf("# N_trial\t= %d\n", N_trial);


	// スレッドの実行
	std::vector<std::future<long>> v;
	for (int t = 0; t < N_thread; t++)
	{
		const long N_trial_per_thread = t == 0 ? N_trial / N_thread + N_trial % N_thread : N_trial / N_thread;
		v.push_back(std::async(std::launch::async, [N_trial_per_thread]
		{
			// 乱数生成器
			std::random_device rd;
			// メルセンヌ・ツイスタ
			std::mt19937 mt(rd());
			// [0:1]の一様分布
			std::uniform_real_distribution<double> ud(0.0, 1.0);

			long cnt = 0;
			for (long t = 0; t < N_trial_per_thread; t++)
			{
				cnt += std::pow(ud(mt), 2.0) + std::pow(ud(mt), 2.0) < 1.0 ? 1 : 0;
			}
			return cnt;
		}));
	}

	// 各スレッドの結果を集計
	long cnt_sum = 0;
	for (int t = 0; t < N_thread; t++)
	{
		cnt_sum += v[t].get();
	}


	// 結果を印字
	printf("# \n");
	printf("# M_PI\t\tMonteCarlo\n");
	printf("%.10f\t%.10f\n", M_PI, 4.0 * cnt_sum / N_trial);
	
	return EXIT_SUCCESS;
}

スレッド数4、試行回数100,000,000回で実行した。

f:id:telltales:20160702050050p:plain

小数第3位までしか合わないのかー、なかなか精度が低い。


ちなみにうちのノート君は論理コアが4つしか無いので4スレッドで実行するとCPU使用率が100％になる。

f:id:telltales:20160702050045j:plain

他のことをしたいのであれば、コア数-1ぐらいにスレッド数を設定するといい感じだった。



参考にしたサイト。

qiita.com

非同期処理について分かりやすくまとめられている。

1.2.3　ベイズ確率（＋BPSKのビット判定）

PRML

最尤推定の説明がこのページだけだと分かりづらいから、実際に推定問題を解いてみて理解してみる。
個人的には信号のビット判定問題が直感的に分かりやすいと思う。（デジタル通信についての知識が必要だけど）

問題

f:id:telltales:20160701005528p:plain

BPSKの信号点は、実軸上に2つ現れる形式が一般的。上の例だと、+1と-1のところ。

位相偏移変調 - Wikipedia

このどちらかを送信すると、受信機は伝搬空間中の雑音の影響で、若干元の信号とは違った点にある信号が受信される。

f:id:telltales:20160701005532p:plain

受信機は元の信号から若干ずれて届いた受信信号を見て、この信号はもともとは+1と-1どちらの信号だったのかを判定する必要がある。
上の例だと近いのは+1だけど、もしかしたら強い雑音のせいで-1の方から飛んできた信号なのかもしれない。

さてどのような判定が最適か？



大体の人がこう考える。

f:id:telltales:20160701004737p:plain

虚軸を境界として判定する。
つまり虚部を無視して、受信信号の実部が正であれば+1、負であれば-1と判定する。要は近い方と判定しましょうという考え方。

結果的に見ると確かにこの考え方は正しい。でも、なんでそういうことが言えるの？となると、なんとなく……としか答えられないことが多い。
今日はこの判定問題についてちゃんと考えてみる。

最尤推定

BPSKは実部だけ考えれば良いので、信号はすべて実数で考える。また雑音は加法性のガウス雑音を仮定する。

受信信号 yを次の式でモデル化。

y=x+n

x：+1か-1のどちらかの値を取る送信信号
n：平均0分散\sigma^2のガウス分布に従う確率変数


最適な判定規範はyを観測した時のxの事後確率を最大にするxを採用することだから、ベイズの定理を用いると

\newcommand{\argmax}{\mathop{\rm arg~max}\limits} \begin{eqnarray} \hat{x} &=& \argmax_{x\in\{+1,-1\}}p(x|y) \\ &=& \argmax_{x\in\{+1,-1\}}\frac{p(y|x)p(x)}{p(y)} \\ &=& \argmax_{x\in\{+1,-1\}}p(y|x)p(x) \end{eqnarray}

ここでもし事前分布 p(x)に関する情報、すなわち xの生起確率についての知識がなければ、それぞれ50%の確率で現れると仮定する。
すると p(x)は定数となるので xの判定に関与せず、最終的に

\newcommand{\argmax}{\mathop{\rm arg~max}\limits} \hat{x} = \argmax_{x\in\{+1,-1\}}p(y|x)

となる。
p(x|y)を最大にする xを推定値とする手法は最大事後確率推定（MAP推定）である。
一方で、今現れた p(y|x)を最大にする方は最尤推定（ML推定）と呼ばれている。

よって、最大事後確率推定と最尤推定の違いは、事前分布 p(x)が推定に関与するかしないかのみである。


さて、じゃあ p(y|x)はどのように書けるのか。
今 xが与えられているので定数と考える。また、nは確率変数なので、xとnの和であるyもまた確率変数となる。ここでnは平均0分散\sigma^2のガウス分布で与えられているので、p(y|x)は平均x分散\sigma^2のガウス分布になる。

よって、
\newcommand{\argmax}{\mathop{\rm arg~max}\limits} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \begin{eqnarray} \hat{x} &=& \argmax_{x\in\{+1,-1\}}p(y|x) \\ &=& \argmax_{x\in\{+1,-1\}}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y-x)^2}{2\sigma^2}\right) \\ &=& \argmin_{x\in\{+1,-1\}}(y-x)^2 \\ &=& \argmin_{x\in\{+1,-1\}}\{y^2-2xy+x^2\} \\ \end{eqnarray}

最後の式の第1項はxに依らないため無視。第3項はxがどちらの値をとっても+1になるのでやはり無視。
というわけで第2項が小さくなる方のxを採用すればよいことがわかる。

故に、最終的にな判定基準として

x= \begin{cases} +1 & (y\geq0) \\ -1 & (y<0) \end{cases}

を得ることができ、これは要するにさっき絵で示した「受信信号の実部が正であれば+1、虚部が負であれば-1と判定する」を表している。



先生は「直感的な判断を信用せず、一度はちゃんと理論的に証明したほうが良い」と何度もおっしゃってる。
実際その通りだと思うし、その計算過程は意外と楽しかったりするから、なんとなくで済ませずちゃんと確認するよう心がけよう。

「int *a」 vs. 「int* a」

C/C++

f:id:telltales:20160630210954j:plain

どっちでもいいってことで決着はついてるけど暇だし蒸し返してみる。
ちなみに僕は「int* a」派。
教科書は「int *a」が多い気がする。僕がC++を勉強した独習C++も「int *a」だった。

独習C++ 第4版

独習C++ 第4版


とりあえず逆の、「int *a」派の意見から見てみる。

「int *a」派の主張

よく見るやつがこれ。

int* foo, bar;

こう書くと、fooはポインタ変数となるが、barは普通のint型変数になる。

本来、

型名　変数1, 変数2;

という形式になっていれば変数1と変数2の両方が同じ型になるはずである。
だから「int*」は型足り得ない。また、

int *foo, *bar;

こうすれば両方ポインタ変数になることから、アスタリスクは変数側につけるべきである、というもの。

あと、「int *a」の読み方は、「変数名の左に*演算子をつけるとint型になる変数a」だというのもあった。

「int* a」派の主張

上の例に対して、じゃあ関数のプロトタイプ宣言で型がポインタ型だったらなんて書くんだよ？という疑問がある。

int* func(int*);

この場合引数の方は変数がないから（実際は変数まで書いてもいいんだけど）、変数の左側に*をつけるというルールが適用できなくなる。

戻り値についても、関数がポインタになっていることを表しているわけではなく、戻り値がint型のポインタ変数であることを示したいのだからやっぱり型側に*をつける方が自然な気がする。

ただこれについてはソース見ると

int *func(int *);

って書いてあるのも見る。
引数の方は「int *a」って書いてあってaが省略されている、と読めなくもないけど、戻り値の方は果たしてどうなんだこれは……。


もう一つ。
「そもそも宣言する変数がポインタ型であることを示す*と、ポインタの先を指す間接参照演算子としての*は別物だ」という主張。

これが結構納得行く説明だと思ってる。
だって宣言の「int *a」の*はaの先を示しているわけではないから。

実際、

int *a = 1;

というコードは書けない。（実体が無いから当たり前だけど）

結局どっち？

f:id:telltales:20160630221422j:plain

書き方が一貫してさえいればそれでいいという結論。あとコーディング規約が決まっているならそれに従いましょう。
あとそもそも生のポインタなんて使ってないで、出来るだけスマートポインタを使うべきか。

うーん、書いてて思い出したけどスマートポインタの使い方についてイマイチちゃんと勉強してないから、今度はそれを勉強しようかな。


ちなみにCに慣れている人は「int *a」、C++に慣れている人は「int* a」で書くことが多いらしい。

www.naturalsoftware.jp


一応比較のためにソースコード検索しようと思ったんだけど、searchcodeもGitHubのコード検索も*とかの記号に対応されてないらしく検索出来んかった……。

GitHubの方は記号は使えんぞってはっきり書いてあった（Searching code - User Documentation）けど、探せばアスタリスク付けて検索できるサービスあるのかな。

続き：多項式曲線フィッティング

PRML

昨日の続き。


二乗和誤差関数を最小にするような係数ベクトル {\bf w}を計算することで、様々な非線型関数へのフィッティングが出来るらしい。
そのためには例の E({\bf w})を w_mで偏微分して0とおいた方程式を M+1個作って、その連立方程式を解けばいいらしい。

頑張ればやれるだろうけど、実はその結果が次のようになることが、演習1.1に載っている。

\sum_{j=0}^MA_{ij}w_j=T_i

ただし、

A_{ij}=\sum_{n=1}^N(x_n)^{i+j} \\ T_i=\sum_{n=1}^N(x_n)^it_n

演習の内容はこうなることを示せというもの。
演習そのものは省略するとして、せっかくだから本当にこれでフィッティング出来るのか、教科書の図1.2と同じシチュエーションで試してみる。

以下、（糞）プログラムコード

#define _USE_MATH_DEFINES
#include <math.h>

#include <iostream>
#include <random>
#include <fstream>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>

typedef boost::numeric::ublas::vector<double> dvector;
typedef boost::numeric::ublas::matrix<double> dmatrix;



int main(int argc, char** argv)
{
	// 入力が足りなければエラーメッセージを出力
	if (argc != 3)
	{
		printf("$N $M\n");
		return EXIT_FAILURE;
	}

	const int N = std::atoi(argv[1]);	// データ数
	const int M = std::atoi(argv[2]);	// 次元数


	// データを生成する関数（ sin(2πx) ）
	auto f = [](double x) { return std::sin(2.0 * M_PI * x); };

	std::mt19937 mt(100);
        // 平均0, 標準偏差0.1のガウス分布に従う乱数生成器
	std::normal_distribution<> noise(0, 0.1);


	dvector x(N); dvector t(N);
	for (int n = 0; n < N; n++)
	{
		x(n) = (double)n / N;
		t(n) = f(x(n)) + noise(mt);	// f(x) + 雑音
	}


	// 係数行列
	dmatrix A(M + 1, M + 1);
	for (int n1 = 0; n1 <= M; n1++)
	{
		for (int n2 = 0; n2 <= M; n2++)
		{
			A(n1, n2) = 0;
			for (int n3 = 0; n3 < N; n3++)
			{
				A(n1, n2) += std::pow(x(n3), n1 + n2);
			}
		}
	}

	// 定数ベクトル
	dvector b(M + 1);
	for (int n1 = 0; n1 <= M; n1++)
	{
		b(n1) = 0;
		for (int n2 = 0; n2 < N; n2++)
		{
			b(n1) += std::pow(x(n2), n1) * t(n2);
		}
	}

	boost::numeric::ublas::permutation_matrix<> pm(A.size1());
	boost::numeric::ublas::lu_factorize(A, pm);		// LU分解
	boost::numeric::ublas::lu_substitute(A, pm, b);		// 前進消去と後退代入


	// csvに書き込み
	std::ofstream ofs("result.csv");
	ofs << "x[n],t[n]" << std::endl;
	for (int n = 0; n < N; n++)
	{
		ofs << std::to_string(x[n]) + "," + std::to_string(t[n]) << std::endl;
	}
	ofs << std::endl << "x,y[x;w]" << std::endl;
	for (int n = 0; n < 100; n++)
	{
		const double x_ = (double)n / 100;
		double y = 0;
		for (int m = 0; m <= M; m++)
		{
			y += b[m] * std::pow(x_, m);
		}
		ofs << std::to_string(x_) + "," + std::to_string(y) << std::endl;
	}

	return EXIT_SUCCESS;
}


うーむ、そのうち暇なときにちゃんとしたコードに直したいなこれは……。

ともあれ結果を見てみる。まずは M=1の場合。

f:id:telltales:20160629224904j:plain

全然推定できてないけど、とりあえず1次関数になってくれたし、それっぽいところ通ってるからプログラムは合ってそう。

次に M=3の場合。

f:id:telltales:20160629224902j:plain

おおおおおおおおおおおお。すごい。かなり綺麗にフィッティングしてる。
データ数10って少なくね？って感じてたけど、以外にもちゃんと学習できるんだなー。

で、最後に M=9の場合。

f:id:telltales:20160629224901j:plain

こ…これ…これは…………過学習だあああああ┗(^o^)┛wwwwww┏(^o^)┓ﾄﾞｺﾄﾞｺﾄﾞｺﾄﾞｺwwwwwww

なんか教科書ほど思いっきり外れてはいないんだけど、それでも全部の観測値を取りつつ、元の関数から大きく外れているという特徴はちゃんと見て取れる。
実際に自分でこういうことを試したことがなかったから結構感動的だ。


さて、この次元数が大きくなると発生する過学習というものは、データ数を多くすることである程度回避出来るらしい。
なんで？って感じだけどその辺詳しく書かれていない（そのうち出てくるのかな？）のでとりあえずそのまま納得するとして、同じ関数で実験してみる。

f:id:telltales:20160629230203j:plain

N=100、 M=9としてフィッティングしてみたところ、確かにデータ数を増やすことで過学習が抑えられていることが見て取れる。
というかむしろ今までで一番良い結果が得られているようにも見える。（所詮1回の試行だからなんとも言えないけど）

うーんやっぱりなんでこうなるのか気になる。
なんかデータ数とモデルのパラメータ数の関係について、「経験則としては」なんて書いてあるし、はっきりとして尺度があるのかな。3章とやらに期待である。


といったところでおしまい。帰ってセキュスペの勉強でもするか。

1.1　例：多項式曲線フィッティング

PRML

１行目「まず最初に単純な回帰問題から始めよう．」

 

……。

回帰問題ってなんやねん！！！！

 

回帰分析

回帰（かいき、英: regression）とは、統計学において、Y が連続値の時にデータに Y = f(X) というモデル(「定量的な関係の構造[1]」)を当てはめる事。別の言い方では、連続尺度の従属変数（目的変数）Y と独立変数（説明変数）X の間にモデルを当てはめること。X が1次元ならば単回帰、X が2次元以上ならば重回帰と言う。Y が離散の場合は分類と言う。

回帰分析（かいきぶんせき、英: regression analysis）とは、回帰により分析する事。

回帰で使われる、最も基本的なモデルは Y = AX + B という形式の線形回帰である。

wikipedia: 回帰分析

 

つまり、入力値 xと、ある関数 f(x)にしたがって生成されたデータ tの組み合わせをたくさん観測したとき、それぞれのベクトル {\bf x}=(x_1,\dots,x_N)^T、 {\bf t}=(t_1,\dots,t_N)^Tを使って f(x)の方を推定しようということらしい。

んで、 f(x)さえ分かってしまえばわざわざ観測しなくても、新しい \hat{x}に対する観測 \hat{t}が予測できるよね、という話か。

ただしここでは観測値にはガウス分布に従うランダムノイズが加算されていることに注意する。こういう雑音がガウス分布でモデル化されることはよくある話。無線通信とかではノイズは複素ガウスの標本としてモデル化されていた。一般に、ホワイトノイズとか呼ばれているあれ。

だから必ずしも t_n=f(x_n)とはならず、ノイズを \epsilonとすると観測値は t_n=f(x_n)+ \epsilonで与えられている。

 

さて、ここでやっと出てくる題目の「曲線フィッティング」とは、 f(x)を

y(x;{\bf w})=w_0+w_1x+w_2x^2+\cdots+w_Mx^M=\sum_{j=0}^Mw_jx^j

という M次多項式で当てはめようというアプローチのこと。

実際の関数が \sinでも \expでも、こういう形式なら線形非線型に限らずそれなりに当てはめられそう。

 

 だからやることは、観測された {\bf t}を使ってなんとかして多項式の係数ベクトル {\bf w}を求めること。そのための方法として、誤差関数を最小化しようという考え方がある。

つまり、

E({\bf w})=\frac{1}{2}\sum_{n=1}^N\{y(x_n,{\bf w}) - t_n\}^2

これが最小になる {\bf w}を採用しようというもの。

 

まあ、予測値と実際に得られた観測値の誤差が一番小さくなるようなモデルを採用しようという考え方は直感的には理解できる。ただ、本当にそれが最適な推定規範なの？という疑問を持つことは、ついつい疎かにしがちだけど、とても大切。

後半に、確かにこの規範は最適であることが述べられるが、直感的に正しいと思われる解釈に流されて思考を止めることが無いよう注意したいところ。

 

さてここまでくれば話は簡単。

E({\bf w})は下に凸の二次関数だから、最小になる {\bf w}を見つけたければ {\bf w}についての1階微分が0になる点を求めればいい。

 

 

というところで眠くなってきた。明日はこれの続きやる。

実際にプログラムを書いてちゃんと関数が推定できていることを確認するぐらいまでやりたいなー。