//+------------------------------------------------------------------+
//|                                                 distribution.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"base.mqh"
#include<Math/Stat/Gamma.mqh>
#include<Arch/Utility/igami.mqh>
#include<Arch/Utility/igam.mqh>
//+------------------------------------------------------------------+
//| constraints container                                            |
//+------------------------------------------------------------------+
struct Constraints
  {
   matrix            _one;
   vector            _two;

                     Constraints(void)
     {
      _one = matrix::Zeros(0,0);
      _two = vector::Zeros(0);
     }
                    ~Constraints(void)
     {
     }
                     Constraints(Constraints &other)
     {
      _one = other._one;
      _two = other._two;
     }
   void              operator=(Constraints &other)
     {
      _one = other._one;
      _two = other._two;
     }
  };

//+------------------------------------------------------------------+
//| base class for distributions                                     |
//+------------------------------------------------------------------+
class CDistribution
  {
protected:
   ulong             m_num_params;
   vector            m_parameters;
   vector            m_vector;
   string            m_name;
   uint              m_seed;
   CHighQualityRandStateShell m_generator;
   ENUM_DISTRIBUTION_MODEL m_dist_model;
   bool              m_initialized;
   vector            _check_constraints(vector &params_)
     {
      matrix bounds_ = bounds(vector::Zeros(0));
      ulong nparams = params_.Size();
      if(nparams!=bounds_.Rows())
        {
         Print(__FUNCTION__, " error ", "parameters must have ",bounds_.Rows()," elements");
         return vector::Zeros(0);
        }
      if(!bounds_.Rows())
         return vector::Zeros(0);
      for(ulong i = 0; i<bounds_.Rows(); ++i)
        {
         if(!(bounds_[i,0]<=params_[i]<=bounds_[i,1]))
           {
            Print(__FUNCTION__, " error ", ". Bounds violated ");
            return vector::Zeros(0);
           }
        }
      return params_;
     }

   bool              _initialize_generator(int rstate, bool use_alglib=true)
     {
      if(use_alglib)
        {
         m_seed = uint(MathAbs(rstate));

         if(rstate)
            CAlglib::HQRndSeed((int)m_seed,(int)m_seed+1,m_generator);
         else
            CHighQualityRand::HQRndRandomize(m_generator.GetInnerObj());
         return true;

        }
      else
        {
         m_seed = uint(MathAbs(rstate));
         MathSrand(m_seed);
        }
      return true;
     }

   virtual bool      _initialize(vector& distribution_params, int seed=0)
     {
      return false;
     }
public:
                     CDistribution(void):m_name(NULL),
                     m_dist_model(WRONG_VALUE),
                     m_num_params(0),
                     m_initialized(false),
                     m_parameters(vector::Zeros(0))
     {
     }

                     CDistribution(vector& params, int seed=0)
     {
      initialize(params,seed);
     }

                     CDistribution(CDistribution& other)
     {
      m_name = other.name();
      m_dist_model = other.distribution_type();
      m_num_params = other.numParams();
      m_parameters = other.get_params();
      m_seed = other.randomState();
      m_initialized = other.is_initialized();
     }

   void              operator=(CDistribution& other)

     {
      m_name = other.name();
      m_dist_model = other.distribution_type();
      m_num_params = other.numParams();
      m_parameters = other.get_params();
      m_seed = other.randomState();
      m_initialized = other.is_initialized();
     }

                    ~CDistribution(void)
     {

     }
   virtual bool      initialize(vector& distribution_params, int seed=0)
     {
      if(is_initialized())
        {
         m_initialized = false;
        }

      return _initialize(distribution_params,seed);
     }
   virtual ENUM_DISTRIBUTION_MODEL distribution_type(void)
     {
      return m_dist_model;
     }
   virtual bool      is_initialized(void) { return m_initialized; }
   virtual CHighQualityRandStateShell* generator(void) { return GetPointer(m_generator); }
   virtual uint      randomState(void) { return m_seed; }
   virtual ulong     numParams(void)   { return m_num_params; }
   virtual string    name(void)        { return m_name; }
   virtual vector    get_params(void)  { return m_parameters; }
   void              set_params(vector &parameters)
     {
      m_parameters = parameters;
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    rng(ulong nsize)
     {
      matrix out = simulator(nsize);
      return out.Col(0);
     }
   virtual matrix    rng(ulong nrows,ulong ncols)
     {
      return simulator(nrows,ncols);
     }
   virtual Constraints constraints(void)
     {
      return Constraints();
     }
   virtual matrix    bounds(vector& resids)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individial = false)
     {
      return vector::Zeros(0);
     }
   virtual vector    startingValues(vector &resids)
     {
      return vector::Zeros(0);
     }
   virtual vector    ppf(vector &pits, vector &parameters)
     {
      return vector::Zeros(0);
     }
   virtual vector    cdf(vector &resids, vector &parameters)
     {
      return vector::Zeros(0);
     }
   virtual double    moment(int n, vector& parameters)
     {
      return EMPTY_VALUE;
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      return EMPTY_VALUE;
     }
  };
//+------------------------------------------------------------------+
//|Standard normal distribution for use with ARCH models             |
//+------------------------------------------------------------------+
class CNormal: public CDistribution
  {
private:

   int               factorial2(int k)
     {
      int out = k;
      for(int i = k-2; i>1; out*=i, i-=2);
      return out;
     }

   double            munp(ulong n)
     {
      if(!n)
         return 1.;
      if(MathMod(double(n),2.)==0)
         return (double)factorial2(int(n)-1);
      else
         return 0.;
     }
   vector            moment_from_stats(ulong n, double mu,double mu2, double g1, double g2)
     {
      vector out(1);
      if(n==0)
        {
         out[0] = 1.0;
         return out;
        }
      else
         if(n==1)
           {
            if(!mu)
               out[0] = munp(1);
            else
               out[0] = mu;
           }
         else
            if(n==2)
              {
               if(!mu2 || !mu)
                  out[0] = munp(2);
               else
                  out[0] = mu2+mu*mu;
              }
            else
               if(n == 3)
                 {
                  if(!g1 || !mu2)
                     out[0] = munp(3);
                  else
                    {
                     double mu3 = g1*pow(mu2,1.5);
                     out[1] = mu3+3.0*mu*mu2+mu*mu*mu;
                    }
                 }
               else
                  if(n==4)
                    {
                     if(!g1 || !g2)
                        out[0] = munp(4);
                     else
                       {
                        double mu3, mu4;
                        mu4 = (g2+3.0)*pow(mu2,2.0);
                        mu3 = g1*pow(mu2, 1.5);
                        out[1] = mu4+4*mu*mu3+6*mu*mu*mu2+mu*mu*mu*mu;
                       }
                    }
                  else
                     out[1] = munp(n);

      return out;

     }
   vector            norm_moment(ulong order, double loc = 0.0, double scale = 1.)
     {
      ulong n = order;
      bool i0 = (scale>0.);
      bool i1 = bool(!loc);
      bool i2 = bool(loc);

      double mu,mu1,g1, g2;
      mu = mu1 = g1 = g2 = 0.0;
      if(n>0 && n<5)
        {
         mu = 0.0;
         mu1 = 1.0;
         g1 = 0.0;
         g2 = 0.0;
        }
      vector val(1);
      val = moment_from_stats(n,mu,mu1,g1,g2);
      val[0] = pow(scale, order)*val[0];
      return val;
     }
   virtual bool      _initialize(vector& distribution_params, int seed=0) override
     {
      m_name = "Normal";
      m_parameters = distribution_params;
      m_num_params = 0;
      m_initialized =  _initialize_generator(seed,true);
      return m_initialized;
     }
public:
                     CNormal(void)
     {
      m_dist_model = DIST_NORMAL;
     }
                     CNormal(vector& params, int seed=0)
     {
      initialize(params,seed);
     }

                    ~CNormal(void)
     {
     }
   virtual Constraints  constraints(void)  override
     {
      Constraints out;
      return out;
     }
   virtual matrix    bounds(vector &resids)
     {
      return matrix::Zeros(0,0);
     }
   virtual vector    loglikelihood(vector &parameters,vector& resids, vector& sigmas,bool individual = false) override
     {
      vector lls = -0.5 * (log(2 * M_PI) + log(sigmas) + pow(resids, 2.0) / sigmas) ;
      if(individual)
         return lls;
      else
        {
         vector out(1);
         out[0] = lls.Sum();
         return out;
        }
     }
   virtual vector    startingValues(vector& resids) override
     {
      return vector::Zeros(0);
     }
   virtual matrix    simulator(ulong nsize,ulong ncols = 1) override
     {
      matrix out(nsize,ncols);
      CAlglib::HQRndNormalM(m_generator,int(nsize),int(ncols),out);
      return out;
     }
   virtual vector    ppf(vector &pits, vector &parameters) override
     {
      _check_constraints(parameters);
      vector out(pits.Size());
      for(ulong i = 0; i<out.Size(); CAlglib::InvNormalCDF(pits[i]), ++i);
      return out;
     }
   virtual vector    cdf(vector &resids, vector &parameters) override
     {
      set_params(parameters);
      _check_constraints(parameters);
      vector out(resids.Size());
      for(ulong i = 0; i<out.Size(); out[i] = CAlglib::NormalCDF(resids[i]), ++i);
      return out;
     }
   virtual double    moment(int n, vector& parameters)
     {
      set_params(parameters);
      if(n<0)
         return double("nan");
      vector out =  norm_moment(ulong(n));
      return out[0];
     }
   virtual double    partialMoment(int n, vector &parameters, double z = 0.0)
     {
      set_params(parameters);
      if(n<0)
         return double("nan");
      else
         if(n == 0)
            return CAlglib::NormalCDF(z);
         else
            if(n == 1)
               return -1.*CAlglib::NormalPDF(z);
            else
               return -pow(z,double(n-1))*CAlglib::NormalPDF(z) + double(n-1)* partialMoment(n-2,parameters,z);
     }
  };