Merge pull request #100 from ccrma/random

updated Random chugin

Random/Random-test.ck

@@ -1,22 +1,234 @@
+@import "Random.chug"
 
-0 => float sum;
-float nums[0];
-10000 => int count;
 
-repeat(count)
-{
-    Random.gaussian(0, 10) => float r;
-    chout <= r <= IO.nl();
-    r +=> sum;
-    nums << r;
+// will sample RVs, calc these stats, and compare to expected
+
+fun float calc_mean(int nums[]) {
+    0 => float sum;
+    for (int num : nums) sum + num => sum;
+    return sum / nums.size();
+}
+
+fun float calc_fmean(float nums[]) {
+    0 => float sum;
+    for (float num : nums) sum + num => sum;
+    return sum / nums.size();
+}
+
+fun float calc_fstd_dev(float nums[]) {
+    calc_fmean(nums) => float mean;
+
+    0 => float sum_sq_diff;
+    for (float num : nums) Math.pow(num - mean, 2) +=> sum_sq_diff;
+    return Math.sqrt(sum_sq_diff/nums.size());
+}
+
+fun float calc_std_dev(int nums[]) {
+    calc_mean(nums) => float mean;
+
+    0 => float sum_sq_diff;
+    for (int num : nums) Math.pow(num - mean, 2) +=> sum_sq_diff;
+    return Math.sqrt(sum_sq_diff/nums.size());
+}
+
+
+
+// sampling functions
+
+fun float[] sample_gaussian(float mean, float std_dev, int num_samples) {
+    float samples[num_samples];
+    for (int i; i < num_samples; i++) Random.gaussian(mean, std_dev) => samples[i];
+    return samples;
+}
+
+fun int[] sample_binom(int n, float p, int num_samples) {
+    int samples[num_samples];
+    for (int i; i < num_samples; i++) Random.binomial(n, p) => samples[i];
+    return samples;
+}
+
+fun float[] sample_exp(float scale, int num_samples) {
+    float samples[num_samples];
+    for (int i; i < num_samples; i++) Random.exponential(scale) => samples[i];
+    return samples;
+}
+
+fun int[] sample_geom(float p, int num_samples) {
+    int samples[num_samples];
+    for (int i; i < num_samples; i++) Random.geometric(p) => samples[i];
+    return samples;
+}
+
+fun int[] sample_poisson(float p, int num_samples) {
+    int samples[num_samples];
+    for (int i; i < num_samples; i++) Random.poisson(p) => samples[i];
+    return samples;
 }
 
-sum/count => float mean;
-0 => float sum_sq_diff;
-// re-iterate through sequence and calculate std deviation
-for(0 => int i; i < nums.size(); i++)
-    (nums[i]-mean)*(nums[i]-mean) +=> sum_sq_diff;
 
-chout <= "mean value: " <= mean <= IO.nl();
-chout <= "std deviation: " <= Math.sqrt(sum_sq_diff/count) <= IO.nl();
 
+// unit testing
+
+fun void compare_stats_gaussian() {
+    <<< "Gaussian Test", "" >>>;
+
+    1. => float mean;
+    5. => float std_dev;
+    100000 => int num_samples;
+
+    sample_gaussian(mean, std_dev, num_samples) @=> float samples[];
+
+    calc_fmean(samples) => float sample_mean;
+    calc_fstd_dev(samples) => float sample_std_dev;
+
+    chout <= "True mean: " <= mean <= IO.nl();
+    chout <= "Obs. mean: " <= sample_mean <= IO.nl();
+
+    chout <= "True std. dev.: " <= std_dev <= IO.nl();
+    chout <= "Obs. std. dev.: " <= sample_std_dev <= IO.nl();
+
+    Math.fabs((mean - sample_mean) / mean) => float mean_off;
+    Math.fabs((std_dev - sample_std_dev) / std_dev) => float std_dev_off;
+
+    if (mean_off < 0.05 && std_dev_off < 0.05) {
+        <<< "=== Test: PASS ===", "\n" >>>;
+    } else {
+        <<< "=== Test: FAIL ===", "\n" >>>;
+    }
+}
+
+
+fun void compare_stats_binom() {
+    <<< "Binomial Test", "" >>>;
+
+    20 => int n;
+    0.4 => float p;
+    100000 => int num_samples;
+
+    sample_binom(n, p, num_samples) @=> int samples[];
+
+    calc_mean(samples) => float sample_mean;
+    calc_std_dev(samples) => float sample_std_dev;
+
+    n * p => float mean;
+    Math.pow(n * p * (1-p), 0.5) => float std_dev;
+
+    chout <= "True mean: " <= mean <= IO.nl();
+    chout <= "Obs. mean: " <= sample_mean <= IO.nl();
+
+    chout <= "True std. dev.: " <= std_dev <= IO.nl();
+    chout <= "Obs. std. dev.: " <= sample_std_dev <= IO.nl();
+
+    Math.fabs((mean - sample_mean) / mean) => float mean_off;
+    Math.fabs((std_dev - sample_std_dev) / std_dev) => float std_dev_off;
+
+    if (mean_off < 0.05 && std_dev_off < 0.05) {
+        <<< "=== Test: PASS ===", "\n" >>>;
+    } else {
+        <<< "=== Test: FAIL ===", "\n" >>>;
+    }
+}
+
+
+fun void compare_stats_exp() {
+    <<< "Exponential Test", "" >>>;
+
+    3. => float scale; // 1 / lambda, not lambda
+    100000 => int num_samples;
+
+    sample_exp(scale, num_samples) @=> float samples[];
+
+    calc_fmean(samples) => float sample_mean;
+    calc_fstd_dev(samples) => float sample_std_dev;
+
+    scale => float mean;
+    scale => float std_dev;
+
+    chout <= "True mean: " <= mean <= IO.nl();
+    chout <= "Obs. mean: " <= sample_mean <= IO.nl();
+
+    chout <= "True std. dev.: " <= std_dev <= IO.nl();
+    chout <= "Obs. std. dev.: " <= sample_std_dev <= IO.nl();
+
+    Math.fabs((mean - sample_mean) / mean) => float mean_off;
+    Math.fabs((std_dev - sample_std_dev) / std_dev) => float std_dev_off;
+
+    if (mean_off < 0.05 && std_dev_off < 0.05) {
+        <<< "=== Test: PASS ===", "\n" >>>;
+    } else {
+        <<< "=== Test: FAIL ===", "\n" >>>;
+    }
+}
+
+
+
+fun void compare_stats_geom() {
+    <<< "Geometric Test", "" >>>;
+
+    0.2 => float p;
+    100000 => int num_samples;
+
+    sample_geom(p, num_samples) @=> int samples[];
+
+    calc_mean(samples) => float sample_mean;
+    calc_std_dev(samples) => float sample_std_dev;
+
+    1. / p => float mean;
+    Math.pow(1 - p, 0.5) / p => float std_dev;
+
+    chout <= "True mean: " <= mean <= IO.nl();
+    chout <= "Obs. mean: " <= sample_mean <= IO.nl();
+
+    chout <= "True std. dev.: " <= std_dev <= IO.nl();
+    chout <= "Obs. std. dev.: " <= sample_std_dev <= IO.nl();
+
+    Math.fabs((mean - sample_mean) / mean) => float mean_off;
+    Math.fabs((std_dev - sample_std_dev) / std_dev) => float std_dev_off;
+
+    if (mean_off < 0.05 && std_dev_off < 0.05) {
+        <<< "=== Test: PASS ===", "\n" >>>;
+    } else {
+        <<< "=== Test: FAIL ===", "\n" >>>;
+    }
+}
+
+
+fun void compare_stats_poisson() {
+    <<< "Poisson Test", "" >>>;
+
+    3 => float lambda;
+    100000 => int num_samples;
+
+    sample_poisson(lambda, num_samples) @=> int samples[];
+
+    calc_mean(samples) => float sample_mean;
+    calc_std_dev(samples) => float sample_std_dev;
+
+    lambda => float mean;
+    Math.pow(lambda, 0.5) => float std_dev;
+
+    chout <= "True mean: " <= mean <= IO.nl();
+    chout <= "Obs. mean: " <= sample_mean <= IO.nl();
+
+    chout <= "True std. dev.: " <= std_dev <= IO.nl();
+    chout <= "Obs. std. dev.: " <= sample_std_dev <= IO.nl();
+
+    Math.fabs((mean - sample_mean) / mean) => float mean_off;
+    Math.fabs((std_dev - sample_std_dev) / std_dev) => float std_dev_off;
+
+    if (mean_off < 0.05 && std_dev_off < 0.05) {
+        <<< "=== Test: PASS ===", "\n" >>>;
+    } else {
+        <<< "=== Test: FAIL ===", "\n" >>>;
+    }
+}
+
+
+
+// If these all print "PASS", you're good-ish
+
+compare_stats_gaussian();
+compare_stats_binom();
+compare_stats_exp();
+compare_stats_geom();
+compare_stats_poisson();

Random/Random.cpp

@@ -28,6 +28,10 @@ static void srandom( unsigned s ) { srand( s ); }
 // example of getter/setter
 CK_DLL_SFUN(Random_seed);
 CK_DLL_SFUN(Random_gaussian);
+CK_DLL_SFUN(Random_binomial);
+CK_DLL_SFUN(Random_exponential);
+CK_DLL_SFUN(Random_geometric);
+CK_DLL_SFUN(Random_poisson);
 
 
 CK_DLL_QUERY( Random )
@@ -44,6 +48,19 @@ CK_DLL_QUERY( Random )
     QUERY->add_arg(QUERY, "float", "mean");
     QUERY->add_arg(QUERY, "float", "stdv");
 
+    QUERY->add_sfun(QUERY, Random_binomial, "int", "binomial");
+    QUERY->add_arg(QUERY, "int", "n");
+    QUERY->add_arg(QUERY, "float", "p");
+
+    QUERY->add_sfun(QUERY, Random_exponential, "float", "exponential");
+    QUERY->add_arg(QUERY, "float", "scale");
+
+    QUERY->add_sfun(QUERY, Random_geometric, "int", "geometric");
+    QUERY->add_arg(QUERY, "float", "p");
+
+    QUERY->add_sfun(QUERY, Random_poisson, "int", "poisson");
+    QUERY->add_arg(QUERY, "float", "lambda");
+
     // end the class definition
     // IMPORTANT: this MUST be called!
     QUERY->end_class(QUERY);
@@ -75,3 +92,52 @@ CK_DLL_SFUN(Random_gaussian)
     RETURN->v_float = mean+Z0*stdv;
 }
 
+CK_DLL_SFUN(Random_binomial)
+{
+    t_CKFLOAT n = GET_NEXT_INT(ARGS);
+    t_CKFLOAT p = GET_NEXT_FLOAT(ARGS);
+
+    t_CKINT count = 0;
+    for (t_CKINT i = 0; i < n; ++i) {
+        t_CKFLOAT U = random()/((float) CK_RANDOM_MAX);
+        if (U < p) count++;
+    }
+    RETURN->v_int = count;
+}
+
+CK_DLL_SFUN(Random_exponential)
+{
+    t_CKFLOAT scale = GET_NEXT_FLOAT(ARGS);
+
+    t_CKFLOAT U = random()/((float) CK_RANDOM_MAX);
+    RETURN->v_float = -1 * scale * log(U);
+}
+
+CK_DLL_SFUN(Random_geometric)
+{
+    t_CKFLOAT p = GET_NEXT_FLOAT(ARGS);
+
+    t_CKINT count = 1;
+    while (1) {
+        t_CKFLOAT U = random()/((float) CK_RANDOM_MAX);
+        if (U < p) break;
+        count++;
+    }
+    RETURN->v_int = count;
+}
+
+CK_DLL_SFUN(Random_poisson)
+{
+    t_CKFLOAT lambda = GET_NEXT_FLOAT(ARGS);
+
+    t_CKINT x = 0;
+    t_CKFLOAT p = exp(- lambda);
+    t_CKFLOAT s = p;
+    t_CKFLOAT U = random()/((float) CK_RANDOM_MAX);
+    while (U > s) {
+        x++;
+        p = p * lambda / x;
+        s = s + p;
+    }
+    RETURN->v_int = x;
+}