samples.nim (2003B)
1 import std/math 2 import std/sugar 3 import std/random 4 import std/sequtils 5 6 randomize() 7 8 ## Distribution functions 9 10 ## Normal 11 ## <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform#Basic_form> 12 proc ur_normal(): float = 13 let u1 = rand(1.0) 14 let u2 = rand(1.0) 15 let z = sqrt(-2.0 * ln(u1)) * sin(2 * PI * u2) 16 return z 17 18 proc normal(mean: float, sigma: float): float = 19 return (mean + sigma * ur_normal()) 20 21 proc lognormal(logmean: float, logsigma: float): float = 22 let answer = pow(E, normal(logmean, logsigma)) 23 return answer 24 25 proc to(low: float, high: float): float = 26 let normal95confidencePoint = 1.6448536269514722 27 let loglow = ln(low) 28 let loghigh = ln(high) 29 let logmean = (loglow + loghigh)/2 30 let logsigma = (loghigh - loglow) / (2.0 * normal95confidencePoint) 31 return lognormal(logmean, logsigma) 32 33 ## echo ur_normal() 34 ## echo normal(10, 20) 35 ## echo lognormal(2, 4) 36 ## echo to(10, 90) 37 38 ## Manipulate samples 39 40 proc mixture(fs: seq[proc (): float{.nimcall.}], ps: seq[float], n: int): seq[float] = 41 42 assert fs.len == ps.len 43 44 var ws: seq[float] 45 var sum = 0.0 46 for i, p in ps: 47 sum = sum + p 48 ws.add(sum) 49 ws = ws.map(w => w/sum) 50 51 var samples: seq[float] 52 let rs = toSeq(1..n).map(_=>rand(1.0)) 53 for i in 0..(n-1): 54 let r = rs[i] 55 var j = ws.len - 1 56 for k, w in ws: 57 if r < w: 58 j = k 59 break 60 ## only occasion when ^ doesn't assign j 61 ## is when r is exactly 1 62 ## which would correspond to choosing the last item in ws 63 ## which is why j is initialized to ws.len - 1 64 let f = fs[j] 65 samples.add(f()) 66 return samples 67 68 69 ## Actual model 70 71 let n = 1000000 72 73 let p_a = 0.8 74 let p_b = 0.5 75 let p_c = p_a * p_b 76 77 let weights = @[ 1.0 - p_c, p_c/2.0, p_c/4.0, p_c/4.0 ] 78 79 let fs = @[ proc (): float = 0.0, proc (): float = 1.0, proc (): float = to(1.0, 3.0), proc (): float = to(2.0, 10.0)] 80 let result = mixture(fs, weights, n) 81 let mean_result = foldl(result, a + b, 0.0) / float(result.len) 82 83 # echo result 84 echo mean_result