Ошибки в библиотеке статистики

Новый комментарий
 
Полагаю нелишним завести отдельную ветку для выявления, локализации и предложений по исправлению ошибок в MQL5-библиотеке статистики.
 
Предлагаю начать с ошибки вычисления квантилей Т-распределения Стьюдента. Уже было об этом на форуме, но повторю ещё раз.
#include <Math\Stat\T.mqh>
//+------------------------------------------------------------------+
void OnStart()
{ int err = 0;
  double qT = MathQuantileT(0.05, 100, err);
  if (err != 0)
    Print("MathQuantileT() error ", err);
  Print(qT);
}

MathQuantileT() error 4
nan


должно быть около -1.660234

 
Aleksey Nikolayev #:
Предлагаю начать с ошибки вычисления квантилей Т-распределения Стьюдента

Собственно вычислениями занимается эта функция из T.mqh

double MathQuantileT(const double probability,const double nu,const bool tail,const bool log_mode,int &error_code)
  {
//--- check NaN
   if(!MathIsValidNumber(probability) || !MathIsValidNumber(nu))
     {
      error_code=ERR_ARGUMENTS_NAN;
      return QNaN;
     }
//--- check nu
   if(nu!=MathRound(nu) || nu<0.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;
     }

//--- check cases when probability==0 or 1
   if(prob==0.0 || prob==1.0)
     {
      error_code=ERR_RESULT_INFINITE;
      //---
      if(prob==0.0)
         return(QNEGINF);
      else
         return(QPOSINF);
     }

   error_code=ERR_OK;
//--- special case nu=1
   if(nu==1.0)
      return MathTan(M_PI*(prob-0.5));
//--- special case
   if(prob==0.5)
      return 0.0;
//---
   int    max_iterations=50;
   int    iterations=0;
//--- initial values
   double h=1.0;
   double h_min=10E-20;
   double x=0.5;
   int    err_code=0;
//--- Newton iterations
   while(iterations<max_iterations)
     {
      //--- check convegence
      if((MathAbs(h)>h_min && MathAbs(h)>MathAbs(h_min*x))==false)
         break;
      //--- calculate pdf and cdf
      double pdf=MathProbabilityDensityT(x,nu,err_code);
      double cdf=MathCumulativeDistributionT(x,nu,err_code);
      //--- calculate ratio
      h=(cdf-prob)/pdf;
      //---
      double x_new=x-h;
      //--- check x
      if(x_new<0.0)
         x_new=x*0.1;
      else
      if(x_new>1.0)
         x_new=1.0-(1.0-x)*0.1;

      x=x_new;

      iterations++;
     }
//--- check convergence
   if(iterations<max_iterations)
      return x;
   else
     {
      error_code=ERR_NON_CONVERGENCE;
      return QNaN;
     }
  }

Ошибка с кодом 4 - это ERR_NON_CONVERGENCE - не сходится метод Ньютона.

 
Aleksey Nikolayev #:
Ошибка с кодом 4 - это ERR_NON_CONVERGENCE - не сходится метод Ньютона.

Имхо, если целью была простота кода, то лучше было сделать поиск корня делением пополам.
Для Ньютона нужно аккуратно выбирать начальную точку, аккуратнее относиться к близости плотности к нулю в хвостах и тд.

С помощью ИИ набросал на скорую руку. На первый взгляд вроде работает - особо тщательно не проверял, поскольку код просто для примера

double MathQuantileT(const double probability, const double nu, const bool tail, const bool log_mode, int &error_code)
{
   // 1. Проверка на NaN
   if(!MathIsValidNumber(probability) || !MathIsValidNumber(nu))
   {
      error_code = ERR_ARGUMENTS_NAN;
      return QNaN;
   }

   // 2. Проверка nu (убрали nu != MathRound(nu), так как df могут быть дробными)
   if(nu <= 0.0)
   {
      error_code = ERR_ARGUMENTS_INVALID;
      return QNaN;
   }

   // 3. Расчет реальной вероятности (учитываем хвосты и логарифмы)
   double prob = TailLogProbability(probability, tail, log_mode);

   // 4. Проверка диапазона вероятности
   if(prob < 0.0 || prob > 1.0)
   {
      error_code = ERR_ARGUMENTS_INVALID;
      return QNaN;
   }

   // 5. Обработка предельных случаев
   if(prob == 0.0) return QNEGINF;
   if(prob == 1.0) return QPOSINF;
   if(prob == 0.5) return 0.0;

   // 6. Спецкейс для nu=1 (распределение Коши)
   if(nu == 1.0) return MathTan(M_PI * (prob - 0.5));

   error_code = ERR_OK;
   int err_code = 0;

   // 7. Метод бисекции (надежный поиск корня)
   double x_low = -1.0;
   double x_high = 1.0;
   
   // Динамически расширяем границы поиска, пока корень не окажется внутри [x_low, x_high]
   while(MathCumulativeDistributionT(x_low, nu, err_code) > prob && x_low > -1e6) 
      x_low *= 2.0;
   while(MathCumulativeDistributionT(x_high, nu, err_code) < prob && x_high < 1e6) 
      x_high *= 2.0;

   double x = 0.0;
   int max_iterations = 100;
   double precision = 1e-12; // Высокая точность

   for(int i = 0; i < max_iterations; i++)
   {
      x = (x_low + x_high) / 2.0;
      double cdf = MathCumulativeDistributionT(x, nu, err_code);

      if(MathAbs(x_high - x_low) < precision) 
         break;

      if(cdf < prob)
         x_low = x;
      else
         x_high = x;
         
      if(i == max_iterations - 1) error_code = ERR_NON_CONVERGENCE;
   }

   return x;
}
 
Aleksey Nikolayev #:
Предлагаю начать с ошибки вычисления квантилей Т-распределения Стьюдента. Уже было об этом на форуме, но повторю ещё раз.

MathQuantileT() error 4
nan


должно быть около -1.660234

Странно, что функция ожидает пять обязательных параметров

double MathQuantileT(const double probability, const double nu, const bool tail, const bool log_mode, int &error_code)

вы передаёте всего три

double qT = MathQuantileT(0.05, 100, err);

И при этом не возникает ошибки компиляции?

 
Roman #:
вы передаёте всего три

Их два варианта. Та что с тремя параметрами вызывает ту что с пятью (в Math\Stat\T.mqh)

double MathQuantileT(const double probability,const double nu,int &error_code)
  {
   return MathQuantileT(probability,nu,true,false,error_code);
  }
 
С помощью встроенной функции для Бета распределения чуть проще

#include <Math\Stat\Beta.mqh> 
#include <Math\Stat\Math.mqh> 

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {

 int err = 0;
  double qT = MathQuantileT(0.05, 100,true,false, err);
  if (err != 0)
    Print("MathQuantileT() error ", err);
  Print(qT);
   
  }
//+------------------------------------------------------------------+

double MathQuantileT(const double probability, const double nu, const bool tail, const bool log_mode, int &error_code)
{
   //--- 1. Проверка на NaN
   if(!MathIsValidNumber(probability) || !MathIsValidNumber(nu))
   {
      error_code = ERR_ARGUMENTS_NAN;
      return QNaN;
   }

   //--- 2. Проверка степеней свободы (nu > 0)
   // Примечание: Стьюдент определен и для нецелых nu (например, 2.5)
   if(nu <= 0.0)
   {
      error_code = ERR_ARGUMENTS_INVALID;
      return QNaN;
   }

   //--- 3. Вычисление реальной вероятности если передали probability как логарифм 
   double prob = TailLogProbability(probability, tail, log_mode);

   //--- 4. Проверка диапазона вероятности
   if(prob < 0.0 || prob > 1.0)
   {
      error_code = ERR_ARGUMENTS_INVALID;
      return QNaN;
   }

   //--- 5. Обработка предельных случаев
   if(prob == 0.0) return QNEGINF;
   if(prob == 1.0) return QPOSINF;
   if(prob == 0.5) return 0.0;

   //--- 6. Спецслучай nu=1 (Распределение Коши)
   if(nu == 1.0)
      return MathTan(M_PI * (prob - 0.5));

   error_code = ERR_OK;

   //--- 7. ОСНОВНОЙ РАСЧЕТ ЧЕРЕЗ BETA
   // Используем симметрию t-распределения
   bool is_lower = (prob < 0.5);
   double p_beta = is_lower ? 2.0 * prob : 2.0 * (1.0 - prob);

   // Вызываем встроенную функцию MQL5 для квантиля Бета-распределения
   int beta_err = 0;
   double x = MathQuantileBeta(p_beta, nu / 2.0, 0.5, beta_err);

   if(beta_err != 0 || !MathIsValidNumber(x))
   {
      error_code = ERR_NON_CONVERGENCE;
      return QNaN;
   }

   // Финальная трансформация: x из Beta в t-значение
   // Формула: t = sqrt( nu * (1-x) / x )
   double t = MathSqrt(nu * (1.0 - x) / x);

   return is_lower ? -t : t;
}
Проверил разные варианты, совпадает со значениями в матлаб.
 
Evgeniy Chernish #:
С помощью встроенной функции для Бета распределения чуть проще

Спасибо!

Пишите если найдёте ещё ошибки в библиотеке.

Вроде какие-то проблемы были с гипергеометрическим и биномиальным, но точно уже не помню. Наверно надо просто все функции протестить на рандомных данных и сравнить с результатом в R.

Библиотека вроде бы есть, а пользоваться страшновато)

 

Спасибо за сообщения, ошибки расчета квантиля исправлены, обновленная версия библиотеки в приложении.

Evgeniy Chernish #:
С помощью встроенной функции для Бета распределения чуть проще
Расчет квантиля T исправлен на этот вариант.
Aleksey Nikolayev #:

Имхо, если целью была простота кода, то лучше было сделать поиск корня делением пополам.
Для Ньютона нужно аккуратно выбирать начальную точку, аккуратнее относиться к близости плотности к нулю в хвостах и тд.

Для других методов при расчете квантилей вместо метода Ньютона теперь используется метод половинного деления.
Aleksey Nikolayev #:

должно быть около -1.660234

Текущая версия дает   -1.6602343260853307


Также исправлены граничные значения, тесты из MQL5\Scripts\UnitTests\Stat\ работают нормально.

В приложении скрипт теста TestStatConsistency, для всех распределений он проверяет внутреннюю согласованность всей библиотеки (scalar/array/log/tail/quantile должны быть взаимно совместимы)

  1. scalar vs array (функции для массивов дают тот же результат, что и одиночные вызовы)
  2. pdf  (плотность/вероятность неотрицательная)
  3. cdf  (CDF лежит в диапазоне [0,1])
  4. lower/upper tail  (условие lower_tail + upper_tail = 1)
  5. log_mode (условие exp(log_pdf) == pdf, exp(log_cdf) == cdf)
  6. quantile inverse  (условие q(cdf(x)) ≈ x, а также q(logcdf(x)) ≈ x)
  7. invalid arguments (некорректные параметры возвращают ошибку и NaN)
Файлы:
Stat.zip  86 kb
TestStatConsistency.mq5  73 kb
 
Спасибо!
Новый комментарий