(-\infty, a]累積分布関数は閉じた形では利用できないため、から確率密度関数を統合したいと思います。しかし、C++でこれを行う方法がわかりません。
このタスクはMathematicaではとても簡単です。私がする必要があるのは、関数を定義することだけです。
f[x_, lambda_, alpha_, beta_, mu_] :=
Module[{gamma},
gamma = Sqrt[alpha^2 - beta^2];
(gamma^(2*lambda)/((2*alpha)^(lambda - 1/2)*Sqrt[Pi]*Gamma[lambda]))*
Abs[x - mu]^(lambda - 1/2)*
BesselK[lambda - 1/2, alpha Abs[x - mu]] E^(beta (x - mu))
];
次に、ルーチンを呼び出しNIntegrateて数値積分します。
F[x_, lambda_, alpha_, beta_, mu_] :=
NIntegrate[f[t, lambda, alpha, beta, mu], {t, -\[Infinity], x}]
今、私はC++で同じことを達成したいと思います。gsl_integration_qagilgsl数値ライブラリのルーチンを使用しています。それは私が欲しいものである半無限の間隔(-\infty, a]で機能を統合するように設計されています。しかし、残念ながら私はそれを機能させることができません。
これはC++の密度関数です。
density(double x)
{
using namespace boost::math;
if(x == _mu)
return std::numeric_limits<double>::infinity();
return pow(_gamma, 2*_lambda)/(pow(2*_alpha, _lambda-0.5)*sqrt(_pi)*tgamma(_lambda))* pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu));
}
次に、gslルーチンを呼び出して、統合してcdfを取得しようとします。
cdf(double x)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
double result, error;
gsl_function F;
F.function = &density;
double epsabs = 0;
double epsrel = 1e-12;
gsl_integration_qagil (&F, x, epsabs, epsrel, 1000, w, &result, &error);
printf("result = % .18f\n", result);
printf ("estimated error = % .18f\n", error);
printf ("intervals = %d\n", w->size);
gsl_integration_workspace_free (w);
return result;
}
ただしgsl_integration_qagil、エラーを返しますnumber of iterations was insufficient。
double mu = 0.0f;
double lambda = 3.0f;
double alpha = 265.0f;
double beta = -5.0f;
cout << cdf(0.01) << endl;
ワークスペースのサイズを大きくすると、ベッセル関数は評価されません。
私の問題について何か洞察を与えてくれる人はいないかと思っていました。上記の対応する数学関数Fを呼び出すと、がx = 0.01返されます0.904384。
密度が非常に小さな間隔に集中している可能性があります（つまり[-0.05, 0.05]、密度の外側はほぼ0です。プロットを以下に示します）。もしそうなら、これについて何ができるか。読んでくれてありがとう。
