ガンマ関数についてのメモ&プログラミング初心者でもガンマ関数をプロットしたい
どうも数学とプログラミング初心者の☆ありゅ☆さんです.
今回はガンマ関数と積分を組んでみた後,ワイエルシュトラス形式のガンマ関数の対数形の計算結果との実関数としての対比を行った後に複素数平面全体に定義して描画してみたいと思います.理論は調べたので描きますし,自分は組んでないしブログ書いてないのでやります.この記事より性能の良さそうな実装については参考文献のシキノートさんのページを参照.
2020/8/10追記 漸近展開の式の途中でが抜けていたのを修正しました.
ガンマ関数ってなーに?
基本事項
ガンマ関数については
aryuaryuaryuryu.hatenablog.com
にも少し出てきたが,階乗の一般化として名高い(?)関数で,基本として
と表すことができる.
これについては部分積分を行うと
となり,また
になるため
といったように階乗を考えることができ,正の整数に対し
となる.
また
として積分をすると
が得られてこれでは0以下の整数以外について定義ができる.
階乗について他にも関数を考えることができるらしいが,対数凸という性質があるものはこのガンマ関数しかないらしく,ボーア・モーレップの定理というらしい.
このガンマ関数の有名値としては
から示せる.
として二乗すると
互いに積分は独立しているので
と変形ができる.
ここで及びと置換するとヤコビアンはになるので
は
についてと置換すれば得られることがわかる.
スターリングの公式
ガンマ関数から得られる近似公式というとスターリングの公式.
だがこれについてはLaplace's Methodという簡単な方法で得られる.
Laplace's Method
Laplace's Methodは
という積分について,がで最大値をとり双方ともを中心に減少し続け負の無限に発散していく場合,最大値周辺が一番積分に影響すると考えることができ,他の部分は0に指数関数的に近づいていくことになるので
と考えることができる.
そこでのを中心としたテイラー展開の先頭三項
と入れ替えるが,最大値を取るということはなので実際のところ2項だけになる.
こうして
と書くことができ,及びは定数なので
になり,これからガウス積分を使うと単調減少なのでこの関数はであり
と書き換え
と置換すると
と近似結果が得られる.
Laplace's Methodによるスターリングの公式の導出
こうしてLaplace's Methodを得られたのでこれをガンマ関数に適用する.
について指数関数に統一するために
として,対数を消すためにと置換すると
になる.
これでLaplace's Methodを用いると
において
になって,になるのはになるときで,常になので上に凸.
こうして及びなので
になる.最後にとして両辺にをかけると
が得られる.これがスターリングの公式だ.
ガンマ関数の関連式
ガンマ関数の無限乗積形
ガンマ関数は対数凸という性質がありが下に凸の形をしている.
そこで
という条件で整数に対し,始点がそれぞれに右にずれていき終点が同じ点についての傾き,~間,~間,~の三種の傾きの比較をする.
そこで凸関数の性質から
になる.ここで式を整理するとなので
と綺麗にまとめることができる.
こうして
と変形し,なので
つまり
とまとまる.
ガンマ関数の性質から
になるので
と変形していくと
になるので
は
と書き換えることができる.
としてやなど中途半端だから
と形式を揃える.
左右の不等式を分割して
揃えた分数の部分単体にする
そして不等式をあわせて
この式についてと極限をとってやると
これではさみうちの原理から
が得られるが,の前提があるので定義域を拡張する.
極限を取る前の式を
として
になるので
と変形しになるので
が得られ,逆に
とすると
が得られるのでこれを繰り返していくと定義域を全ての実数に拡張できて全ての実数について
が説明できる.
ガンマ関数のワイエルシュトラス表示
ガンマ関数は0以下の整数について0除算が発生することは
の式からわかる.これについて0除算を回避するために
と逆数にする.
すると
と変形ができるがこれではが邪魔くさいし相乗の部分は絶対収束しない.
ここで
と絶対収束させるためのバイアスをかけてやり
が収束するを考える.相乗に対数をつけて
ここでのテイラー展開
を使って,について
と書き換えることができるので,今後のためにとなる最小のをとおこう.
するとの範囲の総和はせいぜい有限和なのでについて考えよう.
これでを
と定義してやると
とのテイラー展開の残りの開始位置を調整することで相乗の収束・発散をコントロールすることができそうだ.
この式は絶対値をつけて
となり,現在はという条件で考えているので
と書くことができる.
こうして
について
と分割し,有限和の部分を
とすると
になる.
これで無限和は
これが収束する最小の整数は
及び
なのでのときに収束する.
ここで収束させるためのバイアスは
と定義しているので
になる.
こうして
極限を取る際オイラー定数を使って
といった式が得られる.
これがガンマ関数のワイエルシュトラス表示だ.
ガンマ関数の相反公式
ガンマ関数のワイエルシュトラス表示
にを代入して
こうしてを計算すると
こうしての無限積表示
を使うと
が得られる.
ちなみにを代入すると
が得られる.
漸近展開
ベルヌーイ多項式
ベルヌーイ多項式って何?
ベルヌーイ数によくにた定義で
のを中心としたテイラー展開
で定義される関数をベルヌーイ多項式と言い,言われてみるとわかると思うがのときはベルヌーイ数になる.
ベルヌーイ多項式の関係式
これは
についてテイラー展開に置き換えると
としてこれを
と変形し,と置換して
とつまり
の関係式が綺麗に得られる.
ちなみにのときは
といった式になる.
ベルヌーイ多項式の微分
得られた式
について微分を行うと,の時定数なので
と微分結果が得られる.
オイラー・マクローリンの定理
スターリングの公式はLaplace's Methodを使って出したが,ちゃんとした漸近展開を出すならオイラー・マクローリンの和公式を使うのが比較的楽.
オイラーマクローリンの和公式の最初の段階は,積分と幅1の長方形で積分を近似した場合(右側の区分求積法みたいな感じ)の誤差についての式で,これはまず整数に対して
についての積分を計算する.の間はなので
と書き換えて部分積分を行う.
これをからまで足していくと
が得られる.
次にこの左辺の積分を少し変えて
として
と展開してさっきの結果を入れる
式を整理すると
この結果はベルヌーイ数とベルヌーイ多項式を用いて
と表すことができ,ここで右辺の積分について及びなので
となる.これを繰り返すと正の整数に対し
が得られるので
になる.
漸近展開の計算
このオイラー・マクローリン定理を
に用いるのが良いらしい.
と変形して
と対数で分解する.はのときになるから除外して考えて
になる.
これで
についてオイラーマクローリンの定理
を用いる.
と少しだけ総和を分けて
ここで
になるので,とすると
になる.めっちゃ長い・・・これを整理して
ここで右辺について
は
と変形することで
になり,またベルヌーイ数を用いると
整理したもの
を
に入れる.
式を整理して
最後にを取る際
のとき
に注意すると
が得られる.
これは
として定数項
についてはスターリングの公式
よりのとき
がわかるので
が得られ,
とまとめることができる.これを展開すると以前求めたベルヌーイ数
aryuaryuaryuryu.hatenablog.com
を先頭10個ほど
を使うと
になる.
プログラムを組んでみる
積分形式
ガンマ関数について色々見てきたからそろそろプログラムを組んでみたい.
最初は簡単なガンマ関数の形の積分表示
を組んでみる.
積分の実装というと,直感的な方法は長方形近似
で,これだと図を見ての通り幅のとり方によっては誤差が大きくなるが,台形近似
を用いると少しはマシになる.
ちなみにこの刻み幅は幅が倍になると誤差は倍になる.
今後なんか積分計算したいときに使い回せるようにもしたいからそれも考えて実装すると
/*積分開始,積分終了,被積分関数,被積分関数の引数,積分間隔*/
double Integral(double Min, double Max, std::function<double(double,std::vector<double>)> fn, std::vector<double> v, double interval = 0.01) {
double result = 0, i;
for (i = Min; i < Max; i += interval)
result += ( (fn(i, v) + fn(i + interval, v)) * interval) / 2.;
result += ( (fn(i - interval, v) + fn(Max, v)) * (Max - i + interval)) / 2.;
return result;
}
として
についてを実装すると良いようにした.
これで
#define M_E (2.71828182845904)
double Gamma1(double t,std::vector<double> v) {
double X = v[0];
return pow(t, X - 1.) * pow(M_E, -t);
}
として計算を行うと恐ろしいほど計算時間が遅いガンマ関数が実装できた.
この実装の性能としては
■積分開始位置0.0001,
■積分終了位置10000
■台形の幅0.001
を計算するのにおおよそ1分弱,計算精度はに対し計算結果はとまぁまぁボチボチな結果が得られた.
わかってはいたけど控えめに言ってクソ.
そこでのときの被積分関数のグラフを見てみると
になる.の範囲はすごく急に変化しているのに対しその後はゆったりしての時点でめっちゃゼロ.
じゃぁ積分範囲を小さくして
■積分開始位置0.0001,
■積分終了位置5
■台形の幅0.00001
とすると結果は2.5秒ほどかかって1.76348・・・うーん・・・
ただのときは
■積分開始位置0.0001,
■積分終了位置100
■台形の幅0.001
の設定では2msでが得られた.
これは多分
のようなゆったりとしたカーブを描いているのと,台形近似による誤差で少し多めに計算されたのが積分を途中で打ち切っているのを埋めているのかなと思う.
そこでの場合は
という結構急なカーブを描いているのに720と意外とちゃんとした値が得られている.
じゃぁ多分の範囲が上に跳ね上がるのがいけないのかな?というと
今度は結果は0.99999と結構1に近い結果が得られた.やっぱり0.5のときのスタート地点がいけなかったのかな?
といろいろスタート地点等をいじったけど最終的にはやっぱ刻み幅によりそうなので断念.ただある程度大きい数値に対してはそれなりにアバウトな設定でも速くていい感じの値が得られた.
無限乗積形式
ガンマ関数の無限乗積表示
について,これを実装してみると
double Gamma3(double x,double Max = 50000) {
double result = pow(Max, x) / x;
for (int i = 1; i < 50000; i++)result *= (double)i / (i + x);
return result;
}
となる.
正直ネタのつもりだったからあらかたこの記事を書いた後で実装したら,対数形式とあまり性能は変わらなかったと思う.
正直対数と実装はあまり変わってないからそりゃそうだよね・・・
対数形式
ガンマ関数のワイエルシュトラス表示形式
に対数をつけて
となるので,を計算して後でにぶちこめば良いのでは?
ということで
double Gamma2(double x) {
double result = -log(x) - x * M_GAMMA;
for (int i = 1; i < 50000; i++)result -= log(1. + x / i) - x / i;
return pow(M_E,result);
}
と実装した.
同様にを計算するのにおおよそ0.002秒,計算精度はに対し計算結果はが得られた.さっきの積分形式に比べると恐ろしく性能がよい.
そこでx=5を代入した場合,理論値は24だが23.994が出てきて誤差0.006に.
x=9を代入した結果理論値は362880だが結果は40287.4にと代入値が大きくなると誤差が大きくなっていった.
これについて考えてみると
についてのときのテイラー展開
を利用して
と変形して
と書き換えることができる.
が1未満だとこのの内部の総和が早く収束するので,引数が小さい場合は結構いい感じの値が得やすいのかな?と思う.
そこで1未満の値の場合はそのまま出力して,1超過の場合は
の函数等式を使うのが良さ気で,0未満の場合は
を使えそう.
ただこれは実軸方向への拡張なので,虚数方向に移動するとなると別の方法が必要になる.
漸近展開による方法
ベルヌーイ数は予め求めてたからいきなり値ぶっこんでおくw
そこで漸近展開を用いて計算するなら簡単に書くと
double Gamma4(double x) {
std::vector<double> Bernoulli;
Bernoulli.push_back(1./12.); Bernoulli.push_back(0);
Bernoulli.push_back(-1./360.); Bernoulli.push_back(0);
Bernoulli.push_back(1./1260.); Bernoulli.push_back(0);
Bernoulli.push_back(-1./1680.); Bernoulli.push_back(0);
Bernoulli.push_back(1./1188.); Bernoulli.push_back(0);
Bernoulli.push_back(-691./360360.); Bernoulli.push_back(0);
Bernoulli.push_back(0); Bernoulli.push_back(1./156.);
double Result = sqrt(2. * M_PI / x) * pow(x / M_E, x);
double tmp = 0;
for (int i = 1; i <= Bernoulli.size(); i++) {
double temp=1;
for (int j = 0; j < i; j++)temp *= x;
tmp += Bernoulli[i - 1]/temp;
}
return Result * exp(tmp);
}
になる.
引数が1以上なら結構いい具合の値が出てくる.1以上に対しての計算はもうこれでいいのではないか
絶対値が1未満なら対数形式で,1超過ならこれで十分かも.
複素数平面全体について定義する
基本的に今まで作った型の引数等を『std::complex<double>』に変えてやるだけ.
0以下の整数はとりあえず計算最大値INFTY(適当に決める)を返して,正の整数は階乗返しとけば良いかも.
引数の絶対値が1未満では対数形式,1超過だと漸近展開形式が良さ気.絶対値が1の場合は(1+i)を代入して比較した感じ漸近展開形式がマシだろうはと思う.
ゼータ関数と同じようにDxlibで描画.は0~4の範囲で,及びは-4~4の範囲で描画してみた.
あとは偏角も描画してみた.こっちはHSVのHを偏角として,S=V=1で描画した.