approximating some distributions in C

distributions.c

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

// generate uniformly distributed random number in [0, 1]
float rand_float() {
    return (float)rand() / RAND_MAX;
}

const int NORM_APPROX_STAGES = 12;

// generate normally distributed random number with mean 0 and variance 1,
// using sum of uniform random numbers
float rand_norm() {
    float sum = 0;
    for (int i = 0; i < NORM_APPROX_STAGES; i++) {
        sum += rand_float();
    }
    return sum - (float)(NORM_APPROX_STAGES / 2);
}

// generate normally distributed random number with given mean and variance
float rand_norm_var(float mean, float var) {
    return rand_norm() * var + mean;
}

// generate skew-normal random number with given mean, variance, and skew
float rand_skew_norm(float mean, float var, float skew) {
    float x = rand_norm_var(0, 1);
    float y = rand_norm_var(0, 1);
    float z = (skew * x) + (sqrtf(1 - skew*skew)*y);
    if (x < 0) { z *= -1; }
    return z * var + mean;
}

// generate random number from binomial distribution with given number of trials and probability of success
int rand_binom(int trials, float prob) {
    int count = 0;
    for (int i = 0; i < trials; i++) {
        if (rand_float() < prob) {
            count++;
        }
    }
    return count;
}

// generate random number from Poisson distribution with given mean
int rand_poisson(float mean) {
    const float l = expf(-mean);
    int k = 0;
    float p = 1;
    do {
        k++;
        p *= rand_float();
    } while (p > l);
    return k - 1;
}

// generate random number from geometric distribution with given probability of success
int rand_geom(float prob) {
    return rand_poisson(-logf(1 - prob));
}

// generate random number from gamma distribution with given shape and scale
float rand_gamma(float shape, float scale) {
    float sum = 0;
    for (int i = 0; i < shape; i++) {
        sum += -logf(rand_float());
    }
    return sum * scale;
}

// generate random number from beta distribution with given shape parameters
float rand_beta(float a, float b) {
    float x = rand_gamma(a, 1);
    float y = rand_gamma(b, 1);
    return x / (x + y);
}

/// nah, takes too long
//// give an array of values, return a copy sorted from smallest to largest
//float* sort_values(int n, float* values) {
//    float* sorted = malloc(sizeof(float) * n);
//    for (int i = 0; i < n; i++) {
//        sorted[i] = values[i];
//    }
//    for (int i = 0; i < n - 1; i++) {
//        int min = i;
//        for (int j = i; j < n; j++) {
//            if (sorted[j] < sorted[min]) {
//                min = j;
//            }
//        }
//        float temp = sorted[i];
//        sorted[i] = sorted[min];
//        sorted[min] = temp;
//    }
//    return sorted;
//}

// find the min and max values in a given array of floats
void find_min_max(int n, float* values, float* min, float* max) {
    *min = values[0];
    *max = values[0];
    for (int i = 1; i < n; i++) {
        if (values[i] < *min) {
            *min = values[i];
        } else if (values[i] > *max) {
            *max = values[i];
        }
    }
}

// given an array of test output values and a filename, write the values to a .csv file
void write_values(const char* filename, int n, float* values) {
    FILE* file = fopen(filename, "w");
    for (int i = 0; i < n; i++) {
        fprintf(file, "%f\n", values[i]);
    }
    fclose(file);
}

// given an array of test output values and a filename, compute the probability distribution and write it to the file
void write_distribution(const char* filename, int n, float* values) {
    const int BINS = 100;
    int counts[BINS] = {0};

    float min;
    float max;
    find_min_max(n, values, &min, &max);

    for (int i = 0; i < n; i++) {
        int bin = (int)((values[i] - min) / (max - min) * BINS);
        if (bin < 0) {
            bin = 0;
        } else if (bin >= BINS) {
            bin = BINS - 1;
        }
        counts[bin]++;
    }

    FILE* file = fopen(filename, "w");
    for (int i = 0; i < BINS; i++) {
        fprintf(file, "%d\n", counts[i]);
    }
    fclose(file);
}

///////////////////////////////////////////////
///////////////////////////////////////////////
/// tests

// test normal distribution
void test_norm() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_norm();
    }
    write_values("norm.csv", N, values);
}

// test skew normal with several different skew values
void test_skew_norm() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_skew_norm(0, 1, -0.9);
    }
    write_values("skew_norm_n0.9.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_skew_norm(0, 1, 0);
    }
    write_values("skew_norm_0.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_skew_norm(0, 1, 0.9);
    }
    write_values("skew_norm_0.9.csv", N, values);
}

// test binomial with several different probabilities
void test_binom() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_binom(10, 0.5);
    }
    write_values("binom_10_0.5.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_binom(10, 0.1);
    }
    write_values("binom_10_0.1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_binom(10, 0.9);
    }
    write_values("binom_10_0.9.csv", N, values);
}

// test poisson with several different means
void test_poisson() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_poisson(1);
    }
    write_values("poisson_1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_poisson(5);
    }
    write_values("poisson_5.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_poisson(10);
    }
    write_values("poisson_10.csv", N, values);
}

// test geometric
void test_geom() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_geom(0.1);
    }
    write_values("geom_0.1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_geom(0.5);
    }
    write_values("geom_0.5.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_geom(0.9);
    }
    write_values("geom_0.9.csv", N, values);
}

// test gamma with several different shapes
void test_gamma() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_gamma(1, 1);
    }
    write_values("gamma_1_1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_gamma(2, 1);
    }
    write_values("gamma_2_1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_gamma(3, 1);
    }
    write_values("gamma_3_1.csv", N, values);
}

// test beta with several different shapes
void test_beta() {
    const int N = 100000;
    float values[N];
    for (int i = 0; i < N; i++) {
        values[i] = rand_beta(1, 1);
    }
    write_values("beta_1_1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_beta(2, 1);
    }
    write_values("beta_2_1.csv", N, values);
    for (int i = 0; i < N; i++) {
        values[i] = rand_beta(3, 1);
    }
    write_values("beta_3_1.csv", N, values);
}

int main(const int argc, const char** argv) {
    // run all the tests
    test_norm();
    test_skew_norm();
    test_binom();
    test_poisson();
    test_geom();
    test_gamma();
    test_beta();

    return 0;
}