personal/time-to-botec
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

arrive at working version of squiggle.bc

Browse Source
This commit is contained in:
NunoSempere 2023-11-02 23:34:48 +00:00
parent 8675d98784
commit 90b8804884
3 changed files with 9 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
bc/estimate.bc
View File

@ -23,7 +23,7 @@ define mixture(){
} }
/* n_samples = 1000000 */ /* n_samples = 1000000 */
n_samples = 10000 n_samples = 1000000
sum=0 sum=0
for(i=0; i < n_samples; i++){ for(i=0; i < n_samples; i++){
/* samples[i] = mixture() */ /* samples[i] = mixture() */

2
bc/makefile
View File

@ -1,3 +1,5 @@
SHELL := /bin/bash
compute: compute:
ghbc -l squiggle.bc estimate.bc ghbc -l squiggle.bc estimate.bc

17
bc/squiggle.bc
View File

@ -1,15 +1,14 @@
scale = 8 scale = 16
/* seed = 12345678910 */ /* seed = 12345678910 */
pi = 4 * atan(1) pi = 4 * atan(1)
normal90confidence=1.6448536269514727148638489079916 normal90confidence=1.64485362695
/* 1.6448536269514727148638489079916 */
/* Distribution & sampling functions */ /* Distribution & sampling functions */
/* Unit distributions */ /* Unit distributions */
define sample_unit_uniform(){ define sample_unit_uniform(){
return rand()/maxrand() return rand()/maxrand()
} }
define sample_unit_normal(){ define sample_unit_normal(){
u1=sample_unit_uniform() u1=sample_unit_uniform()
u2=sample_unit_uniform() u2=sample_unit_uniform()
@ -18,19 +17,15 @@ define sample_unit_normal(){
} }
/* Composite distributions */ /* Composite distributions */
define sample_uniform(min, max){ define sample_uniform(min, max){
return (min + sample_unit_uniform()*(max-min)) return (min + sample_unit_uniform()*(max-min))
} }
define sample_normal(mean, sigma){ define sample_normal(mean, sigma){
return (mean + sigma * sample_unit_normal()) return (mean + sigma * sample_unit_normal())
} }
define sample_lognormal(logmean, logstd){ define sample_lognormal(logmean, logstd){
return e(sample_normal(logmean, logstd)) return e(sample_normal(logmean, logstd))
} }
define sample_normal_from_90_confidence_interval(low, high){ define sample_normal_from_90_confidence_interval(low, high){
/* /*
Explanation of key idea: 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, we can set mean = (high + low)/2; the midpoint, and L = high-low,
Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722)) Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722))
*/ */
mean = (high + low) / 2.0 mean = (high + low) / 2
std = (high - low) / (2.0 * normal90confidence) std = (high - low) / (2 * normal90confidence)
return sample_normal(mean, std) return sample_normal(mean, std)
} }
define sample_to(low, high){ define sample_to(low, high){
/* /*
@ -70,3 +64,4 @@ define sample_to(low, high){
loghigh = l(high) loghigh = l(high)
return e(sample_normal_from_90_confidence_interval(loglow, loghigh)) return e(sample_normal_from_90_confidence_interval(loglow, loghigh))
} }