-module(eministat_analysis).

-include("eministat.hrl").

-export([outlier_variance/3]).
-export([relative/3]).

%% -- OUTLIER VARIANCE ------
%%
%% @doc outlier_variance/3 computes the severity of the outlier variance
%%
%% http://www.ellipticgroup.com/misc/article_supplement.pdf
%%
%% @end
outlier_variance(_, Sigma, _) when Sigma < 0.000000000000001 -> {0.0, unaffected};
outlier_variance(Mu, Sigma, A) ->
    MinBy = fun(F, Q, R) -> min(F(Q), F(R)) end,

    MuA = Mu / A,
    MugMin = MuA / 2,
    SigmaG = min(MugMin / 4, Sigma / math:sqrt(A)),
    SigmaG2 = SigmaG * SigmaG,
    Sigma2 = Sigma * Sigma,
    
    VarOut = fun(C) ->
        AC = A - C,
        (AC / A) * (Sigma2 - AC * SigmaG2)
    end,
    
    CMax = fun(X) ->
        K = MuA - X,
        D = K * K,
        AD = A * D,
        K0 = -A * AD,
        K1 = Sigma2 - A * SigmaG2 + AD,
        Det = K1 * K1 - 4 * SigmaG2 * K0,
        trunc(-2 * K0 / (K1 + math:sqrt(Det)))
    end,

    VarOutMin = MinBy(VarOut, 1, (MinBy(CMax, 0, MugMin))) / Sigma2,
    case VarOutMin of
        K when K < 0.01 -> {K, unaffected};
        K when K < 0.1 -> {K, slight};
        K when K < 0.5 -> {K, moderate};
        K when K > 0.5 -> {K, severe}
    end.

%% -- RELATIVE #dataset{} COMPARISONS -----


relative(#dataset { n = DsN } = Ds, #dataset { n = RsN } = Rs, ConfIdx) ->
    I = DsN + RsN - 2,
    T = element(ConfIdx, student_lookup(I)),

    Spool1 = (DsN - 1) * eministat_ds:variance(Ds) + (RsN - 1) * eministat_ds:variance(Rs),
    Spool = math:sqrt(Spool1 / I),

    S = Spool * math:sqrt(1.0 / DsN + 1.0 / RsN),
    D = eministat_ds:mean(Ds) - eministat_ds:mean(Rs),
    E = T * S,

    case abs(D) > E of
        false ->
            {no_difference, element(ConfIdx, student_pct())};
        true ->
            {difference,
             #{ confidence_level => element(ConfIdx, student_pct()),
                difference => {D, E},
                difference_pct => {
                  D * 100 / eministat_ds:mean(Rs),
                  E * 100 / eministat_ds:mean(Rs)},
                pooled_s => Spool
              }}
    end.

%% -- STUDENT's T TABLES -----

%% Constant tables, represented as tuples for O(1) lookup speeds
student_pct() -> {80.0, 90.0, 95.0, 98.0, 99.0, 99.5}.

student_inf() ->
	{	1.282,	1.645,	1.960,	2.326,	2.576,	3.090  }. %%  inf

-define(NSTUDENT, 100). %% Number of elements in the students distribution lookup table
student() ->
  {
	{	3.078,	6.314,	12.706,	31.821,	63.657,	318.313  }, %%  1.
	{	1.886,	2.920,	4.303,	6.965,	9.925,	22.327  }, %%  2.
	{	1.638,	2.353,	3.182,	4.541,	5.841,	10.215  }, %%  3.
	{	1.533,	2.132,	2.776,	3.747,	4.604,	7.173  }, %%  4.
	{	1.476,	2.015,	2.571,	3.365,	4.032,	5.893  }, %%  5.
	{	1.440,	1.943,	2.447,	3.143,	3.707,	5.208  }, %%  6.
	{	1.415,	1.895,	2.365,	2.998,	3.499,	4.782  }, %%  7.
	{	1.397,	1.860,	2.306,	2.896,	3.355,	4.499  }, %%  8.
	{	1.383,	1.833,	2.262,	2.821,	3.250,	4.296  }, %%  9.
	{	1.372,	1.812,	2.228,	2.764,	3.169,	4.143  }, %%  10.
	{	1.363,	1.796,	2.201,	2.718,	3.106,	4.024  }, %%  11.
	{	1.356,	1.782,	2.179,	2.681,	3.055,	3.929  }, %%  12.
	{	1.350,	1.771,	2.160,	2.650,	3.012,	3.852  }, %%  13.
	{	1.345,	1.761,	2.145,	2.624,	2.977,	3.787  }, %%  14.
	{	1.341,	1.753,	2.131,	2.602,	2.947,	3.733  }, %%  15.
	{	1.337,	1.746,	2.120,	2.583,	2.921,	3.686  }, %%  16.
	{	1.333,	1.740,	2.110,	2.567,	2.898,	3.646  }, %%  17.
	{	1.330,	1.734,	2.101,	2.552,	2.878,	3.610  }, %%  18.
	{	1.328,	1.729,	2.093,	2.539,	2.861,	3.579  }, %%  19.
	{	1.325,	1.725,	2.086,	2.528,	2.845,	3.552  }, %%  20.
	{	1.323,	1.721,	2.080,	2.518,	2.831,	3.527  }, %%  21.
	{	1.321,	1.717,	2.074,	2.508,	2.819,	3.505  }, %%  22.
	{	1.319,	1.714,	2.069,	2.500,	2.807,	3.485  }, %%  23.
	{	1.318,	1.711,	2.064,	2.492,	2.797,	3.467  }, %%  24.
	{	1.316,	1.708,	2.060,	2.485,	2.787,	3.450  }, %%  25.
	{	1.315,	1.706,	2.056,	2.479,	2.779,	3.435  }, %%  26.
	{	1.314,	1.703,	2.052,	2.473,	2.771,	3.421  }, %%  27.
	{	1.313,	1.701,	2.048,	2.467,	2.763,	3.408  }, %%  28.
	{	1.311,	1.699,	2.045,	2.462,	2.756,	3.396  }, %%  29.
	{	1.310,	1.697,	2.042,	2.457,	2.750,	3.385  }, %%  30.
	{	1.309,	1.696,	2.040,	2.453,	2.744,	3.375  }, %%  31.
	{	1.309,	1.694,	2.037,	2.449,	2.738,	3.365  }, %%  32.
	{	1.308,	1.692,	2.035,	2.445,	2.733,	3.356  }, %%  33.
	{	1.307,	1.691,	2.032,	2.441,	2.728,	3.348  }, %%  34.
	{	1.306,	1.690,	2.030,	2.438,	2.724,	3.340  }, %%  35.
	{	1.306,	1.688,	2.028,	2.434,	2.719,	3.333  }, %%  36.
	{	1.305,	1.687,	2.026,	2.431,	2.715,	3.326  }, %%  37.
	{	1.304,	1.686,	2.024,	2.429,	2.712,	3.319  }, %%  38.
	{	1.304,	1.685,	2.023,	2.426,	2.708,	3.313  }, %%  39.
	{	1.303,	1.684,	2.021,	2.423,	2.704,	3.307  }, %%  40.
	{	1.303,	1.683,	2.020,	2.421,	2.701,	3.301  }, %%  41.
	{	1.302,	1.682,	2.018,	2.418,	2.698,	3.296  }, %%  42.
	{	1.302,	1.681,	2.017,	2.416,	2.695,	3.291  }, %%  43.
	{	1.301,	1.680,	2.015,	2.414,	2.692,	3.286  }, %%  44.
	{	1.301,	1.679,	2.014,	2.412,	2.690,	3.281  }, %%  45.
	{	1.300,	1.679,	2.013,	2.410,	2.687,	3.277  }, %%  46.
	{	1.300,	1.678,	2.012,	2.408,	2.685,	3.273  }, %%  47.
	{	1.299,	1.677,	2.011,	2.407,	2.682,	3.269  }, %%  48.
	{	1.299,	1.677,	2.010,	2.405,	2.680,	3.265  }, %%  49.
	{	1.299,	1.676,	2.009,	2.403,	2.678,	3.261  }, %%  50.
	{	1.298,	1.675,	2.008,	2.402,	2.676,	3.258  }, %%  51.
	{	1.298,	1.675,	2.007,	2.400,	2.674,	3.255  }, %%  52.
	{	1.298,	1.674,	2.006,	2.399,	2.672,	3.251  }, %%  53.
	{	1.297,	1.674,	2.005,	2.397,	2.670,	3.248  }, %%  54.
	{	1.297,	1.673,	2.004,	2.396,	2.668,	3.245  }, %%  55.
	{	1.297,	1.673,	2.003,	2.395,	2.667,	3.242  }, %%  56.
	{	1.297,	1.672,	2.002,	2.394,	2.665,	3.239  }, %%  57.
	{	1.296,	1.672,	2.002,	2.392,	2.663,	3.237  }, %%  58.
	{	1.296,	1.671,	2.001,	2.391,	2.662,	3.234  }, %%  59.
	{	1.296,	1.671,	2.000,	2.390,	2.660,	3.232  }, %%  60.
	{	1.296,	1.670,	2.000,	2.389,	2.659,	3.229  }, %%  61.
	{	1.295,	1.670,	1.999,	2.388,	2.657,	3.227  }, %%  62.
	{	1.295,	1.669,	1.998,	2.387,	2.656,	3.225  }, %%  63.
	{	1.295,	1.669,	1.998,	2.386,	2.655,	3.223  }, %%  64.
	{	1.295,	1.669,	1.997,	2.385,	2.654,	3.220  }, %%  65.
	{	1.295,	1.668,	1.997,	2.384,	2.652,	3.218  }, %%  66.
	{	1.294,	1.668,	1.996,	2.383,	2.651,	3.216  }, %%  67.
	{	1.294,	1.668,	1.995,	2.382,	2.650,	3.214  }, %%  68.
	{	1.294,	1.667,	1.995,	2.382,	2.649,	3.213  }, %%  69.
	{	1.294,	1.667,	1.994,	2.381,	2.648,	3.211  }, %%  70.
	{	1.294,	1.667,	1.994,	2.380,	2.647,	3.209  }, %%  71.
	{	1.293,	1.666,	1.993,	2.379,	2.646,	3.207  }, %%  72.
	{	1.293,	1.666,	1.993,	2.379,	2.645,	3.206  }, %%  73.
	{	1.293,	1.666,	1.993,	2.378,	2.644,	3.204  }, %%  74.
	{	1.293,	1.665,	1.992,	2.377,	2.643,	3.202  }, %%  75.
	{	1.293,	1.665,	1.992,	2.376,	2.642,	3.201  }, %%  76.
	{	1.293,	1.665,	1.991,	2.376,	2.641,	3.199  }, %%  77.
	{	1.292,	1.665,	1.991,	2.375,	2.640,	3.198  }, %%  78.
	{	1.292,	1.664,	1.990,	2.374,	2.640,	3.197  }, %%  79.
	{	1.292,	1.664,	1.990,	2.374,	2.639,	3.195  }, %%  80.
	{	1.292,	1.664,	1.990,	2.373,	2.638,	3.194  }, %%  81.
	{	1.292,	1.664,	1.989,	2.373,	2.637,	3.193  }, %%  82.
	{	1.292,	1.663,	1.989,	2.372,	2.636,	3.191  }, %%  83.
	{	1.292,	1.663,	1.989,	2.372,	2.636,	3.190  }, %%  84.
	{	1.292,	1.663,	1.988,	2.371,	2.635,	3.189  }, %%  85.
	{	1.291,	1.663,	1.988,	2.370,	2.634,	3.188  }, %%  86.
	{	1.291,	1.663,	1.988,	2.370,	2.634,	3.187  }, %%  87.
	{	1.291,	1.662,	1.987,	2.369,	2.633,	3.185  }, %%  88.
	{	1.291,	1.662,	1.987,	2.369,	2.632,	3.184  }, %%  89.
	{	1.291,	1.662,	1.987,	2.368,	2.632,	3.183  }, %%  90.
	{	1.291,	1.662,	1.986,	2.368,	2.631,	3.182  }, %%  91.
	{	1.291,	1.662,	1.986,	2.368,	2.630,	3.181  }, %%  92.
	{	1.291,	1.661,	1.986,	2.367,	2.630,	3.180  }, %%  93.
	{	1.291,	1.661,	1.986,	2.367,	2.629,	3.179  }, %%  94.
	{	1.291,	1.661,	1.985,	2.366,	2.629,	3.178  }, %%  95.
	{	1.290,	1.661,	1.985,	2.366,	2.628,	3.177  }, %%  96.
	{	1.290,	1.661,	1.985,	2.365,	2.627,	3.176  }, %%  97.
	{	1.290,	1.661,	1.984,	2.365,	2.627,	3.175  }, %%  98.
	{	1.290,	1.660,	1.984,	2.365,	2.626,	3.175  }, %%  99.
	{	1.290,	1.660,	1.984,	2.364,	2.626,	3.174  } %%  100.
  }.

student_lookup(I) when I > ?NSTUDENT -> student_inf();
student_lookup(I) -> element(I, student()).