QuantLib::Mathtool::solver

今日はsolverについてです.
f(x) = 0 を解くsolverです.
クラスには複数種類あってざっとこんな感じ.

Brent,Bisection,Secant,Ridder,Newton,Falsepostion.

●アルゴリズム
Brent法
http://en.wikipedia.org/wiki/Brent's_method
Secant
http://en.wikipedia.org/wiki/Secant_method
Ridder,Falesepoint
後日(ググってもパットはでてこなかった)

どれもコンストラクタに引数はなし.
どのクラスにも2種類のsolve関数が用意されている.
Real solve ( const F& f, Real accuracy , Real guess ,Real step )
Real solve ( const F& f, Real accuracy , Real guess ,Real xMin , Real xMax )
step,xMin,xMaxは名前の通り.accuracyはソルバーのクラスによって意味が異なってくる.コードをみると,たとえばNewtonでは
    if (std::fabs(dx) < xAccuracy)return root_; 
ってあるのでxの変化がaccuracyを越えなければ終わりって感じですね.

f は解きたい関数.Newton以外のソルバーを使う場合,Fは演算子double operator(double) const を持っていれば十分.Newtonを使う時だけは double derivative(double) constも持っていなければならない.(実際solveのコードの部分を見てみるとf.derivative(_root)という形でderivativeを呼び出し,なければエラーを投げる形になっている)

実際にこれを使って円周率を計算させてみます.(sin(x)=0の解をもとめる).
コードはこんな感じになりました.
using namespace QuantLib;

class SquareClass{ public: //operatorとderivativeの定義 double operator()(const double x)const{ return std::sin(x); } double derivative(const double x)const{ return std::cos(x); } }; void testFunction3(){ //Brent法とNewton法を使う Newton newton;Brent brent; //各種パラメータ設定 Real accuracy =0.00001 , guess =2.0; Real min =1.0 , max =4.0; //関数(オブジェト)生成 SquareClass squareinst; //Brent法とNewton法で実際にとく.Brent法だけならderivativeは要らない. Real ans = newton.solve(squareinst,accuracy,guess,min,max); Real ans1 = brent.solve(squareinst,accuracy,guess,min,max); std::cout << "answer is " << ans << std::endl; std::cout << "answer is " << ans1<< std::endl; }
結果
resultもちろんBrent法だけなら,
double f (double x) {return std::sin(x);} 
って関数つくってbrent.solve(f,accuracy,・・・とすることもできる. 


今日から http://quantlib.org/docs.shtml のQuantLib introduction, part IIをもとにQuantLibを扱えるようになる修行を積みたいと思います.C++も全然わかってないので,目標を,

①QuantLibを使いこなせる 
②C++の文法理解
③各種ドキュメントを読めるようになる.

にしてやりたいと思います.本当は何か金融商品のお話を読んで,その中に出てくる計算を実際にやってみるとかのほうが楽しいのですが,どのいい本が良いのかとか全く分からないので上のPDFを読んでみることにしました.

 今回はまずMathmatical Tools のIntegrationです.
この節では積分についてのお話です.まずコードです.
Real normalPdf(const Real& x,const Real& a,const Real& b){
	boost::math::normal_distribution<> d;
	Real t1 = 0.5*(b-a),t2 = 0.5*(b+a);
	return t1*pdf(d,t1*x+t2);
}

void testIntegration(){
	//積分区間の決定
	Real a = -1.96, b = 1.96;
	//関数の作成
	boost::function < double > myPdf;
	myPdf = boost::bind(normalPdf,_1,a,b);
	//積分の準備
	GaussChebyshevIntegration gChebInt(64);
	//結果の表示
	std::cout << gChebInt(myPdf) << std::endl;

}
ざっくりみるとboost::bindとboost::functionを使って関数を持って,GaussChebyshevIntegrationで積分するって感じですね.
細かく見ていきたいと思います.

1.Real 
これは内部的にはdoubleで持っています.QunatLibでは内部的にはdoubleでも型を変えています.これは引数に間違ったものを入れないようにするための親切設計だそうです.(たとえばをボラティリティを引数に取るのにまちがって平均を入れないように,Volatilityという型が用意されていたりします.詳しくは QuantLib introduction, part I にて.)
 
2. 
boost::math::normal_distribution<>
QuantLibではboostライブラリを使用しています. そこでboostのドキュメントの対応する部分(http://www.boost.org/doc/libs/1_47_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/normal_dist.html)を見るとテンプレート引数が省略されたときはdoubleで値を持っています.(もう一つの引数policyは詳しいことは分かっていません)

3.boost::bind
これもドキュメントを見てみます(http://www.boost.org/doc/libs/1_47_0/libs/bind/bind.html).
説明をみると,”boost::bind is a generalization of the standard functions std::bind1st and std::bind2nd”
だそうです.matlabの無名関数みたいなやつですね,たぶん.これは便利._1で元の関数の引数の1番目を残しておき後の引数に値を代入したものをboost::functionでもって置きます.

4.boost::function
ドキュメントを見ます.http://www.boost.org/doc/libs/1_47_0/doc/html/function/tutorial.html
boost::function<T1, (T2)>で返り値T1型,引数T2型です.引数は複数でもよい.T1 operator()(T2)が定義されていればよい.

5.GaussChebyshevIntegration
Gauss-Chebyshev積分法についてはこれを→http://keisan.casio.jp/exec/system/1360721886.コンストラクタの引数は選点数です.積分の種類はいろいろあります(PDF参照)

初回なのでこのくらいで.もうちょっとドキュメントを読み込んで面白い使い方とかみたいですね・・・では.
 

C++のSTL(Standard Template Library)を使いたいなぁと漠然と思ったので演習問題を解いてみました.
なぜかというと,個人的に当面の目標としてQuantLib http://quantlib.org/index.shtml を使えると幸せかなと思い,そのためにtemplateに慣れようということでSTLをいじることにしました.

STLの解説は以下のページと何冊かの本を参考にしました.
・Standard Template Library プログラミング on the Web
http://episteme.wankuma.com/stlprog/

・STL演習
http://plaza.harmonix.ne.jp/~fakira/cppdoc/stl.htm

今回はSTL演習(http://plaza.harmonix.ne.jp/~fakira/cppdoc/stl.htm)の中の演習問題2を手を動かしてみました.
「3×3のマスに1から9までの数字を入れて、縦・横・斜めの和が全て 15になるものを全て求めてください。」

方針としては,
①3×3のマスに1から9を配置する.全ての配置について考える.
②各配置に対して縦・横・ななめの和が15になっているのかチェックする.
という感じです.(方針も何もない)

①について
これはstd::vectorで値を持ちました.std::permutationという関数が1~9までの順列を作ってくれます. ちなみにもとのvecはクラスGenで初期化(?)してますが,これは演算子()を定義して書いてます.これを考え付いた人(過去の偉人?)は偉い気がします.

②について
ナイーブに書きました.このままだとマスがN×Nに対応できない. いい方法があればどなたかお願いします. 全体のコードはこんな感じです.
//3×3のマス目に1-9の値を入れて,各列,各行,各ななめの和がすべて15になるような
//ものを求める

class Gen {
	//1,2,3,という順に要素を代入する
		int	i;
	public:
		Gen(){i = 0;}
		int operator()(){return ++i;}
	};

int _tmain(int argc, _TCHAR* argv[])
{
	std::vector vec(9);
	//1,2.3・・・配列の作成
	std::generate(vec.begin(),vec.end(),Gen());
	//全ての順列の作成
	while(std::next_permutation(vec.begin(),vec.end())){
		int check[8];
		//このwhileではvecが順次置換されていく
		//各行,列,ななめの和の計算
		check[0] = vec[0] + vec[1] + vec[2];
		check[1] = vec[3] + vec[4] + vec[5];
		check[2] = vec[6] + vec[7] + vec[8];
		check[3] = vec[0] + vec[3] + vec[6];
		check[4] = vec[1] + vec[4] + vec[7];
		check[5] = vec[2] + vec[5] + vec[8];
		check[6] = vec[0] + vec[4] + vec[8];
		check[7] = vec[2] + vec[4] + vec[6];

		//先の各和が15になっているかチェック
		int* i = std::find_if( check , check+8 , std::bind1st( std::not_equal_to() , 15 ) );
		if(i == check + 8){
			for(unsigned short j=0;j<9;j++){
				std::cout << vec[j] ;
			}
			std::cout << "\n";
		}
	}

	getchar();
	return 0;
}
なにかアドバイスいただけるととても嬉しいです.レイアウトが汚いのでなんとかします.  次回はもうちょっと高度なものを書きたいです. では!

このページのトップヘ

見出し画像
×