From 90b880488472c95fd392db609d799f9d035909cb Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.sempere@protonmail.com>
Date: Thu, 2 Nov 2023 23:34:48 +0000
Subject: [PATCH] arrive at working version of squiggle.bc

---
 bc/estimate.bc |  2 +-
 bc/makefile    |  2 ++
 bc/squiggle.bc | 17 ++++++-----------
 3 files changed, 9 insertions(+), 12 deletions(-)

diff --git a/bc/estimate.bc b/bc/estimate.bc
index afd15739..574fda5a 100644
--- a/bc/estimate.bc
+++ b/bc/estimate.bc
@@ -23,7 +23,7 @@ define mixture(){
 }
 
 /* n_samples = 1000000 */
-n_samples = 10000
+n_samples = 1000000
 sum=0
 for(i=0; i < n_samples; i++){
   /* samples[i] = mixture() */
diff --git a/bc/makefile b/bc/makefile
index 1370b20b..4b67150b 100644
--- a/bc/makefile
+++ b/bc/makefile
@@ -1,3 +1,5 @@
+SHELL := /bin/bash
+
 compute:
 	ghbc -l squiggle.bc estimate.bc
 
diff --git a/bc/squiggle.bc b/bc/squiggle.bc
index 70f00a9a..b570f2f3 100644
--- a/bc/squiggle.bc
+++ b/bc/squiggle.bc
@@ -1,15 +1,14 @@
-scale = 8
+scale = 16
 /* seed = 12345678910 */
 pi = 4 * atan(1)
-normal90confidence=1.6448536269514727148638489079916
+normal90confidence=1.64485362695
+/* 1.6448536269514727148638489079916 */
 
 /* Distribution & sampling functions */
 /* Unit distributions */
-
 define sample_unit_uniform(){
   return rand()/maxrand()
 }
-
 define sample_unit_normal(){
   u1=sample_unit_uniform()
   u2=sample_unit_uniform()
@@ -18,19 +17,15 @@ define sample_unit_normal(){
 }
 
 /* Composite distributions */
-
 define sample_uniform(min, max){
   return (min + sample_unit_uniform()*(max-min))
 }
-
 define sample_normal(mean, sigma){
   return (mean + sigma * sample_unit_normal())
 }
-
 define sample_lognormal(logmean, logstd){
   return e(sample_normal(logmean, logstd))
 }
-
 define sample_normal_from_90_confidence_interval(low, high){
   /*
   Explanation of key idea:
@@ -52,11 +47,10 @@ define sample_normal_from_90_confidence_interval(low, high){
   we can set mean = (high + low)/2; the midpoint, and L = high-low,
   Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722))
   */
-  mean = (high + low) / 2.0
-  std = (high - low) / (2.0 * normal90confidence)
+  mean = (high + low) / 2
+  std = (high - low) / (2 * normal90confidence)
   return sample_normal(mean, std)
 }
-
 define sample_to(low, high){
 
   /*
@@ -70,3 +64,4 @@ define sample_to(low, high){
   loghigh = l(high)
   return e(sample_normal_from_90_confidence_interval(loglow, loghigh))
 }
+