//+------------------------------------------------------------------+ //| Math.mqh | //| Copyright 2000-2025, MetaQuotes Ltd. | //| https://www.mql5.com | //+------------------------------------------------------------------+ #define ERR_OK 0 #define ERR_ARGUMENTS_NAN 1 #define ERR_ARGUMENTS_INVALID 2 #define ERR_RESULT_INFINITE 3 #define ERR_NON_CONVERGENCE 4 //+------------------------------------------------------------------+ //| Nans | //+------------------------------------------------------------------+ double QNaN =(double)"nan"; // QNaN double QPOSINF=(double)"inf"; // +infinity 1.#INF double QNEGINF=(double)"-inf"; // -infinity -1.#INF //+------------------------------------------------------------------+ //| MathRandom | //+------------------------------------------------------------------+ double MathRandomNonZero(void) { double rnd=0; //--- except 0 and 1 while(rnd==0.0 || rnd==1.0) rnd=MathRand()/32767.0; //--- return(rnd); } //+------------------------------------------------------------------+ //| Computes 4 first moments of the values in array[] | //+------------------------------------------------------------------+ bool MathMoments(const double &array[],double &mean,double &variance,double &skewness,double &kurtosis,const int start=0,const int count=WHOLE_ARRAY) { //--- initial values mean=0.0; variance=0.0; skewness=0.0; kurtosis=0.0; //--- get data size int size=ArraySize(array); int data_count=0; //--- set data count if(count==WHOLE_ARRAY) data_count=size; else data_count=count; //--- check data range if(data_count==0) return(false); if(start+data_count>size) return(false); //--- set indexes int ind1=start; int ind2=ind1+data_count-1; //--- calculate mean for(int i=ind1; i<=ind2; i++) mean+=array[i]; mean=mean/data_count; double S=0.0; //--- if(data_count>1) { //--- calculate variance for(int i=ind1; i<=ind2; i++) variance+=MathPow(array[i]-mean,2); variance=variance/(data_count-1); S=MathSqrt(variance); } else variance=QNaN; //--- if(data_count>2 && S>0.0) { //--- calculate skewness for(int i=ind1; i<=ind2; i++) skewness+=MathPow(array[i]-mean,3); skewness=skewness/(data_count*MathPow(S,3)); } else skewness=QNaN; //--- if(data_count>3 && S>0.0) { //--- calculate kurtosis for(int i=ind1; i<=ind2; i++) kurtosis+=MathPow(array[i]-mean,4); kurtosis=kurtosis/(data_count*MathPow(S,4)); kurtosis-=3; } else kurtosis=QNaN; //--- check values and return calculation result if((!MathIsValidNumber(variance)) || (!MathIsValidNumber(skewness)) || (!MathIsValidNumber(kurtosis))) return(false); else return(true); } #define M_1_SQRT_2PI 1/MathSqrt(2*M_PI) double FactorialsTable[21]= { 1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600, 6227020800,87178291200,1307674368000,20922789888000,355687428096000, 6402373705728000,121645100408832000,2432902008176640000 }; //+------------------------------------------------------------------+ //| MathPowInt | //+------------------------------------------------------------------+ double MathPowInt(const double x,const int power) { //--- check power if(power==0) return(1.0); if(power<0) return(0); //--- calculate product double val=x; for(int i=2; i<=power; i++) val*=x; //--- return result return(val); } //+------------------------------------------------------------------+ //| MathFactorial | //+------------------------------------------------------------------+ double MathFactorial(const int n) { if(n<0) return(0); if(n<=20) //--- use values from the factorials table return(FactorialsTable[n]); else { //--- calculate product starting from 20th element of factorials table double val=FactorialsTable[20]; for(int i=21; i<=n; i++) val*=i; //--- return result return(val); } } //+------------------------------------------------------------------+ //| MathTrunc | //+------------------------------------------------------------------+ double MathTrunc(const double x) { if(x>=0) return(MathFloor(x)); else return(MathCeil(x)); } //+------------------------------------------------------------------+ //| MathRound | //+------------------------------------------------------------------+ //| Round a double number to a given precision | //+------------------------------------------------------------------+ double MathRound(const double x,const int digits) { if(!MathIsValidNumber(x)) return(QNaN); if(x==0.0) return(x); if(digits>0) { double sign=1.0; double xx=x; if(xx<0.0) { xx=-xx; sign=-1.0; } double pwr_factor=MathPow(10,digits); double int_value=MathFloor(xx); double res=(xx-int_value)*pwr_factor; res=MathFloor(MathAbs(res)+0.5); return(sign*(int_value+res/pwr_factor)); } else if(digits==0) return MathRound(x); //--- return(0); } //+------------------------------------------------------------------+ //| Comments from original FORTRAN code: | //| https://www.netlib.org/specfun/gamma | //| | //| This routine calculates the GAMMA function for a real argument X.| //| Computation is based on an algorithm outlined in reference 1. | //| The program uses rational functions that approximate the GAMMA | //| function to at least 20 significant decimal digits. Coefficients | //| for the approximation over the interval (1,2) are unpublished. | //| Those for the approximation for X .GE. 12 are from reference 2. | //| | //| References: | //| 1. "An Overview of Software Development for Special Functions" | //| W. J. Cody, Lecture Notes in Mathematics, 506, | //| Numerical Analysis Dundee, 1975, G. A. Watson (ed.) | //| Springer Verlag, Berlin, 1976. | //| 2. Computer Approximations, Hart, Et. Al., Wiley and sons, | //| New York, 1968. | //| | //| Authors: | //| W. J. Cody and L. Stoltz, Applied Mathematics Division | //| Argonne National Laboratory, Argonne, IL 60439 | //+------------------------------------------------------------------+ double MathGamma(const double x) { //--- mathematical constants const double logsqrt2pi=MathLog(MathSqrt(2*M_PI)); //--- machine dependent parameters const double xminin=2.23e-308; const double xbig=171.624; double eps=2.22e-16; //--- numerator and denominator coefficients for rational minimax //--- approximation over (1,2) const double p[8]= { -1.71618513886549492533811e+0, 2.47656508055759199108314e+1, -3.79804256470945635097577e+2, 6.29331155312818442661052e+2, 8.66966202790413211295064e+2,-3.14512729688483675254357e+4, -3.61444134186911729807069e+4, 6.64561438202405440627855e+4 }; const double q[8]= { -3.08402300119738975254353e+1, 3.15350626979604161529144e+2, -1.01515636749021914166146e+3,-3.10777167157231109440444e+3, 2.25381184209801510330112e+4, 4.75584627752788110767815e+3, -1.34659959864969306392456e+5,-1.15132259675553483497211e+5 }; //--- coefficients for minimax approximation over (12, inf) const double c[7]= { -1.910444077728e-03, 8.4171387781295e-04, -5.952379913043012e-04, 7.93650793500350248e-04, -2.777777777777681622553e-03,8.333333333333333331554247e-02, 5.7083835261e-03 }; int parity=0; double fact=1.0; int n=0; double y=x,y1,res,z,xnum=0,xden=0,ysq=0,sum=0; if(y<=0.0) { //--- argument is negative y=-x; y1=MathTrunc(y); res=y-y1; if(res!=0.0) { if(y1!=MathTrunc(y1*0.5)*2.0) parity=1; fact=-M_PI/MathSin(M_PI*res); y=y+1.0; } else { return(QPOSINF); } } //--- argument is positive if(y=xminin) { res=1.0/y; } else { return(QPOSINF); } } else if(y<12.0) { y1=y; if(y<1.0) { //--- 0.0 < argument < 1.0 z = y; y = y + 1.0; } else { //--- 1.0 < argument < 12.0, reduce argument if necessary n = int(y) - 1; y = y - double(n); z = y - 1.0; } //--- evaluate approximation for 1.0 < argument < 2.0 xnum=0.0; xden=1.0; for(int i=0; i<8; i++) { xnum = (xnum + p[i]) * z; xden = xden * z + q[i]; } res=xnum/xden+1.0; if(y1y) { //--- adjust result for case 2.0 < argument < 12.0 for(int i=0; i 12 is from reference | //| 3, while approximations for X < 12.0 are similar to those in | //| reference 1, but are unpublished. The accuracy achieved depends | //| on the arithmetic system, the compiler, the intrinsic functions, | //| and proper selection of the machine-dependent constants. | //| | //| References: | //| 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for | //| the Natural Logarithm of the Gamma Function,' Math. Comp. 21, | //| 1967, pp. 198-203. | //| 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,| //| 1969. | //| 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New | //| York, 1968. | //| Authors: W. J. Cody and L. Stoltz | //| Argonne National Laboratory | //| | //| https://gcc.gnu.org/ml/fortran/2007-11/msg00061/gamma.diff | //| For negative arguments (where netlib would return +Inf) | //| we use the following identity: | //| lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) | //+------------------------------------------------------------------+ double MathGammaLog(const double x) { //--- mathematical constants const double pnt68=0.6796875e0; const double sqrtpi=0.9189385332046727417803297e0; //--- machine dependent parameters const double xbig=2.55E305; const double xinf=1.79E308; const double eps=2.22E-16; const double frtbig=2.25E76; //--- numerator and denominator coefficients for rational minimax //--- approximation over (0.5,1.5) const double d1=-5.772156649015328605195174e-1; const double p1[8]= { 4.945235359296727046734888e0,2.018112620856775083915565e2, 2.290838373831346393026739e3,1.131967205903380828685045e4, 2.855724635671635335736389e4,3.848496228443793359990269e4, 2.637748787624195437963534e4,7.225813979700288197698961e3 }; const double q1[8]= { 6.748212550303777196073036e1,1.113332393857199323513008e3, 7.738757056935398733233834e3,2.763987074403340708898585e4, 5.499310206226157329794414e4,6.161122180066002127833352e4, 3.635127591501940507276287e4,8.785536302431013170870835e3 }; //--- numerator and denominator coefficients for rational minimax //--- approximation over (1.5,4.0) const double d2=4.227843350984671393993777e-1; const double p2[8]= { 4.974607845568932035012064e0,5.424138599891070494101986e2, 1.550693864978364947665077e4,1.847932904445632425417223e5, 1.088204769468828767498470e6,3.338152967987029735917223e6, 5.106661678927352456275255e6,3.074109054850539556250927e6 }; const double q2[8]= { 1.830328399370592604055942e2,7.765049321445005871323047e3, 1.331903827966074194402448e5,1.136705821321969608938755e6, 5.267964117437946917577538e6,1.346701454311101692290052e7, 1.782736530353274213975932e7,9.533095591844353613395747e6 }; //--- numerator and denominator coefficients for rational minimax //--- approximation over (4.0,12.0) const double d4=1.791759469228055000094023e0; const double p4[8]= { 1.474502166059939948905062e4, 2.426813369486704502836312e6, 1.214755574045093227939592e8, 2.663432449630976949898078e9, 2.940378956634553899906876e10,1.702665737765398868392998e11, 4.926125793377430887588120e11,5.606251856223951465078242e11 }; const double q4[8]= { 2.690530175870899333379843e3, 6.393885654300092398984238e5, 4.135599930241388052042842e7, 1.120872109616147941376570e9, 1.488613728678813811542398e10,1.016803586272438228077304e11, 3.417476345507377132798597e11,4.463158187419713286462081e11 }; //--- coefficients for minimax approximation over (12, inf) const double c[7]= { -1.910444077728e-03, 8.4171387781295e-04, -5.952379913043012e-04, 7.93650793500350248e-04, -2.777777777777681622553e-03,8.333333333333333331554247e-02, 5.7083835261e-03 }; double y=x; double res=0.0; double corr=0; double xm1,xm2,xm4,xden,xnum; double ysq=0; if((y>0 && y<=xbig)) { if(y<=eps) { res=-MathLog(y); } else { if(y<=1.5) { //--- eps < x <= 1.5 if(y=pnt68)) { xden = 1.0; xnum = 0; for(int i=0; i<8; i++) { xnum = xnum*xm1 + p1[i]; xden = xden*xm1 + q1[i]; } res=corr+(xm1 *(d1+xm1*(xnum/xden))); } else { xm2=(y-0.5)-0.5; xden = 1.0; xnum = 0.0; for(int i=0; i<8; i++) { xnum=xnum*xm2+p2[i]; xden=xden*xm2+q2[i]; } res=corr+xm2 *(d2+xm2*(xnum/xden)); } } else { if(y<=4.0) { //--- 1.5 < x .<= 4.0 xm2=y-2.0; xden = 1.0; xnum = 0.0; for(int i=0; i<8; i++) { xnum = xnum*xm2 + p2[i]; xden = xden*xm2 + q2[i]; } res=xm2 *(d2+xm2*(xnum/xden)); } else { if(y<=12.0) { //--- 4.0 < x <= 12.0 xm4=y-4.0; xden = -1.0; xnum = 0.0; for(int i=0; i<8; i++) { xnum = xnum*xm4 + p4[i]; xden = xden*xm4 + q4[i]; } res=d4+xm4*(xnum/xden); } else { //--- evaluate for argument>=12.0, res=0.0; if(y<=frtbig) { res=c[6]; ysq=y*y; for(int i=0; i<6; i++) { res=res/ysq+c[i]; } } res=res/y; corr= MathLog(y); res = res + sqrtpi - 0.5*corr; res = res + y*(corr-1.0); } } } } } else { //--- return for bad arguments //--- res=QPOSINF; //--- for negative arguments (where netlib would return +Inf) //--- we use the following identity: //--- lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) //--- https://gcc.gnu.org/ml/fortran/2007-11/msg00061/gamma.diff if(y<0 && MathFloor(y)!=y) { //--- for abs(y) very close to zero, we use a series expansion to //--- the first order in y to avoid overflow. if(y>-1.e-100) res=-2*MathLog(MathAbs(y))-MathGammaLog(-y); else res=MathLog(M_PI/MathAbs(y*MathSin(M_PI*y)))-MathGammaLog(-y); } else { //--- negative integer values res=QPOSINF; } } //--- return result return(res); } //+------------------------------------------------------------------+ //| MathBeta | //+------------------------------------------------------------------+ double MathBeta(const double a,const double b) { return MathExp(MathBetaLog(a,b)); } //+------------------------------------------------------------------+ //| MathBetaLog | //+------------------------------------------------------------------+ double MathBetaLog(const double a,const double b) { return(MathGammaLog(a)+MathGammaLog(b)-MathGammaLog(a+b)); } //+------------------------------------------------------------------+ //| Incomplete beta function | //+------------------------------------------------------------------+ //| The incomplete beta function ratio is the probability that a | //| random variable from a beta distribution having parameters p | //| and q will be less than or equal to x. | //| | //| Input: | //| x : upper limit of integration. x must be in (0,1) inclusive. | //| p : first beta distribution parameter. p must be>0.0. | //| q : second beta distribution parameter. q must be>0.0. | //| | //| Returns the value of the incomplete Beta function ratio. | //| | //| Original FORTRAN77 version by KL Majumder, GP Bhattacharjee. | //| C version by John Burkardt. | //| | //| Reference: | //| | //| KL Majumder, GP Bhattacharjee, | //| Algorithm AS 63: The incomplete Beta Integral, | //| Applied Statistics,.Volume 22, Number 3, 1973, pages 409-411. | //+------------------------------------------------------------------+ double MathBetaIncomplete(const double x,const double p,const double q) { double acu=0.1E-16; double pp,qq,xx; int ifault=0; double value=x; //--- check the input arguments if(p<=0.0 || q<=0.0) { ifault=1; return(value); } //--- check x range if(x<0.0 || 1.00) | //+------------------------------------------------------------------+ //| Comment from original FORTRAN code: | //| www.nist.gov/sites/default/files/documents/itl/sed/stspac.f | //| ibid .../itl/sed/SED_Note_86-2-2.pdf | //| | //| CDFGAM Written by Charles p. Reeve, Statistical Engineering | //| Division, National Bureau of Standards, Gaithersburg. | //| | //| Computing the cumulative distribution function of the Gamma | //| distribution (also known as the Incomplete Gamma ratio) to a | //| specified accuracy (truncation error in the infinite series). | //| the algorithm, described in Ref. 2, is a modification of | //| the algorithm of Ref. 1. Three features have been added: | //| | //| 1) A precise method of meeting the truncation accuracy, | //| 2) Computation of the upper tail area by decrementing alpha | //| when that method is more efficient, and | //| 3) A constant uflo >= the underflow limit on the computer. | //| | //| References: | //| | //| 1) Lau, Chi-Leung, 'A Simple Series for the Incomplete Gamma | //| Integral', Algorithm AS 147, Applied Statistics, vol. 29, | //| No. 1, 1980, pp. 113-114. | //| | //| 2) Reeve, Charles p., 'An Algorithm for Computing the Gamma | //| C.D.F. to a Specified Accuracy', Statistical Engineering | //| Division,Note 86-2, October 1986. | //+------------------------------------------------------------------+ double MathGammaIncomplete(double x,double alpha) { double eps=10e-20; int iflag=0; bool ll; int imax=5000; double uflo=1.0e-100; double cdfx=0.0; int k=0; //--- check for validity of arguments alpha and eps if(alpha<=uflo || eps<=uflo) { iflag=1; return(QNaN); } iflag=0; //--- check for special case of x if(x<=0.0) return(QNaN); //--- evaluate the gamma p.d.f. and check for underflow double dx=double(x); double pdfl=double(alpha-1.0)*MathLog(dx)-dx-MathGammaLog(alpha); if(pdfl=alpha) cdfx=1.0; } else { double p=alpha; double u=MathExp(pdfl); //--- determine whether to increment or decrement alpha (a.k.a. p) ll=true; if(x>=p) { k=int(p); if(p<=double(k)) k=k-1; double eta=p-double(k); double bl=double((eta-1.0)*MathLog(dx)-dx-MathGammaLog(eta)); ll=bl>MathLog(eps); } //--- double epsx=eps/x; //--- if(ll==true) { //--- increment p for(int i=0; i<=imax; i++) { if(u<=epsx*(p-x)) { return(cdfx); } u=x*u/p; cdfx=cdfx+u; p=p+1.0; } iflag=2; } else { //--- decrement p for(int j=1; j<=k; j++) { p=p-1.0; if(u<=epsx*(x-p)) break; cdfx=cdfx+u; u=p*u/x; } cdfx=1.0-cdfx; } } //--- return result return(cdfx); } //+------------------------------------------------------------------+ //| Computes the binomial coefficient C(n,k)=n!/(k!*(n-k)!) | //| //| The value is calculated in such a way as to avoid overflow and | //| roundoff. The calculation is done in integer arithmetic. | //| | //| Parameters: | //| Input, int N, K, are the values of N and K. | //| Output: the number of combinations of N things taken K at a time.| //| | //| Author: John Burkardt | //| | //| Reference: | //| ML Wolfson, HV Wright, | //| Algorithm 160: Combinatorial of M Things Taken N at a Time, | //| Communications of the ACM, Volume 6, N4, April 1963, page 161. | //+------------------------------------------------------------------+ long MathBinomialCoefficient(const int n,const int k) { int mn,mx; long value=0; int n_k=n-k; //--- if(k0) { value=mx+1; for(int i=2; i<=mn; i++) value=(value*(mx+i))/i; } else { if(mn<0) return(0); if(mn==0) return(1); } //--- return(value); } //+------------------------------------------------------------------+ //| MathBinomialCoefficientLog | //+------------------------------------------------------------------+ double MathBinomialCoefficientLog(const int n,const int k) { int mn,mx; long value=0; int n_k=n-k; //--- if(k0) { value=mx+1; for(int i=2; i<=mn; i++) value=(value*(mx+i))/i; } else { if(mn<0) return(QNEGINF); if(mn==0) return(0); } //--- return MathLog(value); } //+------------------------------------------------------------------+ //| Computes the logarithm of the Binomial coefficient. | //| | //| Log(C(n,k))=Log(n!/(k!*(n-k)!)) | //| | //| Parameters: | //| Input, int n, k, are the values of n and k. | //| Return value: the logarithm of C(n,k). | //+------------------------------------------------------------------+ double MathBinomialCoefficientLog(const double n,const double k) { return(MathGammaLog(n+1)-MathGammaLog(k+1)-MathGammaLog(n-k+1)); } //+------------------------------------------------------------------+ //| Computes the function Hypergeometric_2F2(a,b,c,d,z) using | //| a Taylor method. | //| | //| Author : John Pearson | //| Reference : MSc thesis "Computation of Hypergeometric Functions" | //+------------------------------------------------------------------+ double MathHypergeometric2F2(const double a,const double b,const double c,const double d,const double z) { double a1=1.0; double b1=1.0; double tol=10E-10; // set tolerance //--- direct summation for(int j=1; j<=500; j++) { //--- update current term in terms of previous one a1=(a+j-1)*(b+j-1)/(c+j-1)/(d+j-1)*z/j*a1; //--- update sum of terms computed so far b1=b1+a1; //--- terminate summation if stopping criterion is satisfied if(MathAbs(a1)/MathAbs(b1)>1" is quick division by 2 p_double=array[(first+last)>>1]; while(ip_double) { //--- control the output of the array bounds if(j==0) break; j--; } if(i<=j) { //-- swap elements i and j t_double=array[i]; array[i]=array[j]; array[j]=t_double; //-- swap indices i and j t_int=indices[i]; indices[i]=indices[j]; indices[j]=t_int; i++; //--- control the output of the array bounds if(j==0) break; j--; } } if(first>1" is quick division by 2 p_double=array[(first+last)>>1]; while(ip_double) { //--- control the output of the array bounds if(i==ArraySize(array)-1) break; i++; } while(array[j]0 ascending,otherwise descending) | //| | //| Return value: None | //+------------------------------------------------------------------+ void MathQuickSort(double &array[],int &indices[],int first,int last,int mode) { if(mode>0) MathQuickSortAscending(array,indices,first,last); else MathQuickSortDescending(array,indices,first,last); } //+------------------------------------------------------------------+ //| MathOrder | //+------------------------------------------------------------------+ //| The function calculates the order of elements from the array. | //| | //| Arguments: | //| array[] : Array with values | //| result[] : Array for sorted indices | //| | //| Return value: true if successful, otherwise false. | //+------------------------------------------------------------------+ bool MathOrder(const double &array[],int &result[]) { int size=ArraySize(array); //--- check array size if(size==0) return(false); //--- prepare temporary array double tmp_array[]; ArrayCopy(tmp_array,array,0,0,WHOLE_ARRAY); //--- prepare identical permutation if(ArraySize(result)> logical operation for elements | //| of the array[]. | //| | //| Arguments: | //| array[] : Array with values | //| n : Number of bits to shift | //| result[] : Array for sorted indices | //| | //| Return value: true if successful, otherwise false. | //+------------------------------------------------------------------+ bool MathBitwiseShiftR(const int &array[],const int n,int &result[]) { int size=ArraySize(array); //--- if(size==0) return(false); //--- prepare target array and calculate array>>n if(ArraySize(result)>n; //--- return(true); } //+------------------------------------------------------------------+ //| MathBitwiseShiftR | //+------------------------------------------------------------------+ //| The function calculates >> logical operation for elements | //| of the array[]. | //| | //| Arguments: | //| array[] : Array with values | //| n : Number of bits to shift | //| | //| Return value: true if successful, otherwise false. | //+------------------------------------------------------------------+ bool MathBitwiseShiftR(int &array[],const int n) { int size=ArraySize(array); if(size==0) return(false); //--- calculate array>>n for(int i=0; i>n; //--- return(true); } //+------------------------------------------------------------------+ //| MathCumulativeSum | //+------------------------------------------------------------------+ //| The function calculates cumulative sum of the elements | //| from the array[]. | //| | //| Arguments: | //| array[] : Array with values | //| result[] : Output array with cumulative sum | //| | //| Return value: true if successful, otherwise false | //+------------------------------------------------------------------+ bool MathCumulativeSum(const double &array[],double &result[]) { int size=ArraySize(array); if(size==0) return(false); //--- prepare target array and calculate cumulative sum if(ArraySize(result)=0) result[i]=MathFloor(array[i]); else result[i]=MathCeil(array[i]); } //--- return(true); } //+------------------------------------------------------------------+ //| MathTrunc | //+------------------------------------------------------------------+ //| The function calculates trunc(x) for the elements from the array.| //| | //| Arguments: | //| array[] : Array with values | //| | //| Return value: true if successful, otherwise false | //+------------------------------------------------------------------+ bool MathTrunc(double &array[]) { int size=ArraySize(array); if(size==0) return(false); for(int i=0; i=0) array[i]=MathFloor(array[i]); else array[i]=MathCeil(array[i]); } //--- return(true); } //+------------------------------------------------------------------+ //| MathSqrt | //+------------------------------------------------------------------+ //| The function calculates sqrt(x) for the elements from the array. | //| | //| Arguments: | //| array[] : Array with values | //| result[] : Output array with calculated values | //| | //| Return value: true if successful, otherwise false | //+------------------------------------------------------------------+ bool MathSqrt(const double &array[],double &result[]) { int size=ArraySize(array); if(size==0) return(false); //--- prepare target array if(ArraySize(result)=size) return(false); //--- copy data from initial array if(ArrayCopy(result,array,0,0,WHOLE_ARRAY)!=size) return(false); //--- calculate iterated differences for(int i=0; i=size) return(false); //--- copy data from initial array if(ArrayCopy(result,array,0,0,WHOLE_ARRAY)!=size) return(false); //--- calculate iterated differences for(int i=0; isize) return(false); //--- prepare target array if(ArraySize(result)size) return(false); //--- prepare target array if(ArraySize(result)total samples if(count<=0 || count>size) return(false); //--- probabilities double prob[]; if(ArrayCopy(prob,probabilities,0,0,WHOLE_ARRAY)!=size) return(false); //--- calculate sum and normalize probabilities double sum=0; for(int i=0; itotal samples if(count<=0 || count>size) return(false); //--- probabilities double prob[]; if(ArrayCopy(prob,probabilities,0,0,WHOLE_ARRAY)!=size) return(false); //--- calculate sum and normalize probabilities double sum=0; for(int i=0; i30) return x; if(dig<1) dig=1; double sign=1.0; double xx=x; if(xx<0.0) { sign=-sign; xx=-xx; } //--- calculate log double l10 = MathLog10(xx); double e10 = (int)(dig-1-MathFloor(l10)); double value=0,pwr10; if(e10>0) { pwr10=MathPow(10,e10); value=MathFloor((xx*pwr10)+0.5)/pwr10; } else { pwr10=MathPow(10,-e10); value=MathFloor((xx/pwr10)+0.5)*pwr10; } return(sign*value); } //+------------------------------------------------------------------+ //| MathSignif | //+------------------------------------------------------------------+ //| The function calculates the rounded values (to the specified | //| number of significant digits) for the elements from the array[]. | //| | //| Arguments: | //| array[] : Array with double values | //| digits : Number of significant digits | //| result[] : Output array with calculated values | //| | //| Return value: true if successful, otherwise false | //+------------------------------------------------------------------+ bool MathSignif(const double &array[],int digits,double &result[]) { int size=ArraySize(array); if(size==0) return(false); //--- prepare target array if(ArraySize(result)=values[t-1]) t=1; else { //--- swap tmp=values[k-1]; values[k-1]=values[t-1]; values[t-1]=tmp; tmpi=indices[k-1]; indices[k-1]=indices[t-1]; indices[t-1]=tmpi; t=k; } } i=i+1; } while(i<=size); i=size-1; do { //--- swap tmp=values[i]; values[i]=values[0]; values[0]=tmp; tmpi=indices[i]; indices[i]=indices[0]; indices[0]=tmpi; t=1; while(t!=0) { k=2*t; if(k>i) t=0; else { if(kvalues[k-1]) k++; if(values[t-1]>=values[k-1]) t=0; else { //--- swap tmp=values[k-1]; values[k-1]=values[t-1]; values[t-1]=tmp; tmpi=indices[k-1]; indices[k-1]=indices[t-1]; indices[t-1]=tmpi; t=k; } } } i=i-1; } while(i>=1); } //--- compute tied ranks i=0; while(i=values[t-1]) t=1; else { //--- swap tmp=values[k-1]; values[k-1]=values[t-1]; values[t-1]=tmp; tmpi=indices[k-1]; indices[k-1]=indices[t-1]; indices[t-1]=tmpi; t=k; } } i=i+1; } while(i<=size); i=size-1; do { //--- swap tmp=values[i]; values[i]=values[0]; values[0]=tmp; tmpi=indices[i]; indices[i]=indices[0]; indices[0]=tmpi; t=1; while(t!=0) { k=2*t; if(k>i) t=0; else { if(kvalues[k-1]) k++; if(values[t-1]>=values[k-1]) t=0; else { //--- swap tmp=values[k-1]; values[k-1]=values[t-1]; values[t-1]=tmp; tmpi=indices[k-1]; indices[k-1]=indices[t-1]; indices[t-1]=tmpi; t=k; } } } i=i-1; } while(i>=1); } //--- compute tied ranks i=0; while(i0.0) cnt++; else cnt--; } } } //--- calculate Kendall tau double den=cnt1*cnt2; if(den==0) return(false); tau=cnt/MathSqrt(den); //--- return(true); } //+------------------------------------------------------------------+ //| MathCorrelationKendall | //+------------------------------------------------------------------+ //| The function calculates Kendall's tau correlation coefficient. | //| | //| Arguments: | //| array1[] : First array with integer values | //| array2[] : Second array with integer values | //| tau : Variable for calculated value of the Kendall's tau | //| | //| Return value: true if successful, otherwise false | //+------------------------------------------------------------------+ bool MathCorrelationKendall(const int &array1[],const int &array2[],double &tau) { tau=QNaN; int size=ArraySize(array1); if(size==0 || ArraySize(array2)!=size) return(false); //--- int cnt1=0,cnt2=0,cnt=0; //--- for(int i=0; i0.0) cnt++; else cnt--; } } } //--- calculate Kendall tau double den=cnt1*cnt2; if(den==0) return(false); tau=cnt/MathSqrt(den); //--- return(true); } //+------------------------------------------------------------------+ //| MathQuantile | //+------------------------------------------------------------------+ //| The function produces sample quantiles corresponding to the | //| given probabilities. The smallest observation corresponds | //| to a probability of 0.0 and the largest to a probability of 1.0. | //+------------------------------------------------------------------+ //| The function calculates estimates of underlying distribution | //| quantiles based on one or two order statistics from the supplied | //| elements in x at probabilities in probs. | //| | //| Reference: | //| Hyndman, R. J. and Fan, Y. (1996) | //| "Sample quantiles in statistical packages", | //| American Statistician, vol.50, pp. 361–365. | //+------------------------------------------------------------------+ //| All sample quantiles are defined as weighted averages | //| of consecutive order statistics. | //| Sample quantiles of type i are defined by: | //| Q[i](p) = (1 - gamma)*x[j] + gamma*x[j+1] | //| | //| The function MathQuantile calculates quantiles according to | //| default type in R (type=7) | //| p(k) = (k - 1)/(n - 1). | //| In this case, p(k) = mode[F(x[k])]. This is used by S. | //+------------------------------------------------------------------+ bool MathQuantile(const double &array[],const double &probs[],double &quantile[]) { int size=ArraySize(array); int size_p=ArraySize(probs); if(size==0 || size_p==0) return(false); //--- check probability for(int i=0; i1.0) return(false); } //--- prepare array with values double x_values[]; if(ArrayCopy(x_values,array,0,0,WHOLE_ARRAY)!=size) return(false); //--- prepare output array if(ArraySize(quantile)lo[i]) { double gamma=(ind[i]-lo[i]); quantile[i]=(1.0-gamma)*quantile[i]+gamma*x_values[hi[i]]; } } //--- return(true); } //+------------------------------------------------------------------+ //| MathProbabilityDensityEmpirical | //+------------------------------------------------------------------+ //| The function calculates the empirical probability density | //| function (pdf) for random values from array[]. | //| | //| Arguments: | //| array[] : Array with random values | //| count : Data count, total count pairs (x,pdf(x)) | //| x[] : Output array for x values | //| pdf[] : Output array for empirical pdf(x) values | //| | //| Return value: true if successful, otherwise false | //+------------------------------------------------------------------+ bool MathProbabilityDensityEmpirical(const double &array[],const int count,double &x[],double &pdf[]) { if(count<=1) return(false); int size=ArraySize(array); if(size==0) return(false); //--- check NaN values for(int i=0; i