squiggle.c

commit 555ac7371b8961a282bd0eefb131458821806b51
parent e0df88633c74d03e7ff8f3018e9366bd7a02f663
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Tue,  1 Aug 2023 14:08:19 +0200

add shift for normals and scaling for lognormals

Diffstat:
M
squiggle.c
 | 
73
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------
M
squiggle.h
 | 
4
++--

2 files changed, 59 insertions(+), 18 deletions(-)
diff --git a/squiggle.c b/squiggle.c
@@ -20,7 +20,7 @@ uint64_t xorshift32(uint32_t* seed)
     // See <https://stackoverflow.com/questions/53886131/how-does-xorshift64-works>
     // https://en.wikipedia.org/wiki/Xorshift
     // Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/>
-
+    // 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://stackoverflow.com/questions/53886131/how-does-xorshift64-works>
-    // https://en.wikipedia.org/wiki/Xorshift
-    // Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/>
-
+    // 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
@@ -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);