記事"R で 統計分布を MQL5 に -"についてのディスカッション - ページ 20

新しいコメント
 

これがその結果だ。



どちらもあまりよくないが...。一般的に、非集中ベータ分布での作業は、むしろ非自明なタスクに属する。我々はオープンソースのコードからブーストを 見ることができた。そしてそれは、今SBにあるものよりも複雑になるだろう。

 
ありがとう。
 
Denis Kirichenko #:

同じウィキペディアにこうある：

つまり、2つの確率変数が取られ、1つ目は非心カイ2乗分布から、2つ目は中心分布から 取られる。おそらくコードはこのように修正できるだろう：

//--- 非心カイ二乗を用いて乱数を生成する。
double chi1=MathRandomNoncentralChiSquare(a2,lambda2,error_code);
double chi2=MathRandomChiSquare(b2,error_code);
 result[i]=chi1/(chi1+chi2);

例のグラフを修正したものが以下になります。

ありがとうございます、その通りです、この式で計算する必要があります、しかしもう一つ間違いがありました、ラムダ非中心性パラメータは乗数2なしにすべきです。

//+------------------------------------------------------------------+
//| 非心ベータ分布の確率変分
//+------------------------------------------------------------------+
//| 非心ベータから確率変数を計算する。
//| パラメータa,bとλを持つ分布。
//||
//| 引数：|
//| a :最初の形状パラメータ
//| b :第2形状パラメータ
//| ラムダ ：非中心性パラメータ
//| error_code ：エラーコード用変数
//||
//| 返り値：|
//| 非心ベータ分布による乱数値。
//+------------------------------------------------------------------+
double MathRandomNoncentralBeta(const double a,const double b,const double lambda,int &error_code)
  {
//--- lambda==0の場合、ベータ値を返す
   if(lambda==0.0)
      return MathRandomBeta(a,b,error_code);
//--- パラメータをチェックする
   if(!MathIsValidNumber(a) || !MathIsValidNumber(b) || !MathIsValidNumber(lambda))
     {
      error_code=ERR_ARGUMENTS_NAN;
      return QNaN;
     }
//--- a,b,ラムダは正でなければならない
   if(a<=0.0 || b<=0.0 || lambda<0)
     {
      error_code=ERR_ARGUMENTS_INVALID;
      return QNaN;
     }

   error_code=ERR_OK;
//--- 非心カイ二乗を用いて乱数を生成する。
   double chi1=MathRandomNoncentralChiSquare(2*a,lambda,error_code);
   double chi2=MathRandomChiSquare(2*b,error_code);
   return chi1/(chi1+chi2);
  }
//+------------------------------------------------------------------+
//| 非心ベータ分布の確率変分
//+------------------------------------------------------------------+
//|| 非心ベータ分布から確率変数を生成する。
//| パラメータa,b,ラムダを持つ。|
//||
//| 引数：|
//| a :最初の形状パラメータ
//| b :第2形状パラメータ
//| ラムダ ：非中心性パラメータ
//| data_count : 必要な値の数
//| result : 乱数値を含む出力配列 |.
//||
//| 返り値：|
//| 成功すればtrue、そうでなければfalse。
//+------------------------------------------------------------------+
bool MathRandomNoncentralBeta(const double a,const double b,const double lambda,const int data_count,double &result[])
  {
//--- lambda==0の場合、ベータ値を返す
   if(lambda==0.0)
      return MathRandomBeta(a,b,data_count,result);
//--- パラメータをチェックする
   if(!MathIsValidNumber(a) || !MathIsValidNumber(b) || !MathIsValidNumber(lambda))
      return false;
//--- a,b,ラムダは正でなければならない
   if(a<=0.0 || b<=0.0 || lambda<0)
      return false;

   int error_code=0;
//--- 出力配列を準備し、乱数値を計算する
   ArrayResize(result,data_count);
   double a2=a*2;
   double b2=b*2;
   for(int i=0; i<data_count; i++)
     {
      //--- 非心カイ二乗を用いて乱数を生成する。
      double chi1=MathRandomNoncentralChiSquare(a2,lambda,error_code);
      double chi2=MathRandomChiSquare(b2,error_code);
      result[i]=chi1/(chi1+chi2);
     }
   return true;
  }

テストスクリプト

//+------------------------------------------------------------------+
//|TestNoncentralBeta.mq5
//|著作権 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 double a_par=2;    // ベータ分布の最初のパラメータ (shape1)
input double b_par=5;    // 2番目のベータ分布パラメータ (shape2)
input double l_par=1;    // 非中心性パラメータ (ラムダ)
//+------------------------------------------------------------------+
//| スクリプト番組開始機能|
//+------------------------------------------------------------------+
void OnStart()
  {
//--- 価格チャートを表示しない
   ChartSetInteger(0,CHART_SHOW,false);
//--- 乱数発生器を初期化する 
   MathSrand(GetTickCount());
//--- 確率変数のサンプルを生成する
   long chart=0;
   string name="GraphicNormal";
   int n=1000000;       // サンプル内の値の数
   int ncells=52;       // ヒストグラムの区間数
   double x[];          // ヒストグラム・インターバル・センター
   double y[];          // 区間に入るサンプルの値の数
   double data[];       // ランダムな値のサンプリング 
   double max,min;      // サンプルの最大値と最小値
//--- 非心ベータ分布からサンプルを得る 
   MathRandomNoncentralBeta(a_par,b_par,l_par,n,data);
//--- ヒストグラムを作成するためのデータを計算する
   CalculateHistogramArray(data,x,y,max,min,ncells);
//--- シーケンス境界と理論曲線構築のためのステップを得る
   double step;
   GetMaxMinStepValues(max,min,step);
   step=MathMin(step,(max-min)/ncells);
//--- 区間[min,max]上の理論的に計算されたデータが得られる。
   double x2[];
   double y2[];
   MathSequence(min,max,step,x2);
   MathProbabilityDensityNoncentralBeta(x2,a_par,b_par,l_par,false,y2);
//--- スケーリング
   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;
//--- グラフを表示する
   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);
//--- すべての曲線をプロットする
   graphic.CurveAdd(x,y,CURVE_HISTOGRAM,"Sample").HistogramWidth(6);
//--- そして次に、理論的な分布密度曲線をプロットしてみよう。 
   graphic.CurveAdd(x2,y2,CURVE_LINES,"Theory");
   graphic.CurvePlotAll();
//--- すべての曲線をプロットする
   graphic.Update();
   
   DebugBreak();
  }
//+------------------------------------------------------------------+
//| データセットの度数を計算する|
//+------------------------------------------------------------------+
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);
//--- 区間の中心を設定する
   for(int i=0; i<cells; i++)
     {
      intervals[i]=minv+(i+0.5)*width;
      frequency[i]=0;
     }
//--- 間隔周波数を埋める
   for(int i=0; i<size; i++)
     {
      int ind=int((data[i]-minv)/width);
      if(ind>=cells) ind=cells-1;
      frequency[ind]++;
     }
   return (true);
  }
//+------------------------------------------------------------------+
//| シーケンス生成のための値を計算する
//+------------------------------------------------------------------+
void GetMaxMinStepValues(double &maxv,double &minv,double &stepv)
  {
//--- 正規化精度を得るために、シーケンスの絶対的な大きさを計算する。
   double range=MathAbs(maxv-minv);
   int degree=(int)MathRound(MathLog10(range));
//--- 指定された精度で最大値と最小値を正規化する。
   maxv=NormalizeDouble(maxv,degree);
   minv=NormalizeDouble(minv,degree);
//--- シーケンス生成のステップも、与えられた精度から設定される
   stepv=NormalizeDouble(MathPow(10,-degree),degree);
   if((maxv-minv)/stepv<10)
      stepv/=10.;
  }

結果

非心ベータ(2,5,1)

非心ベータ(2,5,10)

非心ベータ(2,15,100)

非心ベータ(2,15,20)

ファイル:
NoncentralBeta.mqh  41 kb
 

非心ベータ分布の CFD関数もチェックした。 

#include <Math\Stat\NoncentralBeta.mqh>
//+------------------------------------------------------------------+
//| スクリプト番組開始機能|
//+------------------------------------------------------------------+
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}
     };
// プリント
   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);
     }
  }


疑わしい結果もある：

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


別の情報源では，ASA310アルゴリズムが使われています．

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.


例えば，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

これは多かれ少なかれ真実に近い...

1...1314151617181920
新しいコメント