samples.nim (3421B)
1 import std/math 2 import std/sugar 3 import std/random 4 import std/sequtils 5 6 randomize() 7 8 ## Basic math functions 9 proc factorial(n: int): int = 10 if n == 0 or n < 0: 11 return 1 12 else: 13 return n * factorial(n - 1) 14 15 proc sine(x: float): float = 16 let n = 8 17 # ^ Taylor will converge really quickly 18 # notice that the factorial of 17 is 19 # already pretty gigantic 20 var acc = 0.0 21 for i in 0..n: 22 var k = 2*i + 1 23 var taylor = pow(-1, i.float) * pow(x, k.float) / factorial(k).float 24 acc = acc + taylor 25 return acc 26 27 ## Log function 28 ## <https://en.wikipedia.org/wiki/Natural_logarithm#High_precision> 29 30 ## Arithmetic-geomtric mean 31 proc ag(x: float, y: float): float = 32 let n = 32 # just some high number 33 var a = (x + y)/2.0 34 var b = sqrt(x * y) 35 for i in 0..n: 36 let temp = a 37 a = (a+b)/2.0 38 b = sqrt(b*temp) 39 return a 40 41 ## Find m such that x * 2^m > 2^precision/2 42 proc find_m(x:float): float = 43 var m = 0.0; 44 let precision = 64 # bits 45 let c = pow(2.0, precision.float / 2.0) 46 while x * pow(2.0, m) < c: 47 m = m + 1 48 return m 49 50 proc log(x: float): float = 51 let m = find_m(x) 52 let s = x * pow(2.0, m) 53 let ln2 = 0.6931471805599453 54 return ( PI / (2.0 * ag(1, 4.0/s)) ) - m * ln2 55 56 ## Test these functions 57 ## echo factorial(5) 58 ## echo sine(1.0) 59 ## echo log(0.1) 60 ## echo log(2.0) 61 ## echo log(3.0) 62 ## echo pow(2.0, 32.float) 63 64 ## Distribution functions 65 66 ## Normal 67 ## <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform#Basic_form> 68 proc ur_normal(): float = 69 let u1 = rand(1.0) 70 let u2 = rand(1.0) 71 let z = sqrt(-2.0 * log(u1)) * sine(2 * PI * u2) 72 return z 73 74 proc normal(mean: float, sigma: float): float = 75 return (mean + sigma * ur_normal()) 76 77 proc lognormal(logmean: float, logsigma: float): float = 78 let answer = pow(E, normal(logmean, logsigma)) 79 return answer 80 81 proc to(low: float, high: float): float = 82 let normal95confidencePoint = 1.6448536269514722 83 let loglow = log(low) 84 let loghigh = log(high) 85 let logmean = (loglow + loghigh)/2 86 let logsigma = (loghigh - loglow) / (2.0 * normal95confidencePoint); 87 return lognormal(logmean, logsigma) 88 89 ## echo ur_normal() 90 ## echo normal(10, 20) 91 ## echo lognormal(2, 4) 92 ## echo to(10, 90) 93 94 ## Manipulate samples 95 96 proc make_samples(f: () -> float, n: int): seq[float] = 97 result = toSeq(1..n).map(_ => f()) 98 return result 99 100 proc mixture(sxs: seq[seq[float]], ps: seq[float], n: int): seq[float] = 101 102 assert sxs.len == ps.len 103 104 var ws: seq[float] 105 var sum = 0.0 106 for i, p in ps: 107 sum = sum + p 108 ws.add(sum) 109 ws = ws.map(w => w/sum) 110 111 proc get_mixture_sample(): float = 112 let r = rand(1.0) 113 var i = ws.len - 1 114 for j, w in ws: 115 if r < w: 116 i = j 117 break 118 ## only occasion when ^ doesn't assign i 119 ## is when r is exactly 1 120 ## which would correspond to choosing the last item in ws 121 ## which is why i is initialized to ws.len 122 let xs = sxs[i] 123 let l = xs.len-1 124 let k = rand(0..l) 125 return xs[k] 126 127 return toSeq(1..n).map(_ => get_mixture_sample()) 128 129 ## Actual model 130 131 let n = 1000000 132 133 let p_a = 0.8 134 let p_b = 0.5 135 let p_c = p_a * p_b 136 137 let weights = @[ 1.0 - p_c, p_c/2.0, p_c/4.0, p_c/4.0 ] 138 139 let fs = [ () => 0.0, () => 1.0, () => to(1.0, 3.0), () => to(2.0, 10.0) ] 140 let dists = fs.map(f => make_samples(f, n)) 141 let result = mixture(dists, weights, n) 142 let mean_result = foldl(result, a + b, 0.0) / float(result.len) 143 144 # echo result 145 echo mean_result