logp_stats.c

#include <stdio.h>
#include <math.h>
#include "logp_stats.h"

void
stats_init(statistics *stats)
{
    stats->num = 0;
    stats->sum = 0.0;
    stats->sum_sq = 0.0;
    stats->min = 9e99;
    stats->max = 0.0;
}

void
stats_update(statistics *stats, double val)
{
    stats->num++;
    stats->sum += val;
    stats->sum_sq += val * val;
    if (val > stats->max) {
	stats->max = val;
    }
    if (val < stats->min) {
	stats->min = val;
    }
}

double
stats_min(statistics *stats)
{
    return stats->min;
}

double
stats_max(statistics *stats)
{
    return stats->max;
}

double
stats_mean(statistics *stats)
{
    if (stats->num > 0) {
	return stats->sum / stats->num;
    } else {
	return 0.0;
    }
}

double
stats_var(statistics *stats)
{
    if (stats->num > 1) {
	double var;

	var = (stats->sum_sq - (stats->sum * stats->sum) / stats->num) / 
	       (stats->num - 1);
	/* for VAR very close to 0, due to rounding, the difference may
	 * in fact become slightly negative, so:
	 */
	if (var >= 0.0) {
	    return var;
	} else {
	    return 0.0;
	}
    } else {
	return 0.0;
    }
}

/* Following table for the t distribution with n degrees of freedom is
 * derived from a table in TeX code by:
 *   Robert J. MacG. Dawson
 *   Dept.of Math & Computing Science
 *   St. Mary's University, Halifax, NS, Canada B3H 3C3
 */
#define STATS_N_DF 19
static int t_tab_df[STATS_N_DF] = 
    { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12,
      15, 20, 25, 30, 40, 50, 100, 1000 };

static struct {
    double t_1tail;
    double t_2tail;
    double t_val[STATS_N_DF];
} t_tab[] = {
    { 0.5, 1,
      { 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
	0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 }},
    { 0.4, 0.8,
      { 0.325, 0.289, 0.277, 0.271, 0.267, 0.265, 0.263, 0.262, 0.261, 0.260,
	0.259, 0.258, 0.257, 0.256, 0.256, 0.255, 0.255, 0.254, 0.253 }},
    { 0.3, 0.6,
      { 0.727, 0.617, 0.584, 0.569, 0.559, 0.553, 0.549, 0.546, 0.543, 0.542,
	0.539, 0.536, 0.533, 0.531, 0.530, 0.529, 0.528, 0.526, 0.525 }},
    { 0.25, 0.5,
      { 1.000, 0.816, 0.765, 0.741, 0.727, 0.718, 0.711, 0.706, 0.703, 0.700,
	0.695, 0.691, 0.687, 0.684, 0.683, 0.681, 0.679, 0.677, 0.675 }},
    { 0.2, 0.4,
      { 1.376, 1.061, 0.978, 0.941, 0.920, 0.906, 0.896, 0.889, 0.883, 0.879,
	0.873, 0.866, 0.860, 0.856, 0.854, 0.851, 0.849, 0.845, 0.842 }},
    { 0.17, 0.34,
      { 1.691, 1.242, 1.132, 1.082, 1.054, 1.036, 1.024, 1.015, 1.008, 1.002,
	0.994, 0.986, 0.977, 0.973, 0.970, 0.966, 0.963, 0.959, 0.955 }},
    { 0.15, 0.3,
      { 1.963, 1.386, 1.250, 1.190, 1.156, 1.134, 1.119, 1.108, 1.100, 1.093,
	1.083, 1.074, 1.064, 1.058, 1.055, 1.050, 1.047, 1.042, 1.037 }},
    { 0.14, 0.28,
      { 2.125, 1.467, 1.315, 1.248, 1.211, 1.187, 1.171, 1.159, 1.149, 1.142,
	1.131, 1.121, 1.110, 1.104, 1.100, 1.095, 1.092, 1.086, 1.081 }},
    { 0.13, 0.26,
      { 2.311, 1.556, 1.385, 1.311, 1.270, 1.244, 1.226, 1.212, 1.202, 1.194,
	1.182, 1.171, 1.159, 1.153, 1.148, 1.143, 1.139, 1.133, 1.127 }},
    { 0.12, 0.24,
      { 2.526, 1.654, 1.462, 1.379, 1.333, 1.304, 1.284, 1.269, 1.258, 1.249,
	1.236, 1.224, 1.211, 1.204, 1.199, 1.193, 1.189, 1.182, 1.176 }},
    { 0.11, 0.22,
      { 2.778, 1.763, 1.545, 1.453, 1.401, 1.369, 1.347, 1.331, 1.318, 1.308,
	1.294, 1.280, 1.266, 1.258, 1.253, 1.246, 1.242, 1.234, 1.227 }},
    { 0.1, 0.2,
      { 3.078, 1.886, 1.638, 1.533, 1.476, 1.440, 1.415, 1.397, 1.383, 1.372,
	1.356, 1.341, 1.325, 1.316, 1.310, 1.303, 1.299, 1.290, 1.282 }},
    { 0.09, 0.18,
      { 3.442, 2.026, 1.741, 1.623, 1.558, 1.517, 1.489, 1.469, 1.454, 1.442,
	1.424, 1.406, 1.389, 1.379, 1.373, 1.365, 1.360, 1.350, 1.342 }},
    { 0.08, 0.16,
      { 3.895, 2.189, 1.859, 1.723, 1.649, 1.603, 1.572, 1.549, 1.532, 1.518,
	1.498, 1.478, 1.459, 1.448, 1.441, 1.432, 1.426, 1.416, 1.406 }},
    { 0.075, 0.15,
      { 4.165, 2.282, 1.924, 1.778, 1.699, 1.650, 1.617, 1.592, 1.574, 1.559,
	1.538, 1.517, 1.497, 1.485, 1.477, 1.468, 1.462, 1.451, 1.441 }},
    { 0.07, 0.14,
      { 4.474, 2.383, 1.995, 1.838, 1.753, 1.700, 1.664, 1.638, 1.619, 1.603,
	1.580, 1.558, 1.537, 1.524, 1.516, 1.506, 1.500, 1.488, 1.477 }},
    { 0.065, 0.13,
      { 4.829, 2.495, 2.072, 1.902, 1.810, 1.754, 1.715, 1.687, 1.666, 1.650,
	1.626, 1.602, 1.579, 1.566, 1.557, 1.546, 1.539, 1.527, 1.515 }},
    { 0.06, 0.12,
      { 5.242, 2.620, 2.156, 1.971, 1.873, 1.812, 1.770, 1.740, 1.718, 1.700,
	1.674, 1.649, 1.624, 1.610, 1.600, 1.589, 1.582, 1.568, 1.556 }},
    { 0.055, 0.11,
      { 5.730, 2.760, 2.249, 2.048, 1.941, 1.874, 1.830, 1.797, 1.773, 1.754,
	1.726, 1.699, 1.672, 1.657, 1.647, 1.635, 1.627, 1.613, 1.600 }},
    { 0.05, 0.1,
      { 6.314, 2.920, 2.353, 2.132, 2.015, 1.943, 1.895, 1.860, 1.833, 1.812,
	1.782, 1.753, 1.725, 1.708, 1.697, 1.684, 1.676, 1.660, 1.646 }},
    { 0.045, 0.09,
      { 7.026, 3.104, 2.471, 2.226, 2.098, 2.019, 1.966, 1.928, 1.899, 1.877,
	1.844, 1.812, 1.782, 1.764, 1.752, 1.737, 1.729, 1.712, 1.697 }},
    { 0.04, 0.08,
      { 7.916, 3.320, 2.605, 2.333, 2.191, 2.104, 2.046, 2.004, 1.973, 1.948,
	1.912, 1.878, 1.844, 1.825, 1.812, 1.796, 1.787, 1.769, 1.752 }},
    { 0.035, 0.07,
      { 9.058, 3.578, 2.763, 2.456, 2.297, 2.201, 2.136, 2.090, 2.055, 2.028,
	1.989, 1.951, 1.914, 1.893, 1.879, 1.862, 1.852, 1.832, 1.814 }},
    { 0.03, 0.06,
      { 10.57, 3.896, 2.951, 2.601, 2.422, 2.313, 2.241, 2.189, 2.150, 2.120,
	2.076, 2.034, 1.994, 1.970, 1.955, 1.936, 1.924, 1.902, 1.883 }},
    { 0.025, 0.05,
      { 12.70, 4.303, 3.182, 2.776, 2.571, 2.447, 2.365, 2.306, 2.262, 2.228,
	2.179, 2.131, 2.086, 2.060, 2.042, 2.021, 2.009, 1.984, 1.962 }},
    { 0.02, 0.04,
      { 15.89, 4.849, 3.482, 2.999, 2.757, 2.612, 2.517, 2.449, 2.398, 2.359,
	2.303, 2.249, 2.197, 2.167, 2.147, 2.123, 2.109, 2.081, 2.056 }},
    { 0.017, 0.034,
      { 18.71, 5.284, 3.712, 3.166, 2.895, 2.734, 2.628, 2.553, 2.498, 2.454,
	2.392, 2.333, 2.276, 2.243, 2.222, 2.195, 2.180, 2.150, 2.123 }},
    { 0.015, 0.03,
      { 21.21, 5.643, 3.896, 3.298, 3.003, 2.829, 2.715, 2.634, 2.574, 2.527,
	2.461, 2.397, 2.336, 2.301, 2.278, 2.250, 2.234, 2.202, 2.173 }},
    { 0.012, 0.024,
      { 26.51, 6.338, 4.241, 3.541, 3.200, 3.000, 2.870, 2.778, 2.710, 2.658,
	2.582, 2.511, 2.442, 2.403, 2.378, 2.346, 2.328, 2.292, 2.261 }},
    { 0.01, 0.02,
      { 31.82, 6.965, 4.541, 3.747, 3.365, 3.143, 2.998, 2.896, 2.821, 2.764,
	2.681, 2.602, 2.528, 2.485, 2.457, 2.423, 2.403, 2.364, 2.330 }}, 
    { 0.007, 0.014,
      { 45.47, 8.363, 5.175, 4.173, 3.700, 3.428, 3.253, 3.131, 3.041, 2.972,
	2.873, 2.781, 2.693, 2.642, 2.610, 2.570, 2.547, 2.501, 2.462 }},
    { 0.005, 0.01,
      { 63.65, 9.925, 5.841, 4.604, 4.032, 3.707, 3.499, 3.355, 3.250, 3.169,
	3.055, 2.947, 2.845, 2.787, 2.750, 2.704, 2.678, 2.626, 2.581 }},
    { 0.003, 0.006,
      { 106.1, 12.85, 6.994, 5.321, 4.570, 4.152, 3.887, 3.705, 3.573, 3.472,
	3.330, 3.197, 3.073, 3.003, 2.957, 2.902, 2.870, 2.808, 2.754 }},
    { 0.002, 0.004,
      { 159.2, 15.76, 8.053, 5.951, 5.030, 4.524, 4.207, 3.991, 3.835, 3.716,
	3.550, 3.395, 3.251, 3.170, 3.118, 3.055, 3.018, 2.946, 2.885 }}, 
    { 0.001, 0.002,
      { 318.3, 22.32, 10.21, 7.173, 5.893, 5.208, 4.785, 4.501, 4.297, 4.144,
	3.930, 3.733, 3.552, 3.450, 3.385, 3.307, 3.261, 3.174, 3.098 }},
    { .0005, 0.001,
      { 636.6, 31.59, 12.92, 8.610, 6.869, 5.959, 5.408, 5.041, 4.781, 4.587,
	4.318, 4.073, 3.850, 3.725, 3.646, 3.551, 3.496, 3.391, 3.300 }},
    { .0001, .0002,
      { 3183 , 70.70, 22.20, 13.03, 9.678, 8.025, 7.063, 6.442, 6.010, 5.694,
	5.263, 4.880, 4.539, 4.352, 4.234, 4.094, 4.014, 3.862, 3.733 }},
    { 5e-5 , .0001,
      { 6366 , 99.98, 27.99, 15.54, 11.18, 9.082, 7.884, 7.120, 6.594, 6.210,
	5.694, 5.239, 4.837, 4.620, 4.482, 4.320, 4.228, 4.053, 3.906 }}
};

double
stats_t(int df, int tail, double conf)
{
    int max_df_ind, max_conf_ind;
    int df_i, conf_i;

    max_df_ind = sizeof(t_tab_df) / sizeof(t_tab_df[0]);
    max_conf_ind = sizeof(t_tab) / sizeof(t_tab[0]);

    /* find a df that is close;
     * TODO: actually find closest one using binary search.
     */
    for (df_i = 0; df_i < max_df_ind; df_i++) {
	if (t_tab_df[df_i] >= df) {
	    break;
	}
    }
    if (df_i >= max_df_ind) {
	df_i = max_df_ind - 1;
    }

    /* Find the requested tail probability that is close.
     */
    for (conf_i = max_conf_ind - 1; conf_i >= 0; conf_i--) {
	if (tail == 1) {
	    if (t_tab[conf_i].t_1tail >= conf) {
		break;
	    }
	} else {
	    if (t_tab[conf_i].t_2tail >= conf) {
		break;
	    }
	}
    }
    if (conf_i < 0) {
	conf_i = 0;
    }

#if 0
    printf("stats_t(%d, %.4f) = %.4f (df = %d, conf = %.4f)\n",
	   df, conf, t_tab[conf_i].t_val[df_i],
	   t_tab_df[df_i], t_tab[conf_i].t_1tail);
#endif
	   
    return t_tab[conf_i].t_val[df_i];
}

double
stats_conf_int(statistics *stats, int tail, int iters, double confint)
{
#if 0
    printf("iters=%d stats_t=%.6f var=%.6f sqrt(var/iters)=%.6f conf=%.6f\n",
	   iters,
           stats_t(tail, iters - 1, 1.0 - confint),
	   stats_var(stats),
	   sqrt(stats_var(stats) / iters),
	   stats_t(tail, iters - 1, 1.0 - confint) *
		sqrt(stats_var(stats) / iters));
#endif

    return stats_t(tail, iters - 1, 1.0 - confint) *
	sqrt(stats_var(stats) / iters);
}