From 249a1ff4344b294fa3b55d7af52df11917ff2f8f Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Thu, 2 Nov 2023 23:24:36 +0000 Subject: [PATCH] initial attempt on bc buggy because wrong base for log, but it's a start --- bc/estimate.bc | 34 +++++++++++++++++++++++ bc/extra/beta.bc | 51 ++++++++++++++++++++++++++++++++++ bc/{ => extra}/scratchpad.bc | 0 bc/makefile | 5 ++++ bc/notes.md | 2 +- bc/squiggle.bc | 54 +----------------------------------- 6 files changed, 92 insertions(+), 54 deletions(-) create mode 100644 bc/estimate.bc create mode 100644 bc/extra/beta.bc rename bc/{ => extra}/scratchpad.bc (100%) create mode 100644 bc/makefile diff --git a/bc/estimate.bc b/bc/estimate.bc new file mode 100644 index 00000000..afd15739 --- /dev/null +++ b/bc/estimate.bc @@ -0,0 +1,34 @@ +p_a = 0.8 +p_b = 0.5 +p_c = p_a * p_b + +weights[0] = 1 - p_c +weights[1] = p_c / 2 +weights[2] = p_c / 4 +weights[3] = p_c / 4 + +/* We'll have to define the mixture manually */ +define mixture(){ + p = sample_unit_uniform() + if(p <= weights[0]){ + return 0 + } + if(p <= (weights[0] + weights[1])){ + return 1 + } + if(p<= (weights[0] + weights[1] + weights[2])){ + return sample_to(1, 3) + } + return sample_to(2, 10) +} + +/* n_samples = 1000000 */ +n_samples = 10000 +sum=0 +for(i=0; i < n_samples; i++){ + /* samples[i] = mixture() */ + sum += mixture() +} +sum/n_samples + +halt diff --git a/bc/extra/beta.bc b/bc/extra/beta.bc new file mode 100644 index 00000000..62269a96 --- /dev/null +++ b/bc/extra/beta.bc @@ -0,0 +1,51 @@ +define sample_gamma(alpha){ + /* + A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001 + https://dl.acm.org/doi/pdf/10.1145/358407.358414 + see also the references/ folder + Note that the Wikipedia page for the gamma distribution includes a scaling parameter + k or beta + https://en.wikipedia.org/wiki/Gamma_distribution + such that gamma_k(alpha, k) = k * gamma(alpha) + or gamma_beta(alpha, beta) = gamma(alpha) / beta + So far I have not needed to use this, and thus the second parameter is by default 1. + */ + + if (alpha >= 1) { + d = alpha - (1/3); + c = 1.0 / sqrt(9.0 * d); + while (1) { + v=-1 + while(v<=0) { + x = sample_unit_normal(); + v = 1 + c * x; + } + v = v * v * v; + u = sample_unit_uniform(); + if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */ + /* + the 0.0331 doesn't inspire much confidence + however, this isn't the whole story + by knowing that Condition 1 implies condition 2 + we realize that this is just a way of making the algorithm faster + i.e., of not using the logarithms + */ + return d * v; + } + if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */ + return d * v; + } + } + } else { + return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha); + /* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */ + } +} + +define sample_beta(a, b) +{ + /* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */ + gamma_a = sample_gamma(a); + gamma_b = sample_gamma(b); + return gamma_a / (gamma_a + gamma_b); +} diff --git a/bc/scratchpad.bc b/bc/extra/scratchpad.bc similarity index 100% rename from bc/scratchpad.bc rename to bc/extra/scratchpad.bc diff --git a/bc/makefile b/bc/makefile new file mode 100644 index 00000000..1370b20b --- /dev/null +++ b/bc/makefile @@ -0,0 +1,5 @@ +compute: + ghbc -l squiggle.bc estimate.bc + +time: + time ghbc -l squiggle.bc estimate.bc diff --git a/bc/notes.md b/bc/notes.md index 727b8cc4..1a3c2578 100644 --- a/bc/notes.md +++ b/bc/notes.md @@ -6,7 +6,7 @@ https://pubs.opengroup.org/onlinepubs/9699919799.2018edition/utilities/bc.html ## gh-bc To build -./configure +./configure.sh -O3 make sudo cp bin/bc /usr/bin/ghbc diff --git a/bc/squiggle.bc b/bc/squiggle.bc index aab8f921..51f1c287 100644 --- a/bc/squiggle.bc +++ b/bc/squiggle.bc @@ -1,4 +1,4 @@ -scale = 32 +scale = 8 /* seed = 12345678910 */ pi = 4 * atan(1) normal90confidence=1.6448536269514727148638489079916 @@ -70,55 +70,3 @@ define sample_to(low, high){ loghigh = log(high, 2) return e(sample_normal_from_90_confidence_interval(loglow, loghigh)) } - -define sample_gamma(alpha){ - /* - A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001 - https://dl.acm.org/doi/pdf/10.1145/358407.358414 - see also the references/ folder - Note that the Wikipedia page for the gamma distribution includes a scaling parameter - k or beta - https://en.wikipedia.org/wiki/Gamma_distribution - such that gamma_k(alpha, k) = k * gamma(alpha) - or gamma_beta(alpha, beta) = gamma(alpha) / beta - So far I have not needed to use this, and thus the second parameter is by default 1. - */ - - if (alpha >= 1) { - d = alpha - (1/3); - c = 1.0 / sqrt(9.0 * d); - while (1) { - v=-1 - while(v<=0) { - x = sample_unit_normal(); - v = 1 + c * x; - } - v = v * v * v; - u = sample_unit_uniform(); - if (u < (1 - (0.0331 * (x * x * x * x)))) { /* Condition 1 */ - /* - the 0.0331 doesn't inspire much confidence - however, this isn't the whole story - by knowing that Condition 1 implies condition 2 - we realize that this is just a way of making the algorithm faster - i.e., of not using the logarithms - */ - return d * v; - } - if (log(u, 2) < ((0.5 * (x * x)) + (d * (1 - v + log(v, 2))))) { /* Condition 2 */ - return d * v; - } - } - } else { - return sample_gamma(1 + alpha) * p(sample_unit_uniform(), 1 / alpha); - /* see note in p. 371 of https://dl.acm.org/doi/pdf/10.1145/358407.358414 */ - } -} - -define sample_beta(a, b) -{ - /* See: https://en.wikipedia.org/wiki/Gamma_distribution#Related_distributions */ - gamma_a = sample_gamma(a); - gamma_b = sample_gamma(b); - return gamma_a / (gamma_a + gamma_b); -}