time-to-botec/bc/squiggle.bc

scale = 16
/* seed = 12345678910 */
pi = 4 * atan(1)
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()
  z = sqrt(-2 * l(u1)) * sin(2 * pi * u2)
  return z
}

/* 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:
  1. We know that the 90% confidence interval of the unit normal is
  [-1.6448536269514722, 1.6448536269514722]
  see e.g.: https://stackoverflow.com/questions/20626994/how-to-calculate-the-inverse-of-the-normal-cumulative-distribution-function-in-p
  or https://www.wolframalpha.com/input?i=N%5BInverseCDF%28normal%280%2C1%29%2C+0.05%29%2C%7B%E2%88%9E%2C100%7D%5D
  2. So if we take a unit normal and multiply it by
  L / 1.6448536269514722, its new 90% confidence interval will be
  [-L, L], i.e., length 2 * L
  3. Instead, if we want to get a confidence interval of length L,
  we should multiply the unit normal by
  L / (2 * 1.6448536269514722)
  Meaning that its standard deviation should be multiplied by that amount
  see: https://en.wikipedia.org/wiki/Normal_distribution?lang=en#Operations_on_a_single_normal_variable
  4. So we have learnt that Normal(0, L / (2 * 1.6448536269514722))
  has a 90% confidence interval of length L
  5. If we want a 90% confidence interval from high to low,
  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
  std = (high - low) / (2 * normal90confidence)
  return sample_normal(mean, std)
}
define sample_to(low, high){

  /*
  Given a (positive) 90% confidence interval,
  returns a sample from a lognorma with a matching 90% c.i.
  Key idea: If we want a lognormal with 90% confidence interval [a, b]
  we need but get a normal with 90% confidence interval [log(a), log(b)].
  Then see code for sample_normal_from_90_confidence_interval
  */
  loglow = l(low)
  loghigh = l(high)
  return e(sample_normal_from_90_confidence_interval(loglow, loghigh))
}