squiggle.c

Self-contained Monte Carlo estimation in C99
Log | Files | Refs | README

commit 6b2349132b284808e4c6456114d7c2ff7a4876be
parent b80b05ca30e1ae5cecfd27a24a26772425341581
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Sun, 23 Jul 2023 15:43:35 +0200

add tests lognormal, and have them use special tolerances.

Diffstat:
Mtest/makefile | 2+-
Mtest/test | 0
Mtest/test.c | 114++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------
3 files changed, 85 insertions(+), 31 deletions(-)

diff --git a/test/makefile b/test/makefile @@ -37,7 +37,7 @@ run: $(SRC) $(OUTPUT) ./$(OUTPUT) verify: $(SRC) $(OUTPUT) - ./$(OUTPUT) | grep "NOT passed" || true + ./$(OUTPUT) | grep "NOT passed" -A 1 --group-separator='' || true time-linux: @echo "Requires /bin/time, found on GNU/Linux systems" && echo diff --git a/test/test b/test/test Binary files differ. diff --git a/test/test.c b/test/test.c @@ -4,9 +4,8 @@ #include <stdlib.h> #include <stdio.h> -#define N 1000 * 1000 -#define PERCENTAGE_TOLERANCE_UNIFORM 1.0/1000.0 -#define PERCENTAGE_TOLERANCE_NORMAL 5.0/1000.0 +#define PERCENTAGE_TOLERANCE 1.0/1000.0 +#define PERCENTAGE_TOLERANCE_LOGNORMAL 5.0/1000.0 #define MAX_NAME_LENGTH 500 // Structs @@ -28,16 +27,18 @@ void test_array_expectations(struct array_expectations e){ double delta_std = std - e.expected_std; - if(fabs(delta_mean) > e.tolerance){ + if(fabs(delta_mean)/fabs(mean) > e.tolerance){ printf("[-] Mean test for %s NOT passed.\n", e.name); - printf("Mean of %s: %f, vs expected mean: %f, delta: %f\n", e.name, mean, e.expected_mean, delta_mean); + printf("Mean of %s: %f, vs expected mean: %f\n", e.name, mean, e.expected_mean); + printf("delta: %f, relative delta: %f\n", delta_mean, delta_mean/fabs(mean)); }else { printf("[x] Mean test for %s PASSED\n", e.name); } - if(fabs(delta_std) > e.tolerance){ + if(fabs(delta_std)/fabs(std) > e.tolerance){ printf("[-] Std test for %s NOT passed.\n", e.name); - printf("Std of %s: %f, vs expected std: %f, delta: %f\n", e.name, std, e.expected_std, delta_std); + printf("Std of %s: %f, vs expected std: %f\n", e.name, std, e.expected_std); + printf("delta: %f, relative delta: %f\n", delta_std, delta_std/fabs(std)); }else { printf("[x] Std test for %s PASSED\n", e.name); } @@ -48,29 +49,32 @@ void test_array_expectations(struct array_expectations e){ // Test unit uniform void test_unit_uniform(uint64_t* seed){ - double* unit_uniform_array = malloc(sizeof(double) * N); + int n = 1000 * 1000; + double* unit_uniform_array = malloc(sizeof(double) * n); - for(int i=0; i<N; i++){ + for(int i=0; i<n; i++){ unit_uniform_array[i] = sample_unit_uniform(seed); } struct array_expectations expectations = { .array = unit_uniform_array, - .n = N, + .n = n, .name = "unit uniform", .expected_mean = 0.5, .expected_std = sqrt(1.0/12.0), - .tolerance = 1 * PERCENTAGE_TOLERANCE_UNIFORM, + .tolerance = 1 * PERCENTAGE_TOLERANCE, }; test_array_expectations(expectations); + free(unit_uniform_array); } // Test uniforms void test_uniform(double start, double end, uint64_t* seed){ - double* uniform_array = malloc(sizeof(double) * N); + int n = 1000 * 1000; + double* uniform_array = malloc(sizeof(double) * n); - for(int i=0; i<N; i++){ + for(int i=0; i<n; i++){ uniform_array[i] = sample_uniform(start, end, seed); } @@ -78,42 +82,46 @@ void test_uniform(double start, double end, uint64_t* seed){ snprintf(name, MAX_NAME_LENGTH, "[%f, %f] uniform", start, end); struct array_expectations expectations = { .array = uniform_array, - .n = N, + .n = n, .name = name, .expected_mean = (start + end)/2, .expected_std = sqrt(1.0/12.0) * fabs(end-start), - .tolerance = fabs(end -start) * PERCENTAGE_TOLERANCE_UNIFORM, + .tolerance = fabs(end -start) * PERCENTAGE_TOLERANCE, }; test_array_expectations(expectations); free(name); + free(uniform_array); } // Test unit normal void test_unit_normal(uint64_t* seed){ - double* unit_normal_array = malloc(sizeof(double) * N); + int n = 1000 * 1000; + double* unit_normal_array = malloc(sizeof(double) * n); - for(int i=0; i<N; i++){ + for(int i=0; i<n; i++){ unit_normal_array[i] = sample_unit_normal(seed); } struct array_expectations expectations = { .array = unit_normal_array, - .n = N, + .n = n, .name = "unit normal", .expected_mean = 0, .expected_std = 1, - .tolerance = 1 * PERCENTAGE_TOLERANCE_NORMAL, + .tolerance = 1 * PERCENTAGE_TOLERANCE, }; test_array_expectations(expectations); + free(unit_normal_array); } // Test normal void test_normal(double mean, double std, uint64_t* seed){ - double* normal_array = malloc(sizeof(double) * N); + int n = 10 * 1000 * 1000; + double* normal_array = malloc(sizeof(double) * n); - for(int i=0; i<N; i++){ + for(int i=0; i<n; i++){ normal_array[i] = sample_normal(mean, std, seed); } @@ -121,23 +129,50 @@ void test_normal(double mean, double std, uint64_t* seed){ snprintf(name, MAX_NAME_LENGTH, "normal(%f, %f)", mean, std); struct array_expectations expectations = { .array = normal_array, - .n = N, + .n = n, .name = name, .expected_mean = mean, .expected_std = std, - .tolerance = std * PERCENTAGE_TOLERANCE_NORMAL, + .tolerance = std * PERCENTAGE_TOLERANCE, }; test_array_expectations(expectations); free(name); + free(normal_array); +} + +// Test lognormal +void test_lognormal(double logmean, double logstd, uint64_t* seed){ + int n = 10 * 1000 * 1000; + double* lognormal_array = malloc(sizeof(double) * n); + + for(int i=0; i<n; i++){ + lognormal_array[i] = sample_lognormal(logmean, logstd, seed); + } + + char* name = malloc(MAX_NAME_LENGTH * sizeof(char)); + snprintf(name, MAX_NAME_LENGTH, "lognormal(%f, %f)", logmean, logstd); + struct array_expectations expectations = { + .array = lognormal_array, + .n = n, + .name = name, + .expected_mean = exp(logmean + pow(logstd, 2)/2), + .expected_std = sqrt((exp(pow(logstd, 2)) - 1) * exp(2*logmean + pow(logstd, 2))), + .tolerance = exp(logstd) * PERCENTAGE_TOLERANCE_LOGNORMAL, + }; + + test_array_expectations(expectations); + free(name); + free(lognormal_array); } // Test beta void test_beta(double a, double b, uint64_t* seed){ - double* beta_array = malloc(sizeof(double) * N); + int n = 10 * 1000 * 1000; + double* beta_array = malloc(sizeof(double) * n); - for(int i=0; i<N; i++){ + for(int i=0; i<n; i++){ beta_array[i] = sample_beta(a, b, seed); } @@ -145,11 +180,11 @@ void test_beta(double a, double b, uint64_t* seed){ snprintf(name, MAX_NAME_LENGTH, "beta(%f, %f)", a, b); struct array_expectations expectations = { .array = beta_array, - .n = N, + .n = n, .name = name, .expected_mean = a/(a+b), .expected_std = sqrt((a*b)/(pow(a+b, 2) * (a + b + 1))), - .tolerance = PERCENTAGE_TOLERANCE_UNIFORM, + .tolerance = PERCENTAGE_TOLERANCE, }; test_array_expectations(expectations); @@ -160,7 +195,7 @@ int main(){ // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with a seed of 0 - + /* printf("Testing unit uniform\n"); test_unit_uniform(seed); @@ -187,8 +222,8 @@ int main(){ printf("Testing small normals\n"); for(int i=0; i<100; i++){ - double mean = sample_normal(-10, 10, seed); - double std = sample_normal(0, 10, seed); + double mean = sample_uniform(-10, 10, seed); + double std = sample_uniform(0, 10, seed); if ( std > 0){ test_normal(mean, std, seed); } @@ -202,7 +237,25 @@ int main(){ test_normal(mean, std, seed); } } + */ + printf("Testing very small lognormals\n"); + for(int i=0; i<10; i++){ + double mean = sample_uniform(-1, 1, seed); + double std = sample_uniform(0, 1, seed); + if ( std > 0){ + test_lognormal(mean, std, seed); + } + } + printf("Testing small lognormals\n"); + for(int i=0; i<10; i++){ + double mean = sample_uniform(-1, 5, seed); + double std = sample_uniform(0, 5, seed); + if ( std > 0){ + test_lognormal(mean, std, seed); + } + } + /* printf("Testing beta distribution\n"); for(int i=0; i<100; i++){ double a = sample_uniform(0, 1000, seed); @@ -222,5 +275,6 @@ int main(){ } free(seed); + */ }