mql-for-begginers/Include/Math/Stat/Gamma.mqh
2025-07-22 18:30:17 +03:00

764 lines
31 KiB
MQL5

//+------------------------------------------------------------------+
//| Gamma.mqh |
//| Copyright 2000-2025, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
#include "Normal.mqh"
const double DoubleEpsilon=1.11022302462515654042E-16;
const double LogMax=7.09782712893383996732E2;
//+------------------------------------------------------------------+
//| Inverse of the incomplete Gamma integral |
//+------------------------------------------------------------------+
double MathInverseGammaIncomplete(const double a,const double y)
{
//--- bound the solution
double x0 = DBL_MAX;
double yl = 0;
double x1 = 0;
double yh = 1.0;
double dithresh=5.0*DoubleEpsilon;
//--- approximation to inverse function
double d=1.0/(9.0*a);
int err_code=0;
double q_normal=MathQuantileNormal(y,0,1,true,false,err_code);
double yy=(1.0-d-q_normal*MathSqrt(d));
double x=a*yy*yy*yy;
double lgm=MathGammaLog(a);
for(int i=0; i<10; i++)
{
if(x>x0 || x<x1)
break;
yy=1.0-MathGammaIncomplete(x,a);
if(yy<yl || yy>yh)
break;
if(yy<y)
{
x0 = x;
yl = yy;
}
else
{
x1 = x;
yh = yy;
}
//--- compute the derivative of the function at this point
d=(a-1.0)*MathLog(x)-x-lgm;
if(d<-LogMax)
break;
d=-MathExp(d);
//--- compute the step to the next approximation of x
d=(yy-y)/d;
if(MathAbs(d/x)<DoubleEpsilon)
return (x);
x=x-d;
}
//--- resort to interval halving if Newton iteration did not converge.
d=0.0625;
if(x0==DBL_MAX)
{
if(x<=0.0)
x=1.0;
while(x0==DBL_MAX && MathIsValidNumber(x))
{
x=(1.0+d)*x;
yy=1.0-MathGammaIncomplete(x,a);
if(yy<y)
{
x0 = x;
yl = yy;
break;
}
d=d+d;
}
}
d=0.5;
double dir=0;
for(int i=0; i<400; i++)
{
double t=x1+d *(x0-x1);
if(!MathIsValidNumber(t))
break;
x=t;
yy=1.0-MathGammaIncomplete(x,a);
lgm=(x0-x1)/(x1+x0);
if(MathAbs(lgm)<dithresh)
break;
lgm=(yy-y)/y;
if(MathAbs(lgm)<dithresh)
break;
if(x<=0.0)
break;
if(yy>=y)
{
x1 = x;
yh = yy;
if(dir<0)
{
dir=0;
d=0.5;
}
else
if(dir>1)
d=0.5*d+0.5;
else
d=(y-yl)/(yh-yl);
dir+=1;
}
else
{
x0 = x;
yl = yy;
if(dir>0)
{
dir=0;
d=0.5;
}
else
if(dir<-1)
d=0.5*d;
else
d=(y-yl)/(yh-yl);
dir-=1;
}
}
if(x==0.0 || !MathIsValidNumber(x))
{
Print("Errors in an arithmetic, casting, or conversion operation.");
return(QNaN);
}
//---
return(x);
}
//+------------------------------------------------------------------+
//| Gamma probability density function (PDF) |
//+------------------------------------------------------------------+
//| The function returns the probability density function of |
//| of the Gamma distribution with shape parameters a and b. |
//| |
//| Arguments: |
//| x : Random variable |
//| a : Shape |
//| b : Scale |
//| log_mode : Logarithm mode flag, if true it returns Log values |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The probability density evaluated at x. |
//+------------------------------------------------------------------+
double MathProbabilityDensityGamma(const double x,const double a,const double b,const bool log_mode,int &error_code)
{
//--- check parameters
if(!MathIsValidNumber(x) || !MathIsValidNumber(a) || !MathIsValidNumber(b))
{
error_code=ERR_ARGUMENTS_NAN;
return QNaN;
}
//--- a and b must be positive
if(a<=0 || b<=0)
{
error_code=ERR_ARGUMENTS_INVALID;
return QNaN;
}
error_code=ERR_OK;
//--- check negative x
if(x<=0)
return TailLog0(true,log_mode);
//--- calculate log Gamma density
double log_result=(a-1.0)*MathLog(x)-(x/b)-MathGammaLog(a)-a*MathLog(b);
if(log_mode==true)
return(log_result);
//--- return Gamma density
return MathExp(log_result);
}
//+------------------------------------------------------------------+
//| Gamma probability density function (PDF) |
//+------------------------------------------------------------------+
//| The function returns the probability density function of |
//| of the Gamma distribution with shape parameters a and b. |
//| |
//| Arguments: |
//| x : Random variable |
//| a : Shape |
//| b : Scale |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The probability density evaluated at x. |
//+------------------------------------------------------------------+
double MathProbabilityDensityGamma(const double x,const double a,const double b,int &error_code)
{
return MathProbabilityDensityGamma(x,a,b,false,error_code);
}
//+------------------------------------------------------------------+
//| Gamma probability density function (PDF) |
//+------------------------------------------------------------------+
//| The function calculates the Gamma probability density function |
//| with parameters a and b for values in x[] array. |
//| |
//| Arguments: |
//| x : Array with random variables |
//| a : Shape |
//| b : Scale |
//| log_mode : Logarithm mode flag,if true it calculates Log values|
//| result : Output array for calculated values |
//| |
//| Return value: |
//| true if successful, otherwise false. |
//+------------------------------------------------------------------+
bool MathProbabilityDensityGamma(const double &x[],const double a,const double b,const bool log_mode,double &result[])
{
//--- check parameters
if(!MathIsValidNumber(a) || !MathIsValidNumber(b))
return false;
//--- a and b must be positive
if(a<=0 || b<=0)
return false;
int data_count=ArraySize(x);
if(data_count==0)
return false;
int error_code=0;
ArrayResize(result,data_count);
for(int i=0; i<data_count; i++)
{
double x_arg=x[i];
if(!MathIsValidNumber(x_arg))
return false;
if(x_arg>0)
{
//--- calculate log Gamma density
double log_result=(a-1.0)*MathLog(x_arg)-(x_arg/b)-MathGammaLog(a)-a*MathLog(b);
if(log_mode==true)
result[i]=log_result;
else
result[i]=MathExp(log_result);
}
else
result[i]=TailLog0(true,log_mode);
}
return true;
}
//+------------------------------------------------------------------+
//| Gamma probability density function (PDF) |
//+------------------------------------------------------------------+
//| The function calculates the Gamma probability density function |
//| with parameters a and b for values from x[] array. |
//| |
//| Arguments: |
//| x : Array with random variables |
//| a : Shape |
//| b : Scale |
//| result : Array with calculated values |
//| |
//| Return value: |
//| true if successful, otherwise false. |
//+------------------------------------------------------------------+
bool MathProbabilityDensityGamma(const double &x[],const double a,const double b,double &result[])
{
return MathProbabilityDensityGamma(x,a,b,false,result);
}
//+------------------------------------------------------------------+
//| Gamma cumulative distribution function (CDF) |
//+------------------------------------------------------------------+
//| The function returns the cumulative distribution function of the |
//| Gamma distribution with parameters a and b. |
//| |
//| Arguments: |
//| x : The desired quantile |
//| a : Shape |
//| b : Scale |
//| tail : Flag to calculate lower tail |
//| log_mode : Logarithm mode flag,if true it calculates Log values|
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The value of the Gamma cumulative distribution function |
//| with parameters a and b, evaluated at x. |
//+------------------------------------------------------------------+
double MathCumulativeDistributionGamma(const double x,const double a,const double b,const bool tail,const bool log_mode,int &error_code)
{
//--- check NaN
if(!MathIsValidNumber(x) || !MathIsValidNumber(a) || !MathIsValidNumber(b))
{
error_code=ERR_ARGUMENTS_NAN;
return QNaN;
}
//--- a and b must be positive
if(a<=0 || b<=0)
{
error_code=ERR_ARGUMENTS_INVALID;
return QNaN;
}
error_code=ERR_OK;
//--- check x
if(x<=0)
return TailLog0(tail,log_mode);
//--- calculate probability using Incomplete Gamma function and take into account round-off errors
double cdf=MathMin(MathGammaIncomplete(x/b,a),1.0);
return TailLogValue(cdf,tail,log_mode);
}
//+------------------------------------------------------------------+
//| Gamma cumulative distribution function (CDF) |
//+------------------------------------------------------------------+
//| The function returns the cumulative distribution function of the |
//| Gamma distribution with parameters a and b. |
//| |
//| Arguments: |
//| x : The desired quantile |
//| a : Shape |
//| b : Scale |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The value of the Gamma cumulative distribution function |
//| with parameters a and b, evaluated at x. |
//+------------------------------------------------------------------+
double MathCumulativeDistributionGamma(const double x,const double a,const double b,int &error_code)
{
return MathCumulativeDistributionGamma(x,a,b,true,false,error_code);
}
//+------------------------------------------------------------------+
//| Gamma cumulative distribution function (CDF) |
//+------------------------------------------------------------------+
//| The function calculates the values of the Gamma cumulative |
//| distribution function with given a and b for values in x[] array.|
//| |
//| Arguments: |
//| x : Array with random variables |
//| a : Shape |
//| b : Scale |
//| tail : Flag to calculate lower tail |
//| log_mode : Logarithm mode flag,if true it calculates Log values|
//| resut : Output array for calculated values |
//| |
//| Return value: |
//| true if successul, otherwise false. |
//+------------------------------------------------------------------+
bool MathCumulativeDistributionGamma(const double &x[],const double a,const double b,const bool tail,const bool log_mode,double &result[])
{
//--- check NaN
if(!MathIsValidNumber(a) || !MathIsValidNumber(b))
return false;
//--- a and b must be positive
if(a<=0 || b<=0)
return false;
int data_count=ArraySize(x);
if(data_count==0)
return false;
int error_code=0;
ArrayResize(result,data_count);
for(int i=0; i<data_count; i++)
{
double x_arg=x[i];
if(!MathIsValidNumber(x_arg))
return false;
if(x_arg<=0)
result[i]=TailLog0(tail,log_mode);
else
{
//--- calculate probability using Incomplete Gamma function and take into account round-off errors
double cdf=MathMin(MathGammaIncomplete(x_arg/b,a),1.0);
result[i]=TailLogValue(cdf,tail,log_mode);
}
}
return true;
}
//+------------------------------------------------------------------+
//| Gamma cumulative distribution function (CDF) |
//+------------------------------------------------------------------+
//| The function calculates the values of the Gamma cumulative |
//| distribution function with given a and b for values in x[] array.|
//| |
//| Arguments: |
//| x : Array with random variables |
//| a : Shape |
//| b : Scale |
//| result : Output array for calculated values |
//| |
//| Return value: |
//| true if successul, otherwise false. |
//+------------------------------------------------------------------+
bool MathCumulativeDistributionGamma(const double &x[],const double a,const double b,double &result[])
{
return MathCumulativeDistributionGamma(x,a,b,true,false,result);
}
//+------------------------------------------------------------------+
//| Gamma distribution quantile function (inverse CDF) |
//+------------------------------------------------------------------+
//| The function returns the inverse cumulative distribution |
//| function of the Gamma distribution with parameters a and b |
//| for the desired probability. |
//| |
//| Arguments: |
//| probability : The desired probability |
//| a : Shape |
//| b : Scale |
//| tail : Flag to calculate lower tail |
//| log_mode : Logarithm mode, if true it calculates Log values |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The value of the inverse cumulative distribution function |
//| of the Gamma distribution with parameters a and b. |
//+------------------------------------------------------------------+
double MathQuantileGamma(const double probability,const double a,const double b,const bool tail,const bool log_mode,int &error_code)
{
//--- case log probability==-inf
if(log_mode==true && probability==QNEGINF)
{
error_code=ERR_OK;
return 0.0;
}
//--- check NaN
if(!MathIsValidNumber(probability) || !MathIsValidNumber(a) || !MathIsValidNumber(b))
{
error_code=ERR_ARGUMENTS_NAN;
return QNaN;
}
//--- a and b must be positive
if(a<=0 || b<=0)
{
error_code=ERR_ARGUMENTS_INVALID;
return QNaN;
}
//--- calculate real probability
double prob=TailLogProbability(probability,tail,log_mode);
//--- check probability range
if(prob<0.0 || prob>1.0)
{
error_code=ERR_ARGUMENTS_INVALID;
return QNaN;
}
error_code=ERR_OK;
//--- case probability==0
if(prob==0.0)
return 0.0;
//--- case probability==1
if(prob==1.0)
{
error_code=ERR_RESULT_INFINITE;
return QPOSINF;
}
//--- calculate quantile
double quantile=MathInverseGammaIncomplete(a,1.0-prob)*b;
if(!MathIsValidNumber(quantile))
error_code=ERR_NON_CONVERGENCE;
//---
return(quantile);
}
//+------------------------------------------------------------------+
//| Gamma distribution quantile function (inverse CDF) |
//+------------------------------------------------------------------+
//| The function returns the inverse cumulative distribution |
//| function of the Gamma distribution with parameters a and b |
//| for the desired probability. |
//| |
//| Arguments: |
//| probability : The desired probability |
//| a : Shape |
//| b : Scale |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The value of the inverse cumulative distribution function |
//| of the Gamma distribution with parameters a and b. |
//+------------------------------------------------------------------+
double MathQuantileGamma(const double probability,const double a,const double b,int &error_code)
{
return MathQuantileGamma(probability,a,b,true,false,error_code);
}
//+------------------------------------------------------------------+
//| Gamma distribution quantile function (inverse CDF) |
//+------------------------------------------------------------------+
//| The function returns the inverse cumulative distribution |
//| function of the Gamma distribution with parameters a and b |
//| for values from the probability[] array. |
//| |
//| Arguments: |
//| probability : Array with probabilities |
//| a : Shape |
//| b : Scale |
//| tail : Flag to calculate lower tail |
//| log_mode : Logarithm mode,if true it calculates for Log values|
//| result : Output array for calculated values |
//| |
//| Return value: |
//| true if successful, otherwise false. |
//+------------------------------------------------------------------+
bool MathQuantileGamma(const double &probability[],const double a,const double b,const bool tail,const bool log_mode,double &result[])
{
//--- check NaN
if(!MathIsValidNumber(a) || !MathIsValidNumber(b))
return false;
//--- a and b must be positive
if(a<=0 || b<=0)
return false;
int data_count=ArraySize(probability);
if(data_count==0)
return false;
int error_code=0;
ArrayResize(result,data_count);
const double eps=10E-18;
double max_h=MathSqrt(eps);
const int max_iterations=1000;
for(int i=0; i<data_count; i++)
{
//--- calculate real probability
double prob=TailLogProbability(probability[i],tail,log_mode);
if(!MathIsValidNumber(prob))
return false;
//--- check probability range
if(prob<0.0 || prob>1.0)
return false;
//--- case probability==0
if(prob==0.0)
result[i]=0.0;
else
//--- case probability==1
if(prob==1.0)
result[i]=QPOSINF;
else
{
double quantile=MathInverseGammaIncomplete(a,1.0-prob)*b;
if(MathIsValidNumber(quantile))
result[i]=quantile;
else
return false;
}
}
return true;
}
//+------------------------------------------------------------------+
//| Gamma distribution quantile function (inverse CDF) |
//+------------------------------------------------------------------+
//| The function returns the inverse cumulative distribution |
//| function of the Gamma distribution with parameters a and b |
//| for the desired probability. |
//| |
//| Arguments: |
//| probability : The desired probability |
//| a : Shape |
//| b : Scale |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| true if successful, otherwise false. |
//+------------------------------------------------------------------+
bool MathQuantileGamma(const double &probability[],const double a,const double b,double &result[])
{
return MathQuantileGamma(probability,a,b,true,false,result);
}
//+------------------------------------------------------------------+
//| Random variate from the Gamma distribution |
//+------------------------------------------------------------------+
//| Compute the random variable from the Gamma distribution |
//| with parameters a and b. |
//| |
//| Arguments: |
//| a : Shape |
//| b : Scale |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The random value with Gamma distribution. |
//+------------------------------------------------------------------+
//| Author: Robert Kern |
//+------------------------------------------------------------------+
double MathRandomGamma(const double a,const double b)
{
double bb,c,U,V,X=0,Y;
//--- check shape
if(a==1.0)
{
//--- exponential
return -MathLog(1.0-MathRandomNonZero());
}
else
if(a<1.0)
{
for(;;)
{
U=MathRandomNonZero();
//--- exponential
V=-MathLog(1.0-MathRandomNonZero());
if(U<=1.0-a)
{
X=MathPow(U,1.0/a);
if(X<=V)
return b*X;
}
else
{
Y = -MathLog((1-U)/a);
X = MathPow(1.0 - a + a*Y, 1.0/a);
if(X<=(V+Y))
return(b*X);
}
}
}
else
{
bb= a-1.0/3.0;
c = 1.0/MathSqrt(9*bb);
for(;;)
{
do
{
//--- generate normal random variate
double f,x1,x2,r2;
do
{
x1=2.0*MathRandomNonZero()-1.0;
x2=2.0*MathRandomNonZero()-1.0;
r2=x1*x1+x2*x2;
}
while(r2>=1.0 || r2==0.0);
//--- Box-Muller transform
f=MathSqrt(-2.0*MathLog(r2)/r2);
X=f*x2;
V=1.0+c*X;
}
while(V<=0.0);
V = V*V*V;
U = MathRandomNonZero();
if(U<1.0-0.0331*(X*X)*(X*X))
return(bb*V*b);
if(MathLog(U)<0.5*X*X+bb*(1.0-V+MathLog(V)))
return(bb*V*b);
}
}
return(X*b);
}
//+------------------------------------------------------------------+
//| Random variate from the Gamma distribution |
//+------------------------------------------------------------------+
//| Compute the random variable from the Gamma distribution |
//| with parameters a and b. |
//| |
//| Arguments: |
//| a : Shape |
//| b : Scale |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| The random value with Gamma distribution. |
//+------------------------------------------------------------------+
double MathRandomGamma(const double a,const double b,int &error_code)
{
//--- check NaN
if(!MathIsValidNumber(a) || !MathIsValidNumber(b))
{
error_code=ERR_ARGUMENTS_NAN;
return QNaN;
}
//--- a and b must be positive
if(a<=0 || b<=0)
{
error_code=ERR_ARGUMENTS_INVALID;
return QNaN;
}
error_code=ERR_OK;
return MathRandomGamma(a,b);
}
//+------------------------------------------------------------------+
//| Random variate from the Gamma distribution |
//+------------------------------------------------------------------+
//| The function generates random variables from the Gamma |
//| distribution with parameters a and b. |
//| |
//| Arguments: |
//| a : First shape parameter (a>0) |
//| b : Second shape parameter (b>0) |
//| data_count : Number of values needed |
//| result : Output array for random values |
//| |
//| Return value: |
//| true if successful, otherwise false. |
//+------------------------------------------------------------------+
bool MathRandomGamma(const double a,const double b,const int data_count,double &result[])
{
if(data_count<=0)
return false;
//--- check NaN
if(!MathIsValidNumber(a) || !MathIsValidNumber(b))
return false;
//--- a and b must be positive
if(a<=0 || b<=0)
return false;
//--- prepare output array and calculate values
ArrayResize(result,data_count);
for(int i=0; i<data_count; i++)
result[i]=MathRandomGamma(a,b);
return true;
}
//+------------------------------------------------------------------+
//| Gamma distribution moments |
//+------------------------------------------------------------------+
//| The function calculates 4 first moments of Gamma distribution |
//| with parameters a and b. |
//| |
//| Arguments: |
//| a : Shape |
//| b : Scale |
//| mean : Variable for mean value (1st moment) |
//| variance : Variable for variance value (2nd moment) |
//| skewness : Variable for skewness value (3rd moment) |
//| kurtosis : Variable for kurtosis value (4th moment) |
//| error_code : Variable for error code |
//| |
//| Return value: |
//| true if moments calculated successfully, otherwise false. |
//+------------------------------------------------------------------+
bool MathMomentsGamma(const double a,const double b,double &mean,double &variance,double &skewness,double &kurtosis,int &error_code)
{
//--- default values
mean =QNaN;
variance=QNaN;
skewness=QNaN;
kurtosis=QNaN;
//--- check NaN
if(!MathIsValidNumber(a) || !MathIsValidNumber(b))
{
error_code=ERR_ARGUMENTS_NAN;
return false;
}
//--- a and b must be positive
if(a<=0 || b<=0)
{
error_code=ERR_ARGUMENTS_INVALID;
return false;
}
error_code=ERR_OK;
//--- calculate moments
mean =a*b;
variance=a*b*b;
skewness=2/MathSqrt(a);
kurtosis=6/a;
//--- successful
return true;
}
//+------------------------------------------------------------------+