time-to-botec/bc/squiggle.bc

68 lines
2.4 KiB
Plaintext

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))
}