874 lines
34 KiB
MQL5
874 lines
34 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| Binomial.mqh |
|
|
//| Copyright 2000-2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#include "Math.mqh"
|
|
#include "Beta.mqh"
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial probability mass function (PDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the Binomial probability mass function |
|
|
//| with parameters n and p at x. |
|
|
//| |
|
|
//| f(x,n,p)= C(n,x)*(p^x)*(1-p)^(n-x) |
|
|
//| |
|
|
//| where binomial coefficient C(n,k)=n!/(k!*(n-k)!) |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Integer random variable |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| log_mode : Logarithm mode flag, if true it returns Log values |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| The probability mass function evaluated at x. |
|
|
//+------------------------------------------------------------------+
|
|
double MathProbabilityDensityBinomial(const double x,const double n,const double p,const bool log_mode,int &error_code)
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(x) || !MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
{
|
|
error_code=ERR_ARGUMENTS_NAN;
|
|
return QNaN;
|
|
}
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
//--- check p range
|
|
if(p<0.0 || p>1.0)
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
|
|
error_code=ERR_OK;
|
|
//--- case p=0
|
|
if(p==0.0 || p==1.0)
|
|
return TailLog0(true,log_mode);
|
|
//--- check x range
|
|
if(x<0 || x>n)
|
|
return TailLog0(true,log_mode);
|
|
|
|
double log_result=MathGammaLog(n+1.0)-MathGammaLog(x+1.0)-MathGammaLog(n-x+1.0)+x*MathLog(p)+(n-x)*MathLog(1.0-p);
|
|
if(log_mode==true)
|
|
return log_result;
|
|
//--- return probability mass
|
|
return MathExp(log_result);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial probability mass function (PDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the Binomial probability mass function |
|
|
//| with parameters n and p at x. |
|
|
//| |
|
|
//| f(x,n,p)= C(n,x)*(p^x)*(1-p)^(n-x) |
|
|
//| |
|
|
//| where binomial coefficient C(n,k)=n!/(k!*(n-k)!) |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Integer random variable |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| The probability mass function evaluated at x. |
|
|
//+------------------------------------------------------------------+
|
|
double MathProbabilityDensityBinomial(const double x,const double n,const double p,int &error_code)
|
|
{
|
|
return MathProbabilityDensityBinomial(x,n,p,false,error_code);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial probability mass function (PDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function calculates the Binomial probability mass function |
|
|
//| with parameters n and p for values in x[] array. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Array with integer random variables |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| log_mode : Logarithm mode flag,if true it calculates Log values|
|
|
//| result : Array with calculated values |
|
|
//| |
|
|
//| Return value: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathProbabilityDensityBinomial(const double &x[],const double n,const double p,const bool log_mode,double &result[])
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
return false;
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
return false;
|
|
//--- check p range
|
|
if(p<0.0 || p>1.0)
|
|
return false;
|
|
|
|
int data_count=ArraySize(x);
|
|
if(data_count==0)
|
|
return false;
|
|
|
|
int error_code=0;
|
|
ArrayResize(result,data_count);
|
|
//--- case p=0 or p=1
|
|
if(p==0.0 || p==1.0)
|
|
{
|
|
for(int i=0; i<data_count; i++)
|
|
result[i]=TailLog0(true,log_mode);
|
|
return true;
|
|
}
|
|
for(int i=0; i<data_count; i++)
|
|
{
|
|
double x_arg=x[i];
|
|
if(MathIsValidNumber(x_arg) && x_arg==MathRound(x_arg))
|
|
{
|
|
//--- check x range
|
|
if(x_arg<0 || x_arg>n)
|
|
result[i]=TailLog0(true,log_mode);
|
|
else
|
|
{
|
|
double log_result=MathGammaLog(n+1.0)-MathGammaLog(x_arg+1.0)-MathGammaLog(n-x_arg+1.0)+x_arg*MathLog(p)+(n-x_arg)*MathLog(1.0-p);
|
|
if(log_mode==true)
|
|
result[i]=log_result;
|
|
else
|
|
result[i]=MathExp(log_result);
|
|
}
|
|
}
|
|
else
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial probability mass function (PDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function calculates the Binomial probability mass function |
|
|
//| with parameters n and p for values in x[] array. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Array with integer random variables |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| result : Array with calculated values |
|
|
//| |
|
|
//| Return value: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathProbabilityDensityBinomial(const double &x[],const double n,const double p,double &result[])
|
|
{
|
|
return MathProbabilityDensityBinomial(x,n,p,false,result);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial cumulative distribution function (CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the Binomial cumulative |
|
|
//| distribution function with given n and p at the desired x. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Integer random variable |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| 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 cumulative distribution function evaluated at x. |
|
|
//+------------------------------------------------------------------+
|
|
double MathCumulativeDistributionBinomial(const double x,const double n,double p,const bool tail,const bool log_mode,int &error_code)
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(x) || !MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
{
|
|
error_code=ERR_ARGUMENTS_NAN;
|
|
return QNaN;
|
|
}
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
//--- check probability
|
|
if(p<0.0 || p>1.0)
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
|
|
error_code=ERR_OK;
|
|
//--- case p==0
|
|
if(p==0.0)
|
|
{
|
|
if(x>=0)
|
|
return TailLog1(tail,log_mode);
|
|
else
|
|
return TailLog0(tail,log_mode);
|
|
}
|
|
//--- case p==1
|
|
if(p==1.0)
|
|
{
|
|
if(x>n)
|
|
return TailLog1(tail,log_mode);
|
|
else
|
|
return TailLog0(tail,log_mode);
|
|
}
|
|
//--- x must be>=0
|
|
if(x<0)
|
|
return TailLog0(tail,log_mode);
|
|
//--- check x
|
|
if(x>n)
|
|
return TailLog1(tail,log_mode);
|
|
int err_code=0;
|
|
//--- calculate using Beta distribution and correct round-off errors
|
|
double result=MathMin(1.0-MathCumulativeDistributionBeta(p,x+1.0,n-x,err_code),1.0);
|
|
//--- return result depending on arguments
|
|
return TailLogValue(result,tail,log_mode);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial cumulative distribution function (CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the Binomial cumulative |
|
|
//| distribution function with given n and p at the desired x. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Integer random variable |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| The cumulative distribution function evaluated at x. |
|
|
//+------------------------------------------------------------------+
|
|
double MathCumulativeDistributionBinomial(const double x,const double n,double p,int &error_code)
|
|
{
|
|
return MathCumulativeDistributionBinomial(x,n,p,true,false,error_code);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial cumulative distribution function (CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the Binomial cumulative |
|
|
//| distribution function with given n and p at the desired x. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Array with integer random variables |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| 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: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathCumulativeDistributionBinomial(const double &x[],const double n,double p,const bool tail,const bool log_mode,double &result[])
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
return false;
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
return false;
|
|
//--- check probability
|
|
if(p<0.0 || p>1.0)
|
|
return false;
|
|
|
|
int data_count=ArraySize(x);
|
|
if(data_count==0)
|
|
return false;
|
|
|
|
int error_code=0;
|
|
ArrayResize(result,data_count);
|
|
|
|
//--- case p=0 and p==1
|
|
if(p==0.0 || p==1.0)
|
|
{
|
|
if(p==0.0)
|
|
{
|
|
for(int i=0; i<data_count; i++)
|
|
{
|
|
if(x[i]>=0)
|
|
result[i]=TailLog1(tail,log_mode);
|
|
else
|
|
result[i]=TailLog0(tail,log_mode);
|
|
}
|
|
}
|
|
else
|
|
//--- p==1.0
|
|
{
|
|
for(int i=0; i<data_count; i++)
|
|
{
|
|
if(x[i]>n)
|
|
result[i]=TailLog1(tail,log_mode);
|
|
else
|
|
result[i]=TailLog0(tail,log_mode);
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
int err_code=0;
|
|
for(int i=0; i<data_count; i++)
|
|
{
|
|
double x_arg=x[i];
|
|
if(MathIsValidNumber(x_arg))
|
|
{
|
|
if(x_arg<0)
|
|
result[i]=TailLog0(tail,log_mode);
|
|
else
|
|
if(x_arg>n)
|
|
result[i]=TailLog1(tail,log_mode);
|
|
else
|
|
{
|
|
double value=MathMin(1.0-MathCumulativeDistributionBeta(p,x_arg+1.0,n-x_arg,err_code),1.0);
|
|
//--- calculate result depending on arguments
|
|
result[i]=TailLogValue(value,tail,log_mode);
|
|
}
|
|
}
|
|
else
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial cumulative distribution function (CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the Binomial cumulative |
|
|
//| distribution function with given n and p for values |
|
|
//| from x[] array. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| x : Array with integer random variables |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathCumulativeDistributionBinomial(const double &x[],const double n,double p,double &result[])
|
|
{
|
|
return MathCumulativeDistributionBinomial(x,n,p,true,false,result);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial distribution quantile function (inverse CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the inverse Binomial cumulative|
|
|
//| distribution function with parameters n and p for the desired |
|
|
//| probability. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| probability : The desired probability |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| tail : Lower tail flag (lower tail of probability used) |
|
|
//| log_mode : Logarithm mode flag (log probability used) |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| The value of the inverse cumulative distribution function |
|
|
//| of the Binomial distribution with parameters n and p. |
|
|
//+------------------------------------------------------------------+
|
|
double MathQuantileBinomial(const double probability,const double n,const double p,const bool tail,const bool log_mode,int &error_code)
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(probability) || !MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
{
|
|
error_code=ERR_ARGUMENTS_NAN;
|
|
return QNaN;
|
|
}
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
//--- check p range
|
|
if(p<0.0 || p>1.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;
|
|
int iterations=0;
|
|
const int max_iterations=1000;
|
|
//--- direct cdf calculation
|
|
double sum=MathProbabilityDensityBinomial(0,n,p,false,error_code);
|
|
while(sum<prob && iterations<max_iterations)
|
|
{
|
|
sum+=MathProbabilityDensityBinomial(iterations,n,p,false,error_code);
|
|
iterations++;
|
|
}
|
|
//--- check convergence
|
|
if(iterations<max_iterations)
|
|
{
|
|
if(iterations==0)
|
|
return 0.0;
|
|
else
|
|
return iterations-1;
|
|
}
|
|
else
|
|
{
|
|
error_code=ERR_RESULT_INFINITE;
|
|
return QPOSINF;
|
|
}
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial distribution quantile function (inverse CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the inverse Binomial cumulative|
|
|
//| distribution function with parameters n and p for the desired |
|
|
//| probability. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| probability : The desired probability |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| The value of the inverse cumulative distribution function |
|
|
//| of the Binomial distribution with parameters n and p. |
|
|
//+------------------------------------------------------------------+
|
|
double MathQuantileBinomial(const double probability,const double n,const double p,int &error_code)
|
|
{
|
|
return MathQuantileBinomial(probability,n,p,true,false,error_code);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial distribution quantile function (inverse CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function calculates the value of the inverse Binomial |
|
|
//| cumulative distribution function with parameters n and p for |
|
|
//| the probability values from probability[] array. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| probability : Array with probability values |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| tail : Lower tail flag (lower tail of probability used) |
|
|
//| log_mode : Logarithm mode flag (log probability used) |
|
|
//| result : Output array with quantile values |
|
|
//| |
|
|
//| Return value: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathQuantileBinomial(const double &probability[],const double n,const double p,const bool tail,const bool log_mode,double &result[])
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
return false;
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
return false;
|
|
//--- check p range
|
|
if(p<0.0 || p>1.0)
|
|
return false;
|
|
|
|
int data_count=ArraySize(probability);
|
|
if(data_count==0)
|
|
return false;
|
|
|
|
int error_code=0;
|
|
ArrayResize(result,data_count);
|
|
|
|
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))
|
|
{
|
|
//--- check probability range
|
|
if(prob<0.0 || prob>1.0)
|
|
return false;
|
|
|
|
int iterations=0;
|
|
//--- direct cdf calculation
|
|
double sum=MathProbabilityDensityBinomial(0,n,p,false,error_code);
|
|
while(sum<prob && iterations<max_iterations)
|
|
{
|
|
sum+=MathProbabilityDensityBinomial(iterations,n,p,false,error_code);
|
|
iterations++;
|
|
}
|
|
//--- check convergence
|
|
if(iterations<max_iterations)
|
|
{
|
|
if(iterations==0)
|
|
result[i]=0;
|
|
else
|
|
result[i]=iterations-1;
|
|
}
|
|
else
|
|
return false;
|
|
}
|
|
else
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial distribution quantile function (inverse CDF) |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns the value of the inverse Binomial cumulative|
|
|
//| distribution function with parameters n and p for the desired |
|
|
//| probability values from probability[] array. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| probability : Array with probability values |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| result : Output array with quantile values |
|
|
//| |
|
|
//| Return value: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathQuantileBinomial(const double &probability[],const double n,const double p,double &result[])
|
|
{
|
|
return MathQuantileBinomial(probability,n,p,true,false,result);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Random variate from Binomial distribution |
|
|
//+------------------------------------------------------------------+
|
|
//| This procedure generates a single random deviate from Binomial |
|
|
//| distribution whose number of trials is N and whose probability |
|
|
//| of an event in each trial is p. |
|
|
//| |
|
|
//| Input parameters: |
|
|
//| n : Number of binomial trials from which a random deviate |
|
|
//| will be generated. |
|
|
//| p : The probability of an event in each trial of the binomial |
|
|
//| distribution from which a random deviate is to be generated. |
|
|
//| |
|
|
//| Return value: |
|
|
//| The random value with Binomial distribution. |
|
|
//| |
|
|
//| Reference: |
|
|
//| Voratas Kachitvichyanukul, Bruce Schmeiser, |
|
|
//| "Binomial Random Variate Generation", Communications of the ACM, |
|
|
//| Volume 31, Number 2, February 1988, pages 216-222. |
|
|
//| |
|
|
//| Original FORTRAN77 version by Barry Brown, James Lovato. |
|
|
//| C version by John Burkardt. |
|
|
//+------------------------------------------------------------------+
|
|
double MathRandomBinomial(const double n,const double p)
|
|
{
|
|
int ix,ix1,mp;
|
|
double f,g,qn,r,t,u,v,w,w2,x,z;
|
|
int value=0;
|
|
int n1=(int)n;
|
|
|
|
double pp= MathMin(p,1.0-p);
|
|
double q = 1.0-pp;
|
|
double xnp=(double)(n1)*pp;
|
|
|
|
if(xnp<30.0)
|
|
{
|
|
qn= MathPow(q,n1);
|
|
r = pp/q;
|
|
g = r*(double)(n1+1);
|
|
|
|
for(;;)
|
|
{
|
|
ix= 0;
|
|
f = qn;
|
|
u = MathRandomNonZero();
|
|
|
|
for(;;)
|
|
{
|
|
if(u<f)
|
|
{
|
|
if(0.5<p)
|
|
{
|
|
ix=n1-ix;
|
|
}
|
|
value=ix;
|
|
return value;
|
|
}
|
|
if(110<ix)
|
|
break;
|
|
u=u-f;
|
|
ix=ix+1;
|
|
f=f*(g/(double)(ix)-r);
|
|
}
|
|
}
|
|
}
|
|
double ffm=xnp+pp;
|
|
int m=int(ffm);
|
|
double fm=m;
|
|
double xnpq=xnp*q;
|
|
double p1 = (int)(2.195*MathSqrt(xnpq)-4.6*q)+0.5;
|
|
double xm = fm + 0.5;
|
|
double xl = xm - p1;
|
|
double xr = xm + p1;
|
|
double c=0.134+20.5/(15.3+fm);
|
|
double al=(ffm-xl)/(ffm-xl*pp);
|
|
double xll=al*(1.0+0.5*al);
|
|
al=(xr-ffm)/(xr*q);
|
|
double xlr= al*(1.0 + 0.5*al);
|
|
double p2 = p1*(1.0 + c + c);
|
|
double p3 = p2 + c/xll;
|
|
double p4 = p3 + c/xlr;
|
|
//--- generate a variate
|
|
for(;;)
|
|
{
|
|
u = MathRandomNonZero()*p4;
|
|
v = MathRandomNonZero();
|
|
//--- triangle
|
|
if(u<p1)
|
|
{
|
|
ix=int(xm-p1*v+u);
|
|
if(0.5<p)
|
|
ix=n1-ix;
|
|
|
|
value=ix;
|
|
return value;
|
|
}
|
|
//--- parallelogram
|
|
if(u<=p2)
|
|
{
|
|
x = xl+(u - p1)/c;
|
|
v = v*c + 1.0 - MathAbs(xm-x)/p1;
|
|
|
|
if(v<=0.0 || 1.0<v)
|
|
continue;
|
|
ix=int(x);
|
|
}
|
|
else
|
|
if(u<=p3)
|
|
{
|
|
ix=int(xl+MathLog(v)/xll);
|
|
if(ix<0)
|
|
continue;
|
|
v=v*(u-p2)*xll;
|
|
}
|
|
else
|
|
{
|
|
ix=int(xr-MathLog(v)/xlr);
|
|
if(n1<ix)
|
|
continue;
|
|
v=v*(u-p3)*xlr;
|
|
}
|
|
int k=MathAbs(ix-m);
|
|
|
|
if(k<=20 || xnpq/2.0-1.0<=k)
|
|
{
|
|
f = 1.0;
|
|
r = pp/q;
|
|
g = (n1+1)*r;
|
|
|
|
if(m<ix)
|
|
{
|
|
mp=m+1;
|
|
for(int i=mp; i<=ix; i++)
|
|
f=f*(g/i-r);
|
|
}
|
|
else
|
|
if(ix<m)
|
|
{
|
|
ix1=ix+1;
|
|
for(int i=ix1; i<=m; i++)
|
|
f=f/(g/i-r);
|
|
}
|
|
|
|
if(v<=f)
|
|
{
|
|
if(0.5<p)
|
|
ix=n1-ix;
|
|
|
|
value=ix;
|
|
return value;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
double amaxp=(k/xnpq)*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5);
|
|
double ynorm=-double((k*k)/(2.0*xnpq));
|
|
double alv=MathLog(v);
|
|
|
|
if(alv<ynorm-amaxp)
|
|
{
|
|
if(0.5<p)
|
|
ix=n1-ix;
|
|
|
|
value=ix;
|
|
return value;
|
|
}
|
|
|
|
if(ynorm+amaxp<alv)
|
|
continue;
|
|
|
|
double x1 = double(ix+1);
|
|
double f1 = fm + 1.0;
|
|
z = (double)(n1+1) - fm;
|
|
w = (double)(n1-ix+1);
|
|
double z2 = z * z;
|
|
double x2 = x1 * x1;
|
|
double f2 = f1 * f1;
|
|
w2=w*w;
|
|
|
|
t=xm*MathLog(f1/x1)+(n1-m+0.5)*MathLog(z/w)+(double)(ix-m)*MathLog(w*pp/(x1*q))
|
|
+(13860.0 -(462.0 -(132.0 -(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0
|
|
+(13860.0 -(462.0 -(132.0 -(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0
|
|
+(13860.0 -(462.0 -(132.0 -(99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0
|
|
+(13860.0 -(462.0 -(132.0 -(99.0-140.0/w2)/w2)/w2)/w2)/w/166320.0;
|
|
|
|
if(alv<=t)
|
|
{
|
|
if(0.5<p)
|
|
ix=n1-ix;
|
|
|
|
value=ix;
|
|
return value;
|
|
}
|
|
}
|
|
}
|
|
return value;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Random variate from Binomial distribution |
|
|
//+------------------------------------------------------------------+
|
|
//| The function returns random deviate from Binomial distribution |
|
|
//| with parameters n and p. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| error_code : Variable for error code |
|
|
//| |
|
|
//| Return value: |
|
|
//| The random value with Binomial distribution. |
|
|
//+------------------------------------------------------------------+
|
|
double MathRandomBinomial(const double n,const double p,int &error_code)
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
{
|
|
error_code=ERR_ARGUMENTS_NAN;
|
|
return QNaN;
|
|
}
|
|
//--- check n
|
|
if(n<=0 || n!=MathRound(n))
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
//--- check probability
|
|
if(p<=0 || p>=1.0)
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return QNaN;
|
|
}
|
|
|
|
error_code=ERR_OK;
|
|
//--- return binomial random value
|
|
return MathRandomBinomial(n,p);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Random variate from Binomial distribution |
|
|
//+------------------------------------------------------------------+
|
|
//| The function generates random variables from Binomial |
|
|
//| distribution with parameters n and p. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| data_count : Number of values needed |
|
|
//| result : Output array with random values |
|
|
//| |
|
|
//| Return value: |
|
|
//| true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool MathRandomBinomial(const double n,const double p,const int data_count,double &result[])
|
|
{
|
|
if(data_count<=0)
|
|
return false;
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(n) || !MathIsValidNumber(p))
|
|
return false;
|
|
//--- check n
|
|
if(n<=0 || n!=MathRound(n))
|
|
return false;
|
|
//--- check probability
|
|
if(p<=0 || p>=1.0)
|
|
return false;
|
|
//--- prepare output array and calculate random values
|
|
ArrayResize(result,data_count);
|
|
for(int i=0; i<data_count; i++)
|
|
result[i]=MathRandomBinomial(n,p);
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Binomial distriburion moments |
|
|
//+------------------------------------------------------------------+
|
|
//| The function calculates 4 first moments of the Binomial |
|
|
//| distribution with parameters n and p. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| n : Number of trials |
|
|
//| p : Probability of success for each trial |
|
|
//| 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 MathMomentsBinomial(const double n,const double p,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(n) || !MathIsValidNumber(p))
|
|
{
|
|
error_code=ERR_ARGUMENTS_NAN;
|
|
return false;
|
|
}
|
|
//--- check n
|
|
if(n<0 || n!=MathRound(n))
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return false;
|
|
}
|
|
//--- check p range
|
|
if(p<=0.0 || p>=1.0)
|
|
{
|
|
error_code=ERR_ARGUMENTS_INVALID;
|
|
return false;
|
|
}
|
|
|
|
error_code=ERR_OK;
|
|
//--- prepare factors
|
|
double np=n*p;
|
|
double one_mp=(1.0-p);
|
|
//--- calculate moments
|
|
mean =np;
|
|
variance=np*one_mp;
|
|
skewness=(1-2*p)/MathSqrt(variance);
|
|
kurtosis=(1-6*p*one_mp)/variance;
|
|
//--- successful
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|