//+------------------------------------------------------------------+
//|                                              MannWhitneyTest.mqh |
//+------------------------------------------------------------------+
#ifndef MANNWHITNEYTEST_MQH
#define MANNWHITNEYTEST_MQH

//+------------------------------------------------------------------+
//| Helper structure used while ranking the pooled sample            |
//+------------------------------------------------------------------+
struct SRankItem
  {
   double            value;
   int               index;
   double            rank;
  };

//+------------------------------------------------------------------+
//| Class CMannWhitneyTest                                           |
//| Purpose: implements the Mann-Whitney U test for two independent  |
//| samples, using the normal approximation with tie correction.     |
//+------------------------------------------------------------------+
class CMannWhitneyTest
  {
private:
   double            m_group1[];
   double            m_group2[];
   int               m_n1;
   int               m_n2;
   double            m_U1;
   double            m_U2;
   double            m_U;
   double            m_z_score;
   double            m_p_value;
   double            m_tie_sum;

   void              RankPooled(double &ranks[]);
   double            NormalCDF(double z);
   void              SortRankItems(SRankItem &items[]);

public:
                     CMannWhitneyTest(void);
                    ~CMannWhitneyTest(void);

   void              SetGroup1(const double &data[]);
   void              SetGroup2(const double &data[]);
   bool              Run(void);
   double            GetU(void) const;
   double            GetZ(void) const;
   double            GetP(void) const;
   void              PrintResult(double alpha);
   void              ComputeQuartiles(const double &data[],double &q1,double &q2,double &q3);
  };

//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
CMannWhitneyTest::CMannWhitneyTest(void)
  {
   m_n1=0;
   m_n2=0;
   m_U1=0.0;
   m_U2=0.0;
   m_U=0.0;
   m_z_score=0.0;
   m_p_value=1.0;
   m_tie_sum=0.0;
  }

//+------------------------------------------------------------------+
//| Destructor                                                       |
//+------------------------------------------------------------------+
CMannWhitneyTest::~CMannWhitneyTest(void)
  {
  }

//+------------------------------------------------------------------+
//| Copy the first sample into internal storage                      |
//+------------------------------------------------------------------+
void CMannWhitneyTest::SetGroup1(const double &data[])
  {
   m_n1=::ArraySize(data);
   ::ArrayResize(m_group1,m_n1);
   ::ArrayCopy(m_group1,data);
  }

//+------------------------------------------------------------------+
//| Copy the second sample into internal storage                     |
//+------------------------------------------------------------------+
void CMannWhitneyTest::SetGroup2(const double &data[])
  {
   m_n2=::ArraySize(data);
   ::ArrayResize(m_group2,m_n2);
   ::ArrayCopy(m_group2,data);
  }

//+------------------------------------------------------------------+
//| Custom comparison-based sort for SRankItem by value              |
//+------------------------------------------------------------------+
void CMannWhitneyTest::SortRankItems(SRankItem &items[])
  {
   int total=::ArraySize(items);
   for(int i=0; i<total-1; i++)
     {
      int min_idx=i;
      for(int j=i+1; j<total; j++)
        {
         if(items[j].value<items[min_idx].value)
            min_idx=j;
        }
      if(min_idx!=i)
        {
         SRankItem temp=items[i];
         items[i]=items[min_idx];
         items[min_idx]=temp;
        }
     }
  }

//+------------------------------------------------------------------+
//| Pool both groups, rank them with tie averaging, and return       |
//| ranks aligned to the original group order (group1 then group2)   |
//+------------------------------------------------------------------+
void CMannWhitneyTest::RankPooled(double &ranks[])
  {
   int total=m_n1+m_n2;
   SRankItem items[];
   ::ArrayResize(items,total);

   for(int i=0; i<m_n1; i++)
     {
      items[i].value=m_group1[i];
      items[i].index=i;
      items[i].rank=0.0;
     }
   for(int i=0; i<m_n2; i++)
     {
      items[m_n1+i].value=m_group2[i];
      items[m_n1+i].index=m_n1+i;
      items[m_n1+i].rank=0.0;
     }

   SortRankItems(items);

   m_tie_sum=0.0;
   int pos=0;
   while(pos<total)
     {
      int run_end=pos;
      while(run_end+1<total && items[run_end+1].value==items[pos].value)
         run_end++;

      int run_len=run_end-pos+1;
      double avg_rank=(pos+1+run_end+1)/2.0;
      for(int k=pos; k<=run_end; k++)
         items[k].rank=avg_rank;

      if(run_len>1)
        {
         double t=(double)run_len;
         m_tie_sum+=(t*t*t-t);
        }

      pos=run_end+1;
     }

   ::ArrayResize(ranks,total);
   for(int i=0; i<total; i++)
      ranks[items[i].index]=items[i].rank;
  }

//+------------------------------------------------------------------+
//| Standard normal CDF, Abramowitz & Stegun 26.2.17 approximation   |
//+------------------------------------------------------------------+
double CMannWhitneyTest::NormalCDF(double z)
  {
   double sign=1.0;
   if(z<0.0)
     {
      sign=-1.0;
      z=-z;
     }

   double b1= 0.319381530;
   double b2=-0.356563782;
   double b3= 1.781477937;
   double b4=-1.821255978;
   double b5= 1.330274429;
   double p = 0.2316419;
   double c = 0.39894228; // 1 / sqrt(2*pi)

   double t=1.0/(1.0+p*z);
   double poly=t*(b1+t*(b2+t*(b3+t*(b4+t*b5))));
   double pdf=c*::MathExp(-0.5*z*z);
   double cdf_upper_tail=pdf*poly;

   double result=1.0-cdf_upper_tail;
   if(sign<0.0)
      result=1.0-result;

   return(result);
  }

//+------------------------------------------------------------------+
//| Execute the full Mann-Whitney U test                             |
//+------------------------------------------------------------------+
bool CMannWhitneyTest::Run(void)
  {
   if(m_n1<1 || m_n2<1)
      return(false);

   double ranks[];
   RankPooled(ranks);

   double R1=0.0;
   for(int i=0; i<m_n1; i++)
      R1+=ranks[i];

   double R2=0.0;
   for(int i=0; i<m_n2; i++)
      R2+=ranks[m_n1+i];

   m_U1=R1-(m_n1*(m_n1+1))/2.0;
   m_U2=R2-(m_n2*(m_n2+1))/2.0;
   m_U=::MathMin(m_U1,m_U2);

   double n1d=(double)m_n1;
   double n2d=(double)m_n2;
   double N=n1d+n2d;

   double mu_U=n1d*n2d/2.0;
   double sigma2_U=(n1d*n2d/12.0)*((N+1.0)-m_tie_sum/(N*(N-1.0)));

   if(sigma2_U<=0.0)
      return(false);

   double sigma_U=::MathSqrt(sigma2_U);
   m_z_score=(m_U-mu_U)/sigma_U;

   double upper=NormalCDF(::MathAbs(m_z_score));
   m_p_value=2.0*(1.0-upper);

   if(m_p_value>1.0)
      m_p_value=1.0;

   return(true);
  }

//+------------------------------------------------------------------+
//| Accessor: the Mann-Whitney U statistic                           |
//+------------------------------------------------------------------+
double CMannWhitneyTest::GetU(void) const
  {
   return(m_U);
  }

//+------------------------------------------------------------------+
//| Accessor: the z-score of the normal approximation                |
//+------------------------------------------------------------------+
double CMannWhitneyTest::GetZ(void) const
  {
   return(m_z_score);
  }

//+------------------------------------------------------------------+
//| Accessor: the two-tailed p-value                                 |
//+------------------------------------------------------------------+
double CMannWhitneyTest::GetP(void) const
  {
   return(m_p_value);
  }

//+------------------------------------------------------------------+
//| Print a structured summary of the test result                    |
//+------------------------------------------------------------------+
void CMannWhitneyTest::PrintResult(double alpha)
  {
   ::Print("--- Mann-Whitney U Test Result ---");
   ::PrintFormat("Group 1 size (n1): %d",m_n1);
   ::PrintFormat("Group 2 size (n2): %d",m_n2);
   ::PrintFormat("U statistic:       %.4f",m_U);
   ::PrintFormat("z-score:           %.4f",m_z_score);
   ::PrintFormat("p-value:           %.6f",m_p_value);
   ::PrintFormat("Significance level (alpha): %.4f",alpha);

   if(m_p_value<alpha)
     {
      ::PrintFormat("Conclusion: Reject H0 (p=%.6f < alpha=%.4f). "
                    "The two return distributions are statistically distinguishable.",
                    m_p_value,alpha);
     }
   else
     {
      ::PrintFormat("Conclusion: Fail to reject H0 (p=%.6f >= alpha=%.4f). "
                    "No statistically significant distributional difference was detected.",
                    m_p_value,alpha);
     }
  }

//+------------------------------------------------------------------+
//| Compute Q1, median (Q2), and Q3 from an unsorted sample          |
//+------------------------------------------------------------------+
void CMannWhitneyTest::ComputeQuartiles(const double &data[],double &q1,double &q2,double &q3)
  {
   int n=::ArraySize(data);
   double sorted[];
   ::ArrayResize(sorted,n);
   ::ArrayCopy(sorted,data);
   ::ArraySort(sorted);

   if(n==1)
     {
      q1=sorted[0];
      q2=sorted[0];
      q3=sorted[0];
      return;
     }

   if(n%2==0)
      q2=(sorted[n/2-1]+sorted[n/2])/2.0;
   else
      q2=sorted[n/2];

   int lower_count=n/2;
   double lower[];
   ::ArrayResize(lower,lower_count);
   for(int i=0; i<lower_count; i++)
      lower[i]=sorted[i];

   int upper_count=n/2;
   double upper[];
   ::ArrayResize(upper,upper_count);
   int start_upper=(n%2==0)? n/2 : n/2+1;
   for(int i=0; i<upper_count; i++)
      upper[i]=sorted[start_upper+i];

   if(lower_count%2==0)
      q1=(lower[lower_count/2-1]+lower[lower_count/2])/2.0;
   else
      q1=lower[lower_count/2];

   if(upper_count%2==0)
      q3=(upper[upper_count/2-1]+upper[upper_count/2])/2.0;
   else
      q3=upper[upper_count/2];
  }

#endif // MANNWHITNEYTEST_MQH
//+------------------------------------------------------------------+