//+------------------------------------------------------------------+
//|                                                    Inference.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <CategoricalHMM\CategoricalHMM.mqh>
#include <Graphics\Graphic.mqh>
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   int n_states = 2;    // Количество скрытых состояний
   int n_emissions = 6; // Количество возможных исходов (количество классов)

// --- 1. Создаем модель (seed=99 для повторяемости) ---
   CCategoricalHMM model(n_states, n_emissions, 99);

// --- Задаем параметры модели ---
// --- 2. Матрица переходов TR ---
   matrix tr = {{0.95, 0.05},
        {0.1,  0.9}
     };
// --- 3.Матрицу эмиссий (E) ---
   matrix e = {{1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6},
        {0.1,   0.1,   0.1,   0.1,   0.1,   0.5}
     };

// --- 4. Генерируем последовательность наблюдений и скрытых состояний ---
   vector y;
   vector z;
   int n_samples = 300; // Количество генерируемых выборок
   model.Sample(n_samples, tr, e, y, z);

   PrintFirst(z,"First 5 states",5);       // First 5 states:     [1 1 1 1 1]
   PrintFirst(y,"First 5 emissions",5);    // First 5 emissions:  [6 1 5 2 3]

// --- БЛОК ВИЗУАЛИЗАЦИИ ---
   int plotLimit = 300; // Рисуем первые 300 точек
   double x[], y_[];
   ArrayResize(x, plotLimit);
   ArrayResize(y_, plotLimit);
   for(int i = 0; i < plotLimit; i++)
     {
      x[i] = i;
      y_[i] = z[i];
     }
   ChartSetInteger(0,CHART_SHOW,false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0,CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0,CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0,"HMM_States_Plot",0,0,0,int(width),int(height));
   CCurve *curve = graphic.CurveAdd(x, y_, ColorToARGB(clrDodgerBlue, 255), CURVE_LINES,"States");
   curve.LinesStyle(STYLE_SOLID);
   curve.LinesWidth(2);
   graphic.BackgroundMainColor(ColorToARGB(clrBlack,255));
   graphic.BackgroundMainSize(24);
   graphic.BackgroundMain("States over time");
   graphic.XAxis().Name("Time");
   graphic.XAxis().NameSize(18);
   graphic.YAxis().Name("State (1=Fair, 2=Loaded)");
   graphic.YAxis().NameSize(18);
   graphic.YAxis().Min(1);
   graphic.YAxis().Max(2);
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(5*1000);
   ChartSetInteger(0,CHART_SHOW,true);
   graphic.Destroy();
   ChartRedraw(0);

// --- Эмпирическая вероятность наблюдения=6 для симметричной(Fair) и смещенной (Loaded) кости ---
   PrintFormat("Empirical probability '6' (Fair): %.3f", CalculateSixFreq(y, z, 1));
   PrintFormat("Empirical probability '6' (Loaded): %.3f", CalculateSixFreq(y, z, 2));

// --- Inference ---
   matrix Filtering, Smoothing, bs;
   double marginal_loglik;
   vector s;

   if(model.Inference(y,tr,e,Smoothing,marginal_loglik,Filtering,bs,s))
     {
      Print("---Inference---");
      PrintFormat("Логарифм правдоподобия (marginal_loglik): %.10f", marginal_loglik);
      Print("\nМатрица Smoothing  (апостериорные вероятности):");
      PrintMatrixFormatted(Smoothing);
      Print("\nМатрица Filtering :");
      PrintMatrixFormatted(Filtering);
     }
//MatrixToCSV("FilteringDist.csv",Filtering);
//MatrixToCSV("SmoothingDist.csv",Smoothing);
//VectorToCSV("TrueStates.csv", z);

  }
//+------------------------------------------------------------------+

//+-------------------------------------------------------------------+
//|Возвращает частоту выпадения шестёрки                              |                                                                |
//+-------------------------------------------------------------------+
double CalculateSixFreq(const vector &r, const vector &s, int stateTarget)
  {
   int count = 0, total = 0;
   for(int i=0; i<(int)r.Size(); i++)
     {
      if(s[i] == stateTarget)
        {
         total++;
         if(r[i] == 6)
            count++;
        }
     }
   return (total > 0) ? (double)count/total : 0;
  }
//+------------------------------------------------------------------+