diff --git a/squiggle.c b/squiggle.c index 8f77636..51076c4 100644 --- a/squiggle.c +++ b/squiggle.c @@ -20,7 +20,7 @@ uint64_t xorshift32(uint32_t* seed) // See // https://en.wikipedia.org/wiki/Xorshift // Also some drama: , - + // for floats uint64_t x = *seed; x ^= x << 13; x ^= x >> 17; @@ -30,11 +30,7 @@ uint64_t xorshift32(uint32_t* seed) uint64_t xorshift64(uint64_t* seed) { - // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" - // See - // https://en.wikipedia.org/wiki/Xorshift - // Also some drama: , - + // same as above, but for generating doubles uint64_t x = *seed; x ^= x << 13; x ^= x >> 7; @@ -383,6 +379,19 @@ struct box sampler_cdf_double(double cdf(double), uint64_t* seed) return result; } +/* Could also define other variations, e.g., +double sampler_danger(struct box cdf(double), uint64_t* seed) +{ + double p = sample_unit_uniform(seed); + struct box result = inverse_cdf_box(cdf, p); + if(result.empty){ + exit(1); + }else{ + return result.content; + } +} +*/ + // Get confidence intervals, given a sampler struct c_i { @@ -422,15 +431,47 @@ struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* se return result; } -/* Could also define other variations, e.g., -double sampler_danger(struct box cdf(double), uint64_t* seed) +// Do algebra over lognormals and normals +struct normal_parameters { + double mean; + double std; +}; + +struct lognormal_parameters { + double logmean; + double logstd; +}; + +struct normal_parameters algebra_sum_normals(struct normal_parameters a, struct normal_parameters b) { - double p = sample_unit_uniform(seed); - struct box result = inverse_cdf_box(cdf, p); - if(result.empty){ - exit(1); - }else{ - return result.content; - } + struct normal_parameters result = { + .mean = a.mean + b.mean, + .std = sqrt((a.std * a.std) + (b.std * b.std)), + }; + return result; +} +struct normal_parameters algebra_shift_normal(struct normal_parameters a, double shift) +{ + struct normal_parameters result = { + .mean = a.mean + shift, + .std = a.std, + }; + return result; +} + +struct lognormal_parameters algebra_product_lognormals(struct lognormal_parameters a, struct lognormal_parameters b) +{ + struct lognormal_parameters result = { + .logmean = a.logmean + b.logmean, + .logstd = sqrt((a.logstd * a.logstd) + (b.logstd * b.logstd)), + }; + return result; +} +struct lognormal_parameters algebra_scale_lognormal(struct lognormal_parameters a, double k) +{ + struct lognormal_parameters result = { + .logmean = a.logmean + k, + .logstd = a.logstd, + }; + return result; } -*/ diff --git a/squiggle.h b/squiggle.h index 5352069..8d8bd82 100644 --- a/squiggle.h +++ b/squiggle.h @@ -52,8 +52,8 @@ struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed); // Get 90% confidence interval struct c_i { - float low; - float high; + float low; + float high; }; struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed);