Discussion of article "Statistical Distributions in MQL5 - taking the best of R" - page 20

New comment
 

These are the results.



Not very good either... In general, working with non-central beta distribution rather belongs to the class of non-trivial tasks. We managed to see Boost from the open source code. And it will be more complicated than the one we have in SB now.

 
Thanks, we'll check it out.
 
Denis Kirichenko #:

The same Wikipedia has it:

I.e. two random variables are taken, where the first is from a non-central chi-square distribution and the second is from a central distribution. Probably the code can be corrected to this:

//--- generate random number using Noncentral ChiSquare
double chi1=MathRandomNoncentralChiSquare(a2,lambda2,error_code);
double chi2=MathRandomChiSquare(b2,error_code);
 result[i]=chi1/(chi1+chi2);

The modified graphs in the example will be below.

Thank you, you are right, we need to calculate through this formula, however there was another error, the lambda non-centrality parameter should be without multiplier 2.

//+------------------------------------------------------------------+
//| Random variate from the Noncentral Beta distribution |
//+------------------------------------------------------------------+
//| Compute the random variable from the Noncentral Beta |
//| distribution with parameters a,b and lambda. ||
//||
//| Arguments:|
//| a : First shape parameter |
//| b : Second shape parameter |
//| lambda : Noncentrality parameter |
//| error_code : Variable for error code |
//||
//| Return value:|
//| The random value with Noncentral Beta distribution. ||
//+------------------------------------------------------------------+
double MathRandomNoncentralBeta(const double a,const double b,const double lambda,int &error_code)
  {
//--- if lambda==0, return beta random variate
   if(lambda==0.0)
      return MathRandomBeta(a,b,error_code);
//--- check parameters
   if(!MathIsValidNumber(a) || !MathIsValidNumber(b) || !MathIsValidNumber(lambda))
     {
      error_code=ERR_ARGUMENTS_NAN;
      return QNaN;
     }
//--- a,b,lambda must be positive
   if(a<=0.0 || b<=0.0 || lambda<0)
     {
      error_code=ERR_ARGUMENTS_INVALID;
      return QNaN;
     }

   error_code=ERR_OK;
//--- generate random number using Noncentral ChiSquare
   double chi1=MathRandomNoncentralChiSquare(2*a,lambda,error_code);
   double chi2=MathRandomChiSquare(2*b,error_code);
   return chi1/(chi1+chi2);
  }
//+------------------------------------------------------------------+
//| Random variate from the Noncentral Beta distribution |
//+------------------------------------------------------------------+
//|| Generates random variables from the Noncentral Beta distribution |
//| with parameters a,b, lambda.|
//||
//| Arguments:|
//| a : First shape parameter |
//| b : Second shape parameter |
//| lambda : Noncentrality parameter |
//| data_count : Number of values needed |
//| result : Output array with random values |
//||
//| Return value:|
//| true if successful, otherwise false. ||
//+------------------------------------------------------------------+
bool MathRandomNoncentralBeta(const double a,const double b,const double lambda,const int data_count,double &result[])
  {
//--- if lambda==0, return beta random variate
   if(lambda==0.0)
      return MathRandomBeta(a,b,data_count,result);
//--- check parameters
   if(!MathIsValidNumber(a) || !MathIsValidNumber(b) || !MathIsValidNumber(lambda))
      return false;
//--- a,b,lambda must be positive
   if(a<=0.0 || b<=0.0 || lambda<0)
      return false;

   int error_code=0;
//--- prepare output array and calculate random values
   ArrayResize(result,data_count);
   double a2=a*2;
   double b2=b*2;
   for(int i=0; i<data_count; i++)
     {
      //--- generate random number using Noncentral ChiSquare
      double chi1=MathRandomNoncentralChiSquare(a2,lambda,error_code);
      double chi2=MathRandomChiSquare(b2,error_code);
      result[i]=chi1/(chi1+chi2);
     }
   return true;
  }

Test script:

//+------------------------------------------------------------------+
//|TestNoncentralBeta.mq5 |
//|Copyright 2024, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <Graphics\Graphic.mqh>
#include <Math\Stat\NoncentralBeta.mqh>
#include <Math\Stat\Math.mqh>
#property script_show_inputs
//--- input parameters
input double a_par=2;    // first parameter of beta distribution (shape1)
input double b_par=5;    // second beta distribution parameter (shape2)
input double l_par=1;    // non-centrality parameter (lambda)
//+------------------------------------------------------------------+
//| Script programme start function|
//+------------------------------------------------------------------+
void OnStart()
  {
//--- disable price chart display
   ChartSetInteger(0,CHART_SHOW,false);
//--- initialise the random number generator 
   MathSrand(GetTickCount());
//--- generate a sample of the random variable
   long chart=0;
   string name="GraphicNormal";
   int n=1000000;       // number of values in the sample
   int ncells=52;       // number of intervals in the histogram
   double x[];          // histogram interval centres
   double y[];          // number of values from the sample that fall into the interval
   double data[];       // sampling of random values 
   double max,min;      // maximum and minimum values in the sample
//--- get a sample from the non-central beta distribution 
   MathRandomNoncentralBeta(a_par,b_par,l_par,n,data);
//--- calculate data for histogram construction
   CalculateHistogramArray(data,x,y,max,min,ncells);
//--- obtain the sequence boundaries and the step for the theoretical curve construction
   double step;
   GetMaxMinStepValues(max,min,step);
   step=MathMin(step,(max-min)/ncells);
//--- we obtain theoretically calculated data on the interval [min,max]
   double x2[];
   double y2[];
   MathSequence(min,max,step,x2);
   MathProbabilityDensityNoncentralBeta(x2,a_par,b_par,l_par,false,y2);
//--- scaling
   double theor_max=y2[ArrayMaximum(y2)];
   double sample_max=y[ArrayMaximum(y)];
   double k=sample_max/theor_max;
   for(int i=0; i<ncells; i++)
      y[i]/=k;
//--- display graphs
   CGraphic graphic;
   if(ObjectFind(chart,name)<0)
      graphic.Create(chart,name,0,0,0,780,380);
   else
      graphic.Attach(chart,name);
   graphic.BackgroundMain(StringFormat("Noncentral Beta distribution alpha=%G beta=%G lambda=%G",
                          a_par,b_par,l_par));
   graphic.BackgroundMainSize(16);
//--- plot all curves
   graphic.CurveAdd(x,y,CURVE_HISTOGRAM,"Sample").HistogramWidth(6);
//--- and now let's plot the theoretical distribution density curve 
   graphic.CurveAdd(x2,y2,CURVE_LINES,"Theory");
   graphic.CurvePlotAll();
//--- plot all curves
   graphic.Update();
   
   DebugBreak();
  }
//+------------------------------------------------------------------+
//| Calculate frequencies for data set|
//+------------------------------------------------------------------+
bool CalculateHistogramArray(const double &data[],double &intervals[],double &frequency[],
                             double &maxv,double &minv,const int cells=10)
  {
   if(cells<=1) return (false);
   int size=ArraySize(data);
   if(size<cells*10) return (false);
   minv=data[ArrayMinimum(data)];
   maxv=data[ArrayMaximum(data)];
   double range=maxv-minv;
   double width=range/cells;
   if(width==0) return false;
   ArrayResize(intervals,cells);
   ArrayResize(frequency,cells);
//--- set the centres of the intervals
   for(int i=0; i<cells; i++)
     {
      intervals[i]=minv+(i+0.5)*width;
      frequency[i]=0;
     }
//--- fill in the interval frequencies
   for(int i=0; i<size; i++)
     {
      int ind=int((data[i]-minv)/width);
      if(ind>=cells) ind=cells-1;
      frequency[ind]++;
     }
   return (true);
  }
//+------------------------------------------------------------------+
//| Calculates values for sequence generation |
//+------------------------------------------------------------------+
void GetMaxMinStepValues(double &maxv,double &minv,double &stepv)
  {
//--- calculate the absolute magnitude of the sequence to get the normalisation accuracy
   double range=MathAbs(maxv-minv);
   int degree=(int)MathRound(MathLog10(range));
//--- normalise max. and min. values with specified accuracy
   maxv=NormalizeDouble(maxv,degree);
   minv=NormalizeDouble(minv,degree);
//--- step of sequence generation also set from the given accuracy
   stepv=NormalizeDouble(MathPow(10,-degree),degree);
   if((maxv-minv)/stepv<10)
      stepv/=10.;
  }

Result:

NoncentralBeta(2,5,1)

NoncentralBeta(2,5,10)

NoncentralBeta(2,15,100)

NoncentralBeta(2,15,20)

Files:
NoncentralBeta.mqh  41 kb
 

I also checked the CFD function for Noncentral Beta Distribution. 

#include <Math\Stat\NoncentralBeta.mqh>
//+------------------------------------------------------------------+
//| Script programme start function|
//+------------------------------------------------------------------+
void OnStart()
  {
   matrix nc_beta_params_mx =
     {
      // A B LAMBDA X
        {5, 5, 54, 0.864},
        {5, 5, 140, 0.9},
        {5, 5, 170, 0.956},
        {10, 10, 54, 0.8686},
        {10, 10, 140, 0.9},
        {10, 10, 250, 0.9},
        {20, 20, 54, 0.8787},
        {20, 20, 140, 0.9},
        {20, 20, 250, 0.922},
        {10, 20, 150, 0.868},
        {10, 10, 120, 0.9},
        {15, 5, 80, 0.88},
        {20, 10, 110, 0.85},
        {20, 30, 65, 0.66},
        {20, 50, 130, 0.72},
        {30, 20, 80, 0.72},
        {30, 40, 130, 0.8},
        {10, 5, 20, 0.644},
        {10, 10, 54, 0.7},
        {10, 30, 80, 0.78},
        {15, 20, 120, 0.76},
        {10, 5, 55, 0.795},
        {12, 17, 64, 0.56},
        {30, 30, 140, 0.8},
        {35, 30, 20, 0.67}
     };
// print
   PrintFormat("A - B - LAMBDA - X - CDF");
   for(ulong row = 0; row < nc_beta_params_mx.Rows(); row++)
     {
      double a_par, b_par, l_par, x_val;
      a_par = nc_beta_params_mx[row][0];
      b_par = nc_beta_params_mx[row][1];
      l_par = nc_beta_params_mx[row][2];
      x_val = nc_beta_params_mx[row][3];
      int error_code;
      double res_val = MathCumulativeDistributionNoncentralBeta(x_val, a_par,
                       b_par, l_par, error_code);
      ::PrintFormat("%0.0f, %0.0f, %0.0f, %0.2f -->  %0.5f",
                    a_par, b_par, l_par, x_val, res_val);
     }
  }


Some questionable results:

A - B - LAMBDA - X - CDF
5, 5, 54, 0.86 -->  0.74058
5, 5, 140, 0.90 -->  0.40096
5, 5, 170, 0.96 -->  1.00000
10, 10, 54, 0.87 -->  1.00000
10, 10, 140, 0.90 -->  0.99296
10, 10, 250, 0.90 -->  1.00000
20, 20, 54, 0.88 -->  1.00000
20, 20, 140, 0.90 -->  1.00000
20, 20, 250, 0.92 -->  1.00000
10, 20, 150, 0.87 -->  1.00000
10, 10, 120, 0.90 -->  1.00000
15, 5, 80, 0.88 -->  1.00000
20, 10, 110, 0.85 -->  1.00000
20, 30, 65, 0.66 -->  1.00000
20, 50, 130, 0.72 -->  1.00000
30, 20, 80, 0.72 -->  1.00000
30, 40, 130, 0.80 -->  1.00000
10, 5, 20, 0.64 -->  0.05069
10, 10, 54, 0.70 -->  0.18781
10, 30, 80, 0.78 -->  1.00000
15, 20, 120, 0.76 -->  1.00000
10, 5, 55, 0.80 -->  0.13001
12, 17, 64, 0.56 -->  0.00231
30, 30, 140, 0.80 -->  1.00000
35, 30, 20, 0.67 -->  0.88671


Alternative sources use the ASA310 algorithm:

23 July 2024 05:23:28 PM

ASA310_PRB:
  C++ version
  Test the ASA310 library.

TEST01:
  NCBETA computes the noncentral incomplete Beta function.
  Compare to tabulated values.

      A        B     LAMBDA        X          FX                        FX2
                                              (Tabulated)               (NCBETA)            DIFF

        5        5       54       0.864                 0.4563021        0.4563022483390367   1.483 e-07
        5        5      140         0.9                 0.1041337        0.1041334790875982   2.209 e-07
        5        5      170       0.956                 0.6022353        0.6022418158548807   6.516 e-06
       10       10       54      0.8686                  0.918777        0.9187759254330974   1.075 e-06
       10       10      140         0.9                 0.6008106        0.6008067894096832   3.811 e-06
       10       10      250         0.9                  0.090285        0.0902898993260559   4.899 e-06
       20       20       54      0.8787                 0.9998655        0.9998614048611703   4.095 e-06
       20       20      140         0.9                 0.9925997        0.9925954777347362   4.222 e-06
       20       20      250       0.922        0.9641111999999999        0.9641179545701509   6.755 e-06
       10       20      150       0.868              0.9376626573        0.9376626555219731   1.778 e-09
       10       10      120         0.9              0.7306817858         0.730681784479193   1.321 e-09
       15        5       80        0.88              0.1604256918        0.1604256916919781    1.08 e-10
       20       10      110        0.85              0.1867485313         0.186748530948068   3.519 e-10
       20       30       65        0.66        0.6559386874000001        0.6559386873992288   7.713 e-13
       20       50      130        0.72              0.9796881486        0.9796881474202076    1.18 e-09
       30       20       80        0.72              0.1162386423        0.1162386421523797   1.476 e-10
       30       40      130         0.8              0.9930430054        0.9930430042307262   1.169 e-09
       10        5       20       0.644              0.0506899273       0.05068987981355273   4.749 e-08
       10       10       54         0.7              0.1030959706        0.1030959706112684   1.127 e-11
       10       30       80        0.78        0.9978417832000001        0.9978417830424838   1.575 e-10
       15       20      120        0.76              0.2555552369        0.2555552355854933   1.315 e-09
       10        5       55       0.795              0.0668307064        0.0668307064085136   8.514 e-12
       12       17       64        0.56              0.0113601067       0.01136010666591296   3.409 e-11
       30       30      140         0.8              0.7813366615        0.7813366591185901   2.381 e-09
       35       30       20        0.67              0.8867126477        0.8867125569354448   9.076 e-08

ASA310_PRB:
  Normal end of execution.


For example some of the results in Wolfram Mathematica:

Table[CDF[NoncentralBetaDistribution[5,5,140],0.9]]
= 0.1041335
Table[CDF[NoncentralBetaDistribution[15,5,80],0.88]]
= 0.1604257
Table[CDF[NoncentralBetaDistribution[35, 30, 20], 0.67]]
= 0.8867126

Which is more or less similar to the truth...

1...1314151617181920
New comment