multinomial-histograms.c

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

/*
Copyright 2020 Dustin Lang.  Released under BSD 3-term license.

Given a multinomial probability distribution, and a total number
of events, generates all histograms, evaluates their probabilities.
Given a "tail" probability level, computes the total probability mass of
histograms with probability less than the tail.
This lets you give something like a false alarm probability or P value -- you would
expect to get an event with this tail probability or more only once in N draws.

Example output:
Total prob: 1
Number of histograms: 888030
Tail probability: 1e-16
Number of tail probabilities: 64
Total probability in tail: 9.78384e-16
*/

const int FACT_TABLE_N = 100;
double fact_table[FACT_TABLE_N];
double factorial(int k) {
    return fact_table[k];
}

double pmf(double* binprobs, int* counts, int k, int n) {
    double p = 1.0;
    int i;
    p *= factorial(n);
    for (i=0; i<k; i++)
        p *= pow(binprobs[i], counts[i]) / factorial(counts[i]);
    return p;
}

// woo globals
double total_prob = 0.0;
int n_total = 0;

int n_tail = 0;
double p_tail = 0.0;

double tail_cutoff = 0;

void assign(int* counts, double* binprobs, int k, int n, int i, int n_left) {
    // Assign "n" total counts to bins i through k-1 (inclusive).

    // assign last bin with remaining n_left
    if (i == k-1) {
        counts[k-1] = n_left;
        double p = pmf(binprobs, counts, k, n);
        total_prob += p;
        n_total++;
        
        if (p <= tail_cutoff) {
            n_tail++;
            p_tail += p;
        }
        
        return;
    }
    int j;
    for (j=0; j<=n_left; j++) {
        counts[i] = j;
        assign(counts, binprobs, k, n, i+1, n_left - j);
    }
}

int main() {

    tail_cutoff = 1e-16;

    const int n=20;
    const int k=8;

    int counts[k];

    double binprobs[k];
    int i;

    // init
    fact_table[0] = 1.0;
    for (i=1; i<FACT_TABLE_N; i++)
        fact_table[i] = fact_table[i-1] * (double)i;

    for (i=0; i<k; i++)
        binprobs[i] = 1./(double)k;

    assign(counts, binprobs, k, n, 0, n);
    printf("Total prob: %g\n", total_prob);
    printf("Number of histograms: %i\n", n_total);
    printf("Tail probability: %g\n", tail_cutoff);
    printf("Number of tail probabilities: %i\n", n_tail);
    printf("Total probability in tail: %g\n", p_tail);
}