+ nom = pow2 (var_i/n_i + var_j/n_j);
+ denom = pow2 (var_i/n_i) / (n_i - 1) + pow2 (var_j/n_j) / (n_j - 1);
+
+ return nom / denom;
+}
+
+static double lsd_pinv (double std_err, double alpha, double df, int k UNUSED, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+ return std_err * gsl_cdf_tdist_Pinv (1.0 - alpha / 2.0, df);
+}
+
+static double bonferroni_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+ const int m = k * (k - 1) / 2;
+ return std_err * gsl_cdf_tdist_Pinv (1.0 - alpha / (2.0 * m), df);
+}
+
+static double sidak_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i, const struct moments1 *mom_j)
+{
+ const double m = k * (k - 1) / 2;
+ double lp = 1.0 - exp (log (1.0 - alpha) / m ) ;
+ return std_err * gsl_cdf_tdist_Pinv (1.0 - lp / 2.0, df);
+}
+
+static double tukey_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+ return std_err / sqrt (2.0) * qtukey (1 - alpha, 1.0, k, df, 1, 0);
+}
+
+static double scheffe_pinv (double std_err, double alpha, double df, int k, const struct moments1 *mom_i UNUSED, const struct moments1 *mom_j UNUSED)
+{
+ double x = (k - 1) * gsl_cdf_fdist_Pinv (1.0 - alpha, k - 1, df);
+ return std_err * sqrt (x);
+}
+
+static double gh_pinv (double std_err UNUSED, double alpha, double df, int k, const struct moments1 *mom_i, const struct moments1 *mom_j)
+{
+ double n_i, mean_i, var_i;
+ double n_j, mean_j, var_j;
+ double m;
+
+ moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);
+ moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+ m = sqrt ((var_i/n_i + var_j/n_j) / 2.0);
+
+ return m * qtukey (1 - alpha, 1.0, k, df, 1, 0);
+}
+
+
+static double
+multiple_comparison_sig (double std_err,
+ const struct per_var_ws *pvw,
+ const struct descriptive_data *dd_i, const struct descriptive_data *dd_j,
+ const struct posthoc *ph)
+{
+ int k = pvw->n_groups;
+ double df = ph->dff (pvw, dd_i->mom, dd_j->mom);
+ double ts = ph->tsf (k, dd_i->mom, dd_j->mom, std_err);
+ return ph->p1f (ts, k - 1, df);
+}
+
+static double
+mc_half_range (const struct oneway_spec *cmd, const struct per_var_ws *pvw, double std_err, const struct descriptive_data *dd_i, const struct descriptive_data *dd_j, const struct posthoc *ph)
+{
+ int k = pvw->n_groups;
+ double df = ph->dff (pvw, dd_i->mom, dd_j->mom);
+
+ return ph->pinv (std_err, cmd->alpha, df, k, dd_i->mom, dd_j->mom);
+}
+
+static double tukey_1tailsig (double ts, double df1, double df2)
+{
+ double twotailedsig = 1.0 - ptukey (ts, 1.0, df1 + 1, df2, 1, 0);
+
+ return twotailedsig / 2.0;
+}
+
+static double lsd_1tailsig (double ts, double df1 UNUSED, double df2)
+{
+ return ts < 0 ? gsl_cdf_tdist_P (ts, df2) : gsl_cdf_tdist_Q (ts, df2);
+}
+
+static double sidak_1tailsig (double ts, double df1, double df2)
+{
+ double ex = (df1 + 1.0) * df1 / 2.0;
+ double lsd_sig = 2 * lsd_1tailsig (ts, df1, df2);
+
+ return 0.5 * (1.0 - pow (1.0 - lsd_sig, ex));
+}
+
+static double bonferroni_1tailsig (double ts, double df1, double df2)
+{
+ const int m = (df1 + 1) * df1 / 2;
+
+ double p = ts < 0 ? gsl_cdf_tdist_P (ts, df2) : gsl_cdf_tdist_Q (ts, df2);
+ p *= m;
+
+ return p > 0.5 ? 0.5 : p;
+}
+
+static double scheffe_1tailsig (double ts, double df1, double df2)
+{
+ return 0.5 * gsl_cdf_fdist_Q (ts, df1, df2);
+}
+
+
+static double tukey_test_stat (int k, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+ double n_i, mean_i, var_i;
+ double n_j, mean_j, var_j;
+
+ moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);
+ moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+ double ts = (mean_i - mean_j) / std_err;
+ ts = fabs (ts) * sqrt (2.0);
+
+ return ts;
+}
+
+static double lsd_test_stat (int k, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+ double n_i, mean_i, var_i;
+ double n_j, mean_j, var_j;
+
+ moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);
+ moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+ return (mean_i - mean_j) / std_err;
+}
+
+static double scheffe_test_stat (int k, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+ double n_i, mean_i, var_i;
+ double n_j, mean_j, var_j;
+
+ moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);
+ moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+ double t = (mean_i - mean_j) / std_err;
+ t = pow2 (t);
+ t /= k - 1;
+
+ return t;
+}
+
+static double gh_test_stat (int k, const struct moments1 *mom_i, const struct moments1 *mom_j, double std_err)
+{
+ double n_i, mean_i, var_i;
+ double n_j, mean_j, var_j;
+
+ moments1_calculate (mom_i, &n_i, &mean_i, &var_i, 0, 0);
+ moments1_calculate (mom_j, &n_j, &mean_j, &var_j, 0, 0);
+
+ double thing = var_i / n_i + var_j / n_j;
+ thing /= 2.0;
+ thing = sqrt (thing);
+
+ double ts = (mean_i - mean_j) / thing;
+
+ return fabs (ts);
+}
+
+
+
+static const struct posthoc ph_tests [] =
+ {
+ { "LSD", N_("LSD"), df_common, lsd_test_stat, lsd_1tailsig, lsd_pinv},
+ { "TUKEY", N_("Tukey HSD"), df_common, tukey_test_stat, tukey_1tailsig, tukey_pinv},
+ { "BONFERRONI", N_("Bonferroni"), df_common, lsd_test_stat, bonferroni_1tailsig, bonferroni_pinv},
+ { "SCHEFFE", N_("Scheffé"), df_common, scheffe_test_stat, scheffe_1tailsig, scheffe_pinv},
+ { "GH", N_("Games-Howell"), df_individual, gh_test_stat, tukey_1tailsig, gh_pinv},
+ { "SIDAK", N_("Šidák"), df_common, lsd_test_stat, sidak_1tailsig, sidak_pinv}
+ };