defmodule Statistics.Tests.TTest do
  import Statistics
  import Statistics.Math
  alias Statistics.Distributions.T

  @moduledoc """
  Student's t test

  """

  @doc """
  A two-sided test for the null hypothesis that the 
  expected value (mean) of a sample of independent
  observations a is equal to the given population mean, `popmean`.

  Returns the _t_ statistic, and the _p_ value.

  ## Example

      iex> Statistics.Tests.TTest.one_sample([1,2,3,2,1], 3)
      %{p: 0.023206570788795993, t: -3.585685828003181}

  """
  def one_sample(list, popmean) do
    df = length(list) - 1
    t = (mean(list) - popmean) / (stdev(list) / sqrt(length(list)))
    p = get_t_prob(t, df)
    %{t: t, p: p}
  end

  @doc """
  A two-sided test for the null hypothesis that the 
  mean of `list1` is different to the mean of `list2`.

  The variance of the lists should be equal but the 
  sample size of each last can be different.

  Returns the _t_ statistic, and the _p_ value.

  ## Example

      iex> Statistics.Tests.TTest.ind_samples([1,2,3,2,1], [3,2,4,3,5])
      %{p: 0.022802155958137702, t: -2.82842712474619}

      iex> Statistics.Tests.TTest.ind_samples([1,2,3,2,1], [3,2,4,3,5,4,5,6])
      %{p: 0.0044530673387188, t: -3.5858542135407596}

  """
  def ind_samples(list1, list2) do
    df = length(list1) + length(list2) - 2
    mu1 = mean(list1)
    mu2 = mean(list2)
    # calculate pooled standard deviation and
    # sample proportion differently when 
    # sample sizes are unequal
    {sp, sz} =
      case length(list1) == length(list2) do
        true ->
          spt = sqrt((variance(list1) + variance(list2)) / 2)
          szt = sqrt(2 / length(list1))
          {spt, szt}

        false ->
          # weight variances by sample size
          adj_var1 = (length(list1) - 1) * variance(list1)
          adj_var2 = (length(list2) - 1) * variance(list2)
          spf = sqrt((adj_var1 + adj_var2) / df)
          szf = sqrt(1 / length(list1) + 1 / length(list2))
          {spf, szf}
      end

    t = (mu1 - mu2) / (sp * sz)
    p = get_t_prob(t, df)
    %{t: t, p: p}
  end

  defp get_t_prob(t, df) do
    c = T.cdf(df).(t)

    p =
      case t < 0.0 do
        true -> c
        false -> 1 - c
      end

    # two-sided test
    case p < 0.5 do
      true -> 2 * p
      false -> 1.0
    end
  end
end
