Line data Source code
1 : /**
2 : * INFO: a collection of log probability density functions (PDF)
3 : */
4 :
5 : // A number of log probability density functions (PDF). Naming and parameterization
6 : // should match R's, except for that all functions reside in an ld object (
7 : // as in "log density"), so to get a normal log density you would write
8 : // ld.norm(...).
9 : // Most of the code below is directly taken from the great Jstat project
10 : // (https://github.com/jstat/) which includes PDF for many common probability
11 : // distributions. What I have done is only to convert these to log PDFs.
12 :
13 : /*
14 : Original work Copyright (c) 2013 jStat
15 : Modified work Copyright (c) 2015 Rasmus Bååth
16 : Modified work Copyright (c) 2020 Rui Vieira
17 :
18 : Permission is hereby granted, free of charge, to any person obtaining a copy
19 : of this software and associated documentation files (the "Software"), to deal
20 : in the Software without restriction, including without limitation the rights
21 : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
22 : copies of the Software, and to permit persons to whom the Software is
23 : furnished to do so, subject to the following conditions:
24 :
25 : The above copyright notice and this permission notice shall be included in
26 : all copies or substantial portions of the Software.
27 :
28 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
29 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
30 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
31 : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
32 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
33 : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
34 : THE SOFTWARE.
35 :
36 : */
37 1 : import { range } from "../utils.ts";
38 :
39 : interface UnivariateContinuousDistribution {
40 : sample(): number;
41 : sampleN(n: number): Array<number>;
42 : }
43 :
44 : /** Returns a random real number from a normal distribbution defined
45 : * by mean and sd.
46 : * Adapted from https://github.com/jstat/jstat/blob/master/src/special.js */
47 :
48 1 : export function rnorm(mean: number, sd: number) {
49 23001 : var u, v, x, y, q;
50 23001 : do {
51 54345 : u = Math.random();
52 54345 : v = 1.7156 * (Math.random() - 0.5);
53 54345 : x = u - 0.449871;
54 54345 : y = Math.abs(v) + 0.386595;
55 54345 : q = x * x + y * (0.19600 * y - 0.25472 * x);
56 23001 : } while (q > 0.27597 && (q > 0.27846 || v * v > -4 * Math.log(u) * u * u));
57 :
58 23001 : return (v / u) * sd + mean;
59 1 : }
60 :
61 : // Returns the error function erf(x)
62 0 : function erf(x: number): number {
63 0 : const cof = [
64 0 : -1.3026537197817094,
65 0 : 6.4196979235649026e-1,
66 0 : 1.9476473204185836e-2,
67 0 : -9.561514786808631e-3,
68 0 : -9.46595344482036e-4,
69 0 : 3.66839497852761e-4,
70 0 : 4.2523324806907e-5,
71 0 : -2.0278578112534e-5,
72 0 : -1.624290004647e-6,
73 0 : 1.303655835580e-6,
74 0 : 1.5626441722e-8,
75 0 : -8.5238095915e-8,
76 0 : 6.529054439e-9,
77 0 : 5.059343495e-9,
78 0 : -9.91364156e-10,
79 0 : -2.27365122e-10,
80 0 : 9.6467911e-11,
81 0 : 2.394038e-12,
82 0 : -6.886027e-12,
83 0 : 8.94487e-13,
84 0 : 3.13092e-13,
85 0 : -1.12708e-13,
86 0 : 3.81e-16,
87 0 : 7.106e-15,
88 0 : -1.523e-15,
89 0 : -9.4e-17,
90 0 : 1.21e-16,
91 0 : -2.8e-17,
92 0 : ];
93 0 : var j = cof.length - 1;
94 0 : var isneg = false;
95 0 : var d = 0;
96 0 : var dd = 0;
97 0 : var t, ty, tmp, res;
98 :
99 0 : if (x < 0) {
100 0 : x = -x;
101 0 : isneg = true;
102 0 : }
103 :
104 0 : t = 2 / (2 + x);
105 0 : ty = 4 * t - 2;
106 :
107 0 : for (; j > 0; j--) {
108 0 : tmp = d;
109 0 : d = ty * d - dd + cof[j];
110 0 : dd = tmp;
111 0 : }
112 :
113 0 : res = t * Math.exp(-x * x + 0.5 * (cof[0] + ty * d) - dd);
114 0 : return isneg ? res - 1 : 1 - res;
115 0 : }
116 :
117 : // Returns the complmentary error function erfc(x)
118 0 : function erfc(x: number): number {
119 0 : return 1 - erf(x);
120 0 : }// Returns the inverse of the complementary error function
121 :
122 0 : function erfcinv(p: number): number {
123 0 : var j = 0;
124 0 : var x, err, t, pp;
125 0 : if (p >= 2) {
126 0 : return -100;
127 0 : }
128 0 : if (p <= 0) {
129 0 : return 100;
130 0 : }
131 0 : pp = (p < 1) ? p : 2 - p;
132 0 : t = Math.sqrt(-2 * Math.log(pp / 2));
133 0 : x = -0.70711 * ((2.30753 + t * 0.27061) /
134 0 : (1 + t * (0.99229 + t * 0.04481)) - t);
135 0 : for (; j < 2; j++) {
136 0 : err = erfc(x) - pp;
137 0 : x += err / (1.12837916709551257 * Math.exp(-x * x) - x * err);
138 0 : }
139 0 : return (p < 1) ? x : -x;
140 0 : }
141 :
142 : // Distribution
143 1 : export class Normal implements UnivariateContinuousDistribution {
144 : mean: number;
145 : std: number;
146 1 : constructor(mean: number, std: number) {
147 3 : this.mean = mean;
148 3 : this.std = std;
149 1 : }
150 0 : pdf(x: number): number {
151 0 : return Math.exp(
152 0 : -0.5 * Math.log(2 * Math.PI) -
153 0 : Math.log(this.std) -
154 0 : Math.pow(x - this.mean, 2) / (2 * this.std * this.std),
155 0 : );
156 0 : }
157 0 : cdf(x: number): number {
158 0 : return 0.5 *
159 0 : (1 + erf((x - this.mean) / Math.sqrt(2 * this.std * this.std)));
160 0 : }
161 :
162 0 : inv(p: number): number {
163 0 : return -1.41421356237309505 * this.std * erfcinv(2 * p) + this.mean;
164 0 : }
165 1 : sample(): number {
166 11001 : return rnorm(this.mean, this.std);
167 1 : }
168 1 : sampleN(n: number): Array<number> {
169 3 : return range(n).map((_) => this.sample());
170 1 : }
171 0 : variance() {
172 0 : return this.std * this.std;
173 0 : }
174 1 : }
175 :
176 : ////////// Helper functions //////////
177 : //////////////////////////////////////
178 :
179 0 : export function lgamma(x: number) {
180 0 : var j = 0;
181 0 : var cof = [
182 0 : 76.18009172947146,
183 0 : -86.50532032941677,
184 0 : 24.01409824083091,
185 0 : -1.231739572450155,
186 0 : 0.1208650973866179e-2,
187 0 : -0.5395239384953e-5,
188 0 : ];
189 0 : var ser = 1.000000000190015;
190 0 : var xx, y, tmp;
191 0 : tmp = (y = xx = x) + 5.5;
192 0 : tmp -= (xx + 0.5) * log(tmp);
193 0 : for (; j < 6; j++) {
194 0 : ser += cof[j] / ++y;
195 0 : }
196 0 : return log(2.5066282746310005 * ser / xx) - tmp;
197 0 : }
198 :
199 0 : export function lfactorial(n: number) {
200 0 : return n < 0 ? NaN : lgamma(n + 1);
201 0 : }
202 :
203 0 : export function lchoose(n: number, k: number) {
204 0 : return lfactorial(n) - lfactorial(k) - lfactorial(n - k);
205 0 : }
206 :
207 0 : export function lbeta(a: number, b: number) {
208 0 : return lgamma(a) + lgamma(b) - lgamma(a + b);
209 0 : }
210 :
211 1 : let log = Math.log;
212 1 : let abs = Math.abs;
213 1 : let pow = Math.pow;
214 1 : let sqrt = Math.sqrt;
215 1 : let pi = Math.PI;
216 :
217 : ////////// Continous distributions //////////
218 : /////////////////////////////////////////////
219 :
220 0 : export function beta(x: number, shape1: number, shape2: number) {
221 0 : if (x > 1 || x < 0) {
222 0 : return -Infinity;
223 0 : }
224 0 : if (shape1 === 1 && shape2 === 1) {
225 0 : return 0;
226 0 : } else {
227 0 : return (shape1 - 1) * log(x) + (shape2 - 1) * log(1 - x) -
228 0 : lbeta(shape1, shape2);
229 0 : }
230 0 : }
231 :
232 0 : export function cauchy(x: number, location: number, scale: number) {
233 0 : return log(scale) - log(pow(x - location, 2) + pow(scale, 2)) - log(pi);
234 0 : }
235 :
236 1 : export function norm(x: number, mean: number, sd: number) {
237 258292 : return -0.5 * log(2 * pi) - log(sd) - pow(x - mean, 2) / (2 * sd * sd);
238 1 : } // A bivariate Normal distribution parameterized by arrays of two means and SDs, and
239 : // the correlation.
240 :
241 0 : export function bivarnorm(x: any, mean: any, sd: any, corr: number) {
242 0 : var z = pow(x[0] - mean[0], 2) / pow(sd[0], 2) +
243 0 : pow(x[1] - mean[1], 2) / pow(sd[1], 2) -
244 0 : (2 * corr * (x[0] - mean[0]) * (x[1] - mean[1])) / (sd[0] * sd[1]);
245 0 : var normalizing_factor = -(log(2) + log(pi) + log(sd[0]) + log(sd[1]) +
246 0 : 0.5 * log(1 - pow(corr, 2)));
247 0 : var bivar_log_dens = normalizing_factor - z / (2 * (1 - pow(corr, 2)));
248 0 : return bivar_log_dens;
249 0 : }
250 :
251 0 : export function laplace(x: number, location: number, scale: number) {
252 0 : return (-abs(x - location) / scale) - log(2 * scale);
253 0 : }
254 :
255 1 : let dexp = laplace;
256 :
257 0 : export function gamma(x: number, shape: number, rate: number) {
258 0 : var scale = 1 / rate;
259 0 : if (x < 0) {
260 0 : return -Infinity;
261 0 : }
262 0 : if ((x === 0 && shape === 1)) {
263 0 : return -log(scale);
264 0 : } else {
265 0 : return (shape - 1) * log(x) - x / scale - lgamma(shape) -
266 0 : shape * log(scale);
267 0 : }
268 0 : }
269 :
270 0 : export function invgamma(x: number, shape: number, scale: number) {
271 0 : if (x <= 0) {
272 0 : return -Infinity;
273 0 : }
274 0 : return -(shape + 1) * log(x) - scale / x - lgamma(shape) + shape * log(scale);
275 0 : }
276 :
277 0 : export function lnorm(x: number, meanlog: number, sdlog: number) {
278 0 : if (x <= 0) {
279 0 : return -Infinity;
280 0 : }
281 0 : return -log(x) - 0.5 * log(2 * pi) - log(sdlog) -
282 0 : pow(log(x) - meanlog, 2) / (2 * sdlog * sdlog);
283 0 : }
284 :
285 0 : export function pareto(x: number, scale: number, shape: number) {
286 0 : if (x < scale) {
287 0 : return -Infinity;
288 0 : }
289 0 : return log(shape) + shape * log(scale) - (shape + 1) * log(x);
290 0 : }
291 :
292 0 : export function t(x: number, location: number, scale: number, df: number) {
293 0 : df = df > 1e100 ? 1e100 : df;
294 0 : return lgamma((df + 1) / 2) - lgamma(df / 2) - log(sqrt(pi * df) * scale) +
295 0 : log(pow(1 + (1 / df) * pow((x - location) / scale, 2), -(df + 1) / 2));
296 0 : } // This is a direct javascript translation of the R code used to evaluate
297 : // the log density of a weibull distribution:
298 : // https://github.com/wch/r-source/blob/b156e3a711967f58131e23c1b1dc1ea90e2f0c43/src/nmath/dweibull.c
299 :
300 0 : export function weibull(x: number, shape: number, scale: number) {
301 0 : if (x < 0) return -Infinity;
302 0 : if (x === 0 && shape < 1) return Infinity;
303 0 : var tmp1 = pow(x / scale, shape - 1);
304 0 : var tmp2 = tmp1 * (x / scale);
305 0 : return -tmp2 + log(shape * tmp1 / scale);
306 0 : } // This is a direct javascript translation of the R code used to evaluate
307 : // the log density of a logistic distribution:
308 : // https://github.com/wch/r-source/blob/b156e3a711967f58131e23c1b1dc1ea90e2f0c43/src/nmath/dlogis.c
309 :
310 0 : export function logis(x: number, location: number, scale: number) {
311 0 : x = abs((x - location) / scale);
312 0 : var e = Math.exp(-x);
313 0 : var f = 1.0 + e;
314 0 : return -(x + log(scale * f * f));
315 0 : }
316 :
317 0 : export function dirichlet(x: any, alpha: any) {
318 0 : var sum_alpha = 0;
319 0 : var sum_lgamma_alpha = 0;
320 0 : var sum_alpha_sub_1_log_x = 0;
321 0 : var n = alpha.length;
322 0 : for (var i = 0; i < n; i++) {
323 0 : sum_alpha += alpha[i];
324 0 : sum_lgamma_alpha += lgamma(alpha[i]);
325 0 : sum_alpha_sub_1_log_x += (alpha[i] - 1) * log(x[i]);
326 0 : }
327 0 : return lgamma(sum_alpha) - sum_lgamma_alpha + sum_alpha_sub_1_log_x;
328 0 : }
329 :
330 0 : export function exp(x: number, rate: number) {
331 0 : return x < 0 ? -Infinity : log(rate) - rate * x;
332 0 : }
333 :
334 1 : export function unif(x: number, min: number, max: number) {
335 0 : return (x < min || x > max) ? -Infinity : log(1 / (max - min));
336 1 : } ////////// Discrete distributions //////////
337 : ////////////////////////////////////////////
338 :
339 0 : export function bern(x: number, prob: number) {
340 0 : return !(x === 0 || x === 1)
341 0 : ? -Infinity
342 0 : : log(x * prob + (1 - x) * (1 - prob));
343 0 : }
344 :
345 0 : export function cat(x: number, probs: any) {
346 0 : if (x < 1 || x > probs.length) {
347 0 : return -Infinity;
348 0 : } else {
349 0 : return log(probs[x - 1]);
350 0 : }
351 0 : }
352 :
353 0 : export function binom(x: number, size: number, prob: number) {
354 0 : if (x > size || x < 0) {
355 0 : return -Infinity;
356 0 : }
357 0 : if (prob === 0 || prob === 1) {
358 0 : return (size * prob) === x ? 0 : -Infinity;
359 0 : }
360 0 : return lchoose(size, x) + x * log(prob) + (size - x) * log(1 - prob);
361 0 : }
362 :
363 0 : export function ultinom(x: any, probs: any) {
364 0 : var n = x.length;
365 0 : var size = 0;
366 0 : var tmp_term = 0;
367 0 : for (var i = 0; i < n; i++) {
368 0 : if (probs[i] === 0) {
369 0 : if (x[i] !== 0) {
370 0 : return -Infinity;
371 0 : }
372 0 : } else {
373 0 : size += x[i];
374 0 : tmp_term += x[i] * log(probs[i]) - lgamma(x[i] + 1);
375 0 : }
376 0 : }
377 0 : return lgamma(size + 1) + tmp_term;
378 0 : }
379 :
380 0 : export function nbinom(x: number, size: number, prob: number) {
381 0 : if (x < 0) {
382 0 : return -Infinity;
383 0 : }
384 0 : return lchoose(x + size - 1, size - 1) + x * log(1 - prob) + size * log(prob);
385 0 : }
386 :
387 0 : export function hyper(x: number, m: number, n: number, k: number) {
388 0 : if (x < 0 || x > k) {
389 0 : return -Infinity;
390 0 : } else {
391 0 : return lchoose(m, x) + lchoose(n, k - x) - lchoose(m + n, k);
392 0 : }
393 0 : }
394 :
395 0 : export function pois(x: number, lambda: number) {
396 0 : return x < 0 ? -Infinity : log(lambda) * x - lambda - lfactorial(x);
397 0 : }
|
