From 2023e0ba823960565349c12430ea65aeb067ba6c Mon Sep 17 00:00:00 2001 From: kellycochran Date: Sun, 26 Jan 2025 14:26:31 -0800 Subject: [PATCH] updated Random chugin --- Random/Random-test.ck | 244 +++++++++++++++++++++++++++++++++++++++--- Random/Random.cpp | 66 ++++++++++++ 2 files changed, 294 insertions(+), 16 deletions(-) diff --git a/Random/Random-test.ck b/Random/Random-test.ck index a6f74eab..4d1b3030 100644 --- a/Random/Random-test.ck +++ b/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(); \ No newline at end of file diff --git a/Random/Random.cpp b/Random/Random.cpp index 1e88626a..8d0a8a0f 100644 --- a/Random/Random.cpp +++ b/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; +}