２０１９年やりたいこと
「年初の誓い(ﾟДﾟ)ﾉ」はしなかったんですが，「今年こんなことができたらいいな～」ってことを考えてみました．なんだかおっさんの割には好奇心だけは人一倍強くて，無数にやりたいこと・知りたいことが出てくるんですが，どうも容量が悪いのか頭の回転数が足りてないのか（笑）なかなか思ったように進まず．「あれもこれも」的な感じになっていつも思ったように進みません．
が，今年もあえて絞らずに欲望のままにいろいろとやってみたいと思います．
主題は変わらず「ロボット」「画像処理」「ソフトウェア」の３つです．
（そもそもこの切り口が広すぎるんですが．．．）
１．ロボット
SLAMの勉強
・Particle Filterを用いた SLAM (FastSLAM, GMapping) の理解・コードいじり．
・Scan Matching/最適化を用いた SLAM (GraphSLAM, CartoGrapher) の理解・コードいじり．
カメラを使った障害物検知・回避
・SGBM / Stixel を用いた障害物検知と回避の実装．
・Daimler 6D Visionの実装とロボットへの適用．
つくばチャレンジ２０１９への参加
・めざせ今年こそは完走！今年は Lidar で自己位置推定，カメラで障害物検知・回避するところまで持っていきます！
２．画像処理
３次元復元プロジェクト再開
・最適化・バンドル調整の理解・コードいじり．
・画像特徴量（Fast, SIFT etc...）の理解・コードいじり．
・OpenSFM/OpenMVGの理解・コードいじり．
・Visual SLAM の理解・コードいじり．
画像のみを用いた自己位置推定の調査
・2020年のつくばチャレンジで，カメラを使って自己位置推定を達成したいと思っているんですが，このための下準備です．
ディープラーニング
・おそらくこれは2020年に持ち越しですね．．．
３．ソフトウェア
ソフトウェアに関しては，今年は「C++」と「Python」の年にします！仕事でドンピシャのことをやっているわけではないので，基本的に本を読んで知識を蓄えていくしかないんですが，サンプルコードを書いたり，オープンソースのプロジェクトを見て・いじってモノにしていく方向でスキルを上げていけたらと思っています．
C++を（ある程度）マスター！
Pythonを（ある程度）マスター！
ソフトウェアアーキテクチャの勉強
ちなみに，自分は本をたくさん買って文字通り積読する人なんですが（笑），本棚がいっぱいになってきたこともあって「これではいかん！」と最近思いまして，，，課題図書を毎月決めて，読んでいこうと思います．目指せ一か月に２冊読破！
ということで，結局年初の誓いになってしまいましたが，，，，
１．寄り道しながら，やりたいことをやる！（ロボット・画像処理・ソフトウェアに関して）
２．月単位で目標・やることを決めてブログで宣言する．
３．月に２冊の課題図書の設定と読破！
でやっていきます！これだけできたら，年末には相当レベルアップしているはず！！実り多い２０１９年にします！！
帰省からの，，，２０１９年あけましておめでとうございます！
あけましておめでとうございます！
年末もっとブログ（雑記）をアップするつもりだったんですが，実家に帰っていたこともあり更新できずじまいでした．．．
帰省
割と時間があったのでダイエットついでに結構歩き回ったんですが，大きな橋ができたりイオンモールができたりと，風景が変わっているのを見て，「年取ったなー」なんてちょっとセンチメンタルな感じになってました（笑）
おどろいたことに，我らが「徳島阿波踊り空港」が国際空港になってました．週に２回程度ですが，徳島 <-> 香港間の定期便をキャセイパシフィック航空子会社のキャセイドラゴン航空が運営しているみたいです．
２０１９年！
ということで，年初の誓いを立てる(ﾟДﾟ)ﾉと思ったんですが，今年はちょっと岐路の年になりそうで，まだ黙ってます（笑）が，
All life is an experiment. The more experiments you make the better.
Ralph Waldo Emerson
をモットーに，今年もビビりながらいろんな実験をしていけたらと思います！！！ではでは，今年もよろしくお願いします．
EC++ #50, #51, #52 プレースメント new/delete
Effective C++ の 50, 51, 52項の内容です．本的には第8章をまるごと new/delete の話にさいていましたが，これを機会にちょっと new/delete の定義を見てみることにしました．
new/delete とは？
下記，自分の環境にあった new/delete の宣言です．
void* operator new(std::size_t) _GLIBCXX_THROW (std::bad_alloc) __attribute__((__externally_visible__)); void* operator new[](std::size_t) _GLIBCXX_THROW (std::bad_alloc) __attribute__((__externally_visible__)); void operator delete(void*) _GLIBCXX_USE_NOEXCEPT __attribute__((__externally_visible__)); void operator delete[](void*) _GLIBCXX_USE_NOEXCEPT __attribute__((__externally_visible__)); void* operator new(std::size_t, const std::nothrow_t&) _GLIBCXX_USE_NOEXCEPT __attribute__((__externally_visible__)); void* operator new[](std::size_t, const std::nothrow_t&) _GLIBCXX_USE_NOEXCEPT __attribute__((__externally_visible__)); void operator delete(void*, const std::nothrow_t&) _GLIBCXX_USE_NOEXCEPT __attribute__((__externally_visible__)); void operator delete[](void*, const std::nothrow_t&) _GLIBCXX_USE_NOEXCEPT __attribute__((__externally_visible__)); // Default placement versions of operator new. inline void* operator new(std::size_t, void* __p) _GLIBCXX_USE_NOEXCEPT { return __p; } inline void* operator new[](std::size_t, void* __p) _GLIBCXX_USE_NOEXCEPT { return __p; } // Default placement versions of operator delete. inline void operator delete (void*, void*) _GLIBCXX_USE_NOEXCEPT { } inline void operator delete[](void*, void*) _GLIBCXX_USE_NOEXCEPT { }
で，C言語の malloc/free を使っていた時は malloc(sizeof(Hoge)) って感じで構造体のサイズを指定してたので，嫌でもサイズを意識していたと思うんですが，new/delete 場合も結局サイズを指定するんですね．なので，new Hoge で実際に処理されている内容としては，，，
１．new演算子の探索を行う．
３．newが返したポインタの指す位置をthisポインタとしてインスタンスを生成する。
です．
new/delete のカスタマイズ
ということで，プレースメント new/delete の話ですが，演算子を定義しなおせばカスタマイズできます．
下記サンプルコードですが，Object というクラスに対して new 演算子を定義しました．std::ostream を引数としてとっているので，new が呼ばれるとコンソール出力されることがわかります．で，プレースメント new を定義したせいで，グローバルスコープのデフォルト new 演算子は探索対象から外れ，そのままでは呼び出すことができなくなってます．
#include <iostream> class Object { public: static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); }; void* Object::operator new(std::size_t size, std::ostream &stream) noexcept(false) { stream << "operator new" << std::endl; return ::operator new(size); } int main(int, char**) { std::cout << "Hello World!" << std::endl; // デフォルトの new operator は隠蔽される． // Object* obj = new Object; // ostream を引数として取る new operator Object* obj = new(std::cout) Object; return 0; }
で，デフォルトスコープの new 演算子を隠蔽するためにはどうすればいいか．Effective CPPでは，
１．クラスの中で演算子を定義して，デフォルトを呼んであげる．
２．基底クラスにまとめて定義して，派生クラスで using を使って可視にする．
という２つのアプローチでデフォルト new/delete 演算子を定義してました．なお，デフォルトで定義されている new 演算子は次の３つがあるので，すべて使えるようにしてあげないといけません．
static void* operator new(std::size_t size) noexcept(false); static void* operator new(std::size_t size, void* pnt) noexcept(false); static void* operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true);
#include <iostream> class Object { public: static void* operator new(std::size_t size) noexcept(false); static void* operator new(std::size_t size, void* pnt) noexcept(false); static void* operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true); static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); }; void* Object::operator new(std::size_t size) noexcept(false) { return ::operator new(size); } void* Object::operator new(std::size_t size, void* pnt) noexcept(false) { return ::operator new(size, pnt); } void* Object::operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true) { return ::operator new(size, nothrow); } void* Object::operator new(std::size_t size, std::ostream &stream) noexcept(false) { stream << "operator new" << std::endl; return ::operator new(size); } class ObjectChild : public Object{ public: // 基底クラスの new 関数も探索対象に入れる． using Object::operator new; static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); }; void* ObjectChild::operator new(std::size_t size, std::ostream &stream) noexcept(false) { return ::operator new(size); } int main(int, char**) { std::cout << "Hello World!" << std::endl; // デフォルトの new operator も定義した． Object* obj1 = new Object; delete obj1; // ostream を引数として取る new operator Object* obj2 = new(std::cout) Object; delete obj2; // 派生クラスの new operator も定義した． ObjectChild* objc1 = new ObjectChild; delete objc1; // 派生クラスで独自定義した new 演算子． ObjectChild* objc2 = new(std::cout) ObjectChild; delete objc2; return 0; }
プレースメント new に対応するプレースメント delete の定義
で，プレースメント new を定義した時に，もうひとつ気をつけることがありました．プレースメント new(std::ostream &stream) を定義した場合に．同じシグニチャを持つ delete を定義してあげないと行けないのでした．今回の例では，下記のペアの new/delete を作ってあげる必要があります．
static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); static void operator delete(void* pnt, std::ostream &stream) noexcept(false);
で，Hoge のコンストラクタが例外を発生させた場合に，new で既に獲得済のリソースを解放する必要がありますが，この時に引数が対応するバージョンの delete がないとメモリがリークしてしまうのでした．該当するサンプルコードを作ってみましたが，確かにメモリリークしてますね．
#include <iostream> class Object { public: static void* operator new(std::size_t size) noexcept(false); static void* operator new(std::size_t size, void* pnt) noexcept(false); static void* operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true); static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); }; void* Object::operator new(std::size_t size) noexcept(false) { return ::operator new(size); } void* Object::operator new(std::size_t size, void* pnt) noexcept(false) { return ::operator new(size, pnt); } void* Object::operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true) { return ::operator new(size, nothrow); } void* Object::operator new(std::size_t size, std::ostream &stream) noexcept(false) { stream << "operator new" << std::endl; return ::operator new(size); } class ObjectChild : public Object{ public: // 基底クラスの new 関数も探索対象に入れる． ObjectChild() { throw std::exception(); } using Object::operator new; static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); }; void* ObjectChild::operator new(std::size_t size, std::ostream &stream) noexcept(false) { return ::operator new(size); } int main(int, char**) { std::cout << "Hello World!" << std::endl; ObjectChild *obj2; try { // ostream を引数として取る new operator obj2 = new(std::cout) ObjectChild; delete obj2; } catch (std::exception e) { std::cout << e.what() << std::endl; } return 0; }
上記のコードを valgrind で実行すると，確かにメモリリークしていることがわかります．ここで，対応する引数の delete を作ってみましょう．サンプルコードは下記ですが，Exceptionが送出された時点で対応する delete がよばれてメモリリークが避けられていることがわかります．
#include <iostream> class Object { public: static void* operator new(std::size_t size) noexcept(false); static void* operator new(std::size_t size, void* pnt) noexcept(false); static void* operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true); static void operator delete(void* pnt) noexcept(false); static void operator delete(void* pnt1, void* pnt2); static void operator delete(void* pnt, const std::nothrow_t ¬hrow) noexcept(true); static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); static void operator delete(void* pnt, std::ostream &stream) noexcept(false); }; void* Object::operator new(std::size_t size) noexcept(false) { return ::operator new(size); } void* Object::operator new(std::size_t size, void* pnt) noexcept(false) { return ::operator new(size, pnt); } void* Object::operator new(std::size_t size, const std::nothrow_t ¬hrow) noexcept(true) { return ::operator new(size, nothrow); } void* Object::operator new(std::size_t size, std::ostream &stream) noexcept(false) { stream << "placement new" << std::endl; return ::operator new(size); } void Object::operator delete(void* pnt) noexcept(false) { ::operator delete(pnt); } void Object::operator delete(void* pnt1, void* pnt2) noexcept(false) { ::operator delete(pnt1, pnt2); } void Object::operator delete(void* pnt1, const std::nothrow_t ¬hrow) noexcept(true) { ::operator delete(pnt1, nothrow); } void Object::operator delete(void* pnt, std::ostream &stream) noexcept(false) { stream << "placement delete" << std::endl; ::operator delete(pnt); } class ObjectChild : public Object{ public: // 基底クラスの new 関数も探索対象に入れる． ObjectChild() { throw std::exception(); } using Object::operator new; using Object::operator delete; static void* operator new(std::size_t size, std::ostream &stream) noexcept(false); }; void* ObjectChild::operator new(std::size_t size, std::ostream &stream) noexcept(false) { return ::operator new(size); } int main(int, char**) { std::cout << "Hello World!" << std::endl; ObjectChild *obj2; try { // ostream を引数として取る new operator obj2 = new(std::cout) ObjectChild; delete obj2; } catch (std::exception e) { std::cout << e.what() << std::endl; } return 0; }
#10 代入演算子は *this の参照を戻す
ということで，Effective C++の１０項の内容です．
コンストラクタ，コピーコンストラクタの暗黙呼び出し
ということで，ちょっと内容とずれてしまうのですが寄り道です．そもそも，コンストラクタとかコピーコンストラクタってコンパイラが暗黙の呼び出しをしてくれます．この辺をちゃんと理解すべく，簡単にまとめておきます．サンプルケースとして，下記のクラスを作りました．
#include <iostream> class IntNumber { public: IntNumber(int val) : num(val) { std::cout << "IntNumber(int val)" << std::endl; } ~IntNumber() { std::cout << "~IntNumber" << std::endl; } IntNumber(const IntNumber& obj) : num(obj.num) { std::cout << "IntNumber(const IntNumber& obj)" << std::endl; } IntNumber& operator=(const IntNumber& val) { std::cout << "IntNumber& operator=(const IntNumber& val)" << std::endl; if (this != &val) { this->num = val.num; } return *this; } private: int num; }; int main(int, char**) { // １．直接的なコンストラクタのコール std::cout << "1. Explicit Constructor Call" << std::endl; IntNumber num1(1); std::cout << std::endl; // 2．暗黙的なコンストラクタのコール std::cout << "2. Implicit Constructor Call" << std::endl; IntNumber num2 = 2; std::cout << std::endl; // ３．明示的なコピーコンストラクタのコール std::cout << "3. Explicit Copy Constructor Call" << std::endl; IntNumber num3(num1); std::cout << std::endl; // ４．暗黙的なコピーコンストラクタのコール std::cout << "4. Implicit Copy Constructor Call" << std::endl; IntNumber num4 = num1; std::cout << std::endl; // ５．明示的なコピー代入演算子のコール std::cout << "5. Explicit operator= Call" << std::endl; num1.operator=(num2.operator=(num3.operator=(num4))); std::cout << std::endl; // ６．コピー代入演算子のコール std::cout << "6. Implicit operator= Call" << std::endl; num1 = num2 = num3 = num4; std::cout << std::endl; return 0; }
コンストラクタコール
上記サンプルコードですが，main関数の下記の部分に関しては定義と全く同じ形でコンストラクタを読んでいるので，想像どおりコンストラクタが呼ばれます．
// １．直接的なコンストラクタのコール std::cout << "1. Explicit Constructor Call" << std::endl; IntNumber num1(1); std::cout << std::endl;
では，次のケースはどうでしょうか？
// 2．暗黙的なコンストラクタのコール std::cout << "2. Implicit Constructor Call" << std::endl; IntNumber num2 = 2; std::cout << std::endl;
このケースでは，”コピー代入演算子が呼ばれる？？”ってことはなく，コンストラクタが呼ばれます．おそらく，コンパイラのなかで下記のような順番でフィットする関数を探しているのだと思います．
１．num2は今までに作成されていないオブジェクト．なので，オブジェクトを新規作成する必要あり．
２．= の右辺が int だけど，int から直接的にオブジェクトを作れるのは，，，，"IntNumber(int val)"．ということで，コンストラクタをコール．
コピーコンストラクタコール
続いて同じくコピーコンストラクタのケースを考えます．下記のケースでは，定義と同じ形でコピーコンストラクタを読んでいるので，コピーコンストラクタが呼ばれます．
// ３．明示的なコピーコンストラクタのコール std::cout << "3. Explicit Copy Constructor Call" << std::endl; IntNumber num3(num1); std::cout << std::endl;
では，次のケースでは？こちらも暗黙的なコンストラクタコールのときと同様に，，，
１．num4は今までに作成されていないオブジェクト．なので，オブジェクトを新規作成する必要あり．
２．= の右辺が IntNumber だけど，IntNumber から直接的にオブジェクトを作れるのは，，，，"IntNumber(const IntNumber& obj)"，ということで，コピーコンストラクタをコール．
// ４．暗黙的なコピーコンストラクタのコール std::cout << "4. Implicit Copy Constructor Call" << std::endl; IntNumber num4 = num1; std::cout << std::endl;
コピー代入演算子コール
続いてコピー代入演算子ですが， "operatorxx" って普通の関数と見た目が違っててちょっとわかりにくいです．ただ "operatorxx" って明示的に関数を呼ぶと，わりとすっきりわかりやすくなります．サンプルコードの例では，，，，
// ５．明示的なコピー代入演算子のコール std::cout << "5. Explicit operator= Call" << std::endl; num1.operator=(num2.operator=(num3.operator=(num4))); std::cout << std::endl; // ６．コピー代入演算子のコール std::cout << "6. Implicit operator= Call" << std::endl; num1 = num2 = num3 = num4; std::cout << std::endl;
上記の５と６は等価です．やっとEffective C++10項の内容になりますが，”代入演算子は基本的に (a = b = c) っていう書き方ができないといけない！”ってことを言ってます．
"a = b = c"
って書くと，組み込み型では c を b に代入して，（その結果が反映された） b を a に代入する．っていう意味になると思いますがこれって結局明示的に代入演算子を呼ぶと
a.operator=(b.operator=(c))
ってことなんですよね．だから，operator= が参照を返す必要があるということになります．
#6 コンパイラが自動生成する関数を使用禁止する方法
ということで，Effective C++６項のお話です．
コンパイラが自動生成する関数を禁止する方法
ということで，「自作しなければコンパイラが自動で生成してしまう．」というおはなしだったので，自分で作ってコンパイラを黙らせます．で，自作関数をpriavateスコープにしておけば，それらの関数がユーザから使えなくなります．例えば，コピーコンストラクタ・コピー代入演算子の使用を禁止したいとします．この時は下記のように「privateスコープに関数宣言をして，定義を書かない」という対応を取ればOKです．
下記の例では，コピーコンストラクタ・コピー代入演算子は「無事？」コンパイルエラーになります．
#include <iostream> #include <string> #include <vector> class UnCopyableClass { public: // コンストラクタは必要なので，宣言＆定義を書く． UnCopyableClass() {}; // デストラクタは必要なので，宣言＆定義を書く． ~UnCopyableClass() {}; private: // コピーコンストラクタは使わせないようにしたいので，宣言だけにして，praivteスコープとする． UnCopyableClass(const UnCopyableClass& obj); // 代入演算子は使わせないようにしたいので，宣言だけにして，praivteスコープとする． UnCopyableClass& operator=(const UnCopyableClass &obj); }; int main(int, char**) { UnCopyableClass objA, objB; // コンパイルエラー！ コピーコンストラクタのコールはダメ！ //UnCopyableClass objC(objA); // コンパイルエラー！ コピー代入演算子のコールはダメ！ //objA = objB; return 0; }
ちなみに，継承関係を使ってコピーコンストラクタ・コピー代入演算子のコールを禁止することもできます．
#include <iostream> #include <string> #include <vector> class Uncopyable { public: Uncopyable() {} ~Uncopyable() {} private: Uncopyable(const Uncopyable& obj); Uncopyable& operator=(const Uncopyable& obj); }; class ChildUncopyable : private Uncopyable { public: }; int main(int, char**) { ChildUncopyable objA, objB; // コンパイルエラー！ コピーコンストラクタのコールはダメ！ // ChildUncopyable objC(objA); // コンパイルエラー！ コピー代入演算子のコールはダメ！ // objA = objB; return 0; }
確率分布からの乱数の生成２ 分布に従う乱数の発生，Box Muller法
ということで前回のエントリでは，下記の条件が既知の場合に変換後変数の確率密度分布を求める方法を考えました．
１．変換前の確率変数Xの確率密度関数がわかっている．
２．確率変数Xと確率変数Yの変換式Y=Φ(X)
で，結果的に変換後の確率変数が従う分布 g(y) と変換前の確率変数が従う分布 f(x) の間には下記の関係があることがわかりました．
一様分布に従う確率変数からの変換式
ここで，PCでは一様乱数を生成するのは簡単であることを述べました．なので，変換前の確率変数 X を 0 から 1 の確率変数に絞って考えます．このことによって，上記の式は下記のようになります．
いまの前提条件としては，Y = Φ(X) がわかっているという前提で話を進めてきました．なので，上記の変換式から変換後確率変数の分布を計算するのは簡単でした．（エントリ１の例では，y = Φ(x) を y = -log(x) として計算しましたが，x = exp(-y) とできるので，確率密度関数はすんなり求まりました．）
Y = Φ(X) の変換式（X = Φ-1(Y) ）がわからないとき
で，実用上興味があるのはたいていこちらのパターンだと思います．つまり，前提条件として，，，，
１．変換前の確率変数Xの確率密度関数がわかっている（＝一様分布）
２．変換後の確率変数Yの確率密度関数がわかっている（＝乱数を生成したい分布）
で，この二つがわかった状態で，「変換後の確率密度関数に従う乱数を，一様乱数から生成したい」というのが今回の主目的でした．なので，要は Y = Φ(X) が知りたいわけです．で，もう一度下記の数式に戻りたいのですが，
変換後の確率密度関数 p(y) を求めるには dx/dy の微分方程式を解く必要があります．が，ちょっと見方を変えてみると．．．p(y) は確率密度関数なので，その積分は累積分布関数 F(y) になります．なので，この微分方程式の解は x = F(y) と言えそうです．この時，累積分布関数 x = F(y) の逆関数 y = F-1(x) がわかれば，一様乱数から所望の分布に従う乱数を発生させる方法を手にしたことになります．この逆関数はすんなり求まるものもあれば，そうでないものもあるようです．
y = F-1(x) の幾何学的解釈
x が y の累積分布関数 F(y) という見地から見ると，一様乱数 x を入力として， y = F-1(x) を求めるということは，-∞からyまでの確率積分がxと等しくなるようなyを返しているということになります．ちょっとこんがらがってきたので，下記イラストです．
ということで，一様分布から所望の分布に従う乱数を生成する方法でした．．．．としたいところですが，具体的なところにまだ全く入ってないんですよね（笑）
Box Muller 法
上で，一様分布に従う確率変数 X から所望の分布に従う確率変数 Y を生成する方法を書きましたが，要は次の条件が満たされてないとだめでした．
１．所望の分布に従う確率変数Yを発生させるためには，Y=Φ(X) の変換式が必要．
２．所望の分布の累積分布関数（x = F(y)）の逆関数（y = F-1(x)）がわかれば，そいつが変換式になる．
で，今回乱数のケースでは，所望の分布が正規分布です．で，正規分布の累積分布関数が解析的に求まるかというと．．．．求まりません（笑）なので，結局微分方程式にうまく"ハマる"解を見つけてあげなければならず，これを見つけたのが BoxさんとMulllerさんなんです．結果だけを言うと，下記の変換をすれば，y1 と y2 はそれぞれ平均０，分散１の正規分布に従う乱数になります．
上記の考察は，，，あきらめます．代わりに，Pythonで実験です(ﾟДﾟ)ﾉ
Pythonを使って，Box-Muller 法の実験
Box-Muller法って，一様乱数が２つ必要な代わりに，生成する正規乱数も２つできるんですね．Numerical Recipes in Cには，三角関数を使わない高速バージョンが書かれてたんですが，夜も更けてきたので，一番簡単なものだけ実装しておきます．
import matplotlib.pyplot as plt import numpy as np import random import math if __name__ == '__main__': print(__file__) gen_num = 10001 bin_num = 1001 draw_range = (-5, 5) bin_width = float(((draw_range[1] - draw_range[0]) / float(bin_num))) x1 = np.array([random.random() for i in range(gen_num)]) x2 = np.array([random.random() for i in range(gen_num)]) y1 = np.array([math.sqrt(-2*math.log(x1[i])) * math.cos(2 * math.pi * x2[i]) for i in range(gen_num)]) y2 = np.array([math.sqrt(-2*math.log(x1[i])) * math.sin(2 * math.pi * x2[i]) for i in range(gen_num)]) hist1 = np.histogram(y1, bins=bin_num, range=draw_range) histelem1 = np.array([hist1[0][i] / bin_width for i in range(bin_num)]) hist2 = np.histogram(y2, bins=bin_num, range=draw_range) histelem2 = np.array([hist2[0][i] / bin_width for i in range(bin_num)]) exp_x = np.array([bin_width * i for i in range(-bin_num/2, bin_num/2)]) exp_vars = 1 / math.sqrt(2 * math.pi) * np.exp(-np.power(exp_x, 2) / 2.0) print(exp_x) #plt.plot(histelem) plt.plot(exp_x, histelem1) plt.plot(exp_x, histelem2) plt.plot(exp_x, exp_vars * gen_num) plt.show()
実行結果
参考文献
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ
- 作者: William H. Press,William T. Vetterling,Saul A. Teukolsky,Brian P. Flannery,丹慶勝市,佐藤俊郎,奥村晴彦,小林誠
- 出版社/メーカー: 技術評論社
- 発売日: 1993/06/01
- メディア: 単行本
- 購入: 2人 クリック: 102回
- この商品を含むブログ (33件) を見る
確率分布からの乱数の生成１ 確率変数の変換
ということで，下記のエントリのGMappingの解剖実験進めているんですが，Motion Model の計算のところで正規分布に従う変数を生成する必要があったので簡単にまとめておくことにしました．
確率分布と分布に従う乱数の生成
平均値μ，分散σの正規分布が下記のように書けることは確率統計の本を書けば書いてあると思います．
で上式なんですが，「考えている確率変数Xが平均μ，分散σに従うとき，取り出した値がxになる確率はいくらか？」ってことには答えてくれます．ですが，「正規分布に従う乱数を発生させる」ことにはそのまま使えません．
ここがこのエントリのテーマなんですが，ロボットのオドメトリとかセンサの計測モデルとかを考えるうえで，適当な誤差をエイやっと載せて計測の不確実性をモデルに取り込みたいときが多々あるわけです．正規分布に従う乱数発生器は Python とかならライブラリとして提供されていますが，C/C++だとどこかから探してこないといけません．
パソコンで発生させられる乱数
で，C/C++だと乱数発生器は全くないのか？という話なんですが，一様乱数を発生させる関数としていくつか関数が標準ライブラリに用意されてます．下記，0から32ビット最大値までの乱数を生成する関数のサンプルプログラムです．
#include <iostream> #include <random> int main(int, char**) { std::cout << __FILE__ << std::endl; std::mt19937 mt((int)time(0)); for (int i = 0; i < 10; i++) { std::cout << mt() << std::endl; } return 0; }
Pythonを使って実験＆プロット
C++だと可視化するのがめんどくさいので，ここからはPythonを使って実験します．上記の一様乱数をPythonを使って生成し，ヒストグラムを作ってみます．上のC++の関数は32ビットの最大値まで生成しましたが，ここでは 0 ~ 1 までの乱数を生成します．
import matplotlib.pyplot as plt import numpy as np from scipy.stats import norm import random def sample_normal_distribution(avg, var): vars = np.array([random.uniform(-var, var) for i in range(12)]) return 0.5 * np.sum(vars) if __name__ == '__main__': print(__file__) n = np.linspace(-10, 10, 201) random_nums = np.array([random.random() for i in range(1000)]) print(random_nums) plt.xlim((0, 1)) plt.ylim((0, 20)) plt.hist(random_nums, bins=101) plt.show()
一様乱数をつかった他の乱数の生成
で，問題はここからで，生成したい乱数は「正規分布に従う乱数」でした．どうやるかという話なんですが，一様乱数（＝一様分布に従う乱数）から他の分布に従う乱数を生成します．ちょっとここから数学の復習ですが．．．
確率変数Xから変換した確率変数Yの分布を求める．
分布がわかっているある確率変数をXとします．この確率変数Xにある変換をかけてY=Φ(X)を計算するとします．そうすると，このYも確率変数とみなせますが，分布が既にわかっているXの分布を使ってYの分布を計算することはできるでしょうか？XとYの確率密度関数（＝分布）をそれぞれ f(x), g(y) と表すとします．この時，，，
ということで，確率分布が既知の確率変数Xを変換してできた確率変数Yの分布を求める方法がわかりました．
変数変換の例（指数分布）
ひとつ変数変換の実際の例を考えてみたいと思います．y = Φ(x) とした変換を具体的に y = -log(x) とすると，y の分布 g(y) は先ほど求めた関係式から下記のように計算されます．
Pythonを使って確認．
ということで，上記の変数変換の結果が本当にあっているかどうか見てみます．
import matplotlib.pyplot as plt import numpy as np import random if __name__ == '__main__': print(__file__) gen_num = 10001 bin_num = 1001 draw_range = (0, 10) bin_width = float(draw_range[1] / float(bin_num)) random_nums = np.array([random.random() for i in range(gen_num)]) trans_vars = -np.log(random_nums) hist = np.histogram(trans_vars, bins=bin_num, range=draw_range) histelem = np.array([hist[0][i] / bin_width for i in range(bin_num)]) exp_x = np.array([bin_width * i for i in range(bin_num)]) exp_vars = np.exp(-exp_x) plt.plot(exp_x, histelem) plt.plot(exp_x, exp_vars * gen_num) plt.show()
上記コードで生成した結果が，下記のプロットになります．ヒストグラムの結果と確率密度関数 x 生成個数の値がいい感じで一致してます．
このエントリでわかったことは，
確率変数Xの確率密度関数 f(x) と確率変数XからYへの変換式 Y = Φ(X) がわかっている時に，どうやって確率変数Yの確率密度関数を求めるか？
でした．最終的にやりたい・知りたいことは，下記の状況で
どうやって g(y) に従う乱数（=正規乱数）を生成するか？（=一様分布からの変換式 Y = Φ(X) を求めるか）でした．
１．確率変数Xの確率密度関数 f(x) がわかっている．（一様分布）
２． f(x) に従う乱数を生成する方法をもっている．（PCによる一様乱数生成）
参考文献
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ
- 作者: William H. Press,William T. Vetterling,Saul A. Teukolsky,Brian P. Flannery,丹慶勝市,佐藤俊郎,奥村晴彦,小林誠
- 出版社/メーカー: 技術評論社
- 発売日: 1993/06/01
- メディア: 単行本
- 購入: 2人 クリック: 102回
- この商品を含むブログ (33件) を見る
