squiggle.bc (2480B)
1 scale = 16 2 /* seed = 12345678910 */ 3 pi = 4 * atan(1) 4 normal90confidence=1.64485362695 5 /* 1.6448536269514727148638489079916 */ 6 7 /* Distribution & sampling functions */ 8 /* Unit distributions */ 9 define sample_unit_uniform(){ 10 return rand()/maxrand() 11 } 12 define sample_unit_normal(){ 13 u1=sample_unit_uniform() 14 u2=sample_unit_uniform() 15 z = sqrt(-2 * l(u1)) * sin(2 * pi * u2) 16 return z 17 } 18 19 /* Composite distributions */ 20 define sample_uniform(min, max){ 21 return (min + sample_unit_uniform()*(max-min)) 22 } 23 define sample_normal(mean, sigma){ 24 return (mean + sigma * sample_unit_normal()) 25 } 26 define sample_lognormal(logmean, logstd){ 27 return e(sample_normal(logmean, logstd)) 28 } 29 define sample_normal_from_90_confidence_interval(low, high){ 30 /* 31 Explanation of key idea: 32 1. We know that the 90% confidence interval of the unit normal is 33 [-1.6448536269514722, 1.6448536269514722] 34 see e.g.: https://stackoverflow.com/questions/20626994/how-to-calculate-the-inverse-of-the-normal-cumulative-distribution-function-in-p 35 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 36 2. So if we take a unit normal and multiply it by 37 L / 1.6448536269514722, its new 90% confidence interval will be 38 [-L, L], i.e., length 2 * L 39 3. Instead, if we want to get a confidence interval of length L, 40 we should multiply the unit normal by 41 L / (2 * 1.6448536269514722) 42 Meaning that its standard deviation should be multiplied by that amount 43 see: https://en.wikipedia.org/wiki/Normal_distribution?lang=en#Operations_on_a_single_normal_variable 44 4. So we have learnt that Normal(0, L / (2 * 1.6448536269514722)) 45 has a 90% confidence interval of length L 46 5. If we want a 90% confidence interval from high to low, 47 we can set mean = (high + low)/2; the midpoint, and L = high-low, 48 Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722)) 49 */ 50 mean = (high + low) / 2 51 std = (high - low) / (2 * normal90confidence) 52 return sample_normal(mean, std) 53 } 54 define sample_to(low, high){ 55 56 /* 57 Given a (positive) 90% confidence interval, 58 returns a sample from a lognorma with a matching 90% c.i. 59 Key idea: If we want a lognormal with 90% confidence interval [a, b] 60 we need but get a normal with 90% confidence interval [log(a), log(b)]. 61 Then see code for sample_normal_from_90_confidence_interval 62 */ 63 loglow = l(low) 64 loghigh = l(high) 65 return e(sample_normal_from_90_confidence_interval(loglow, loghigh)) 66 } 67