Line data Source code
1 : /**
2 : * INFO: a port of `bayes.js` to Deno
3 : */
4 1 : import { rnorm } from "../mentat/stats/distributions.ts";
5 : ////////// Helper Functions //////////
6 : //////////////////////////////////////
7 : /** Returns a random real number between min and max */
8 0 : export function runif(min: number, max: number) {
9 0 : return Math.random() * (max - min) + min;
10 0 : } /** Returns a random integer between min and max */
11 :
12 0 : export function runif_discrete(min: number, max: number) {
13 0 : return Math.floor(Math.random() * (max - min + 1)) + min;
14 0 : } /** Returns a deep clone of src, sort of... It only copies a limited
15 : * number of types and, for example, function are not copied.
16 : * From http://davidwalsh.name/javascript-clone
17 : */
18 :
19 1 : export function deep_clone(src: any): any {
20 7 : function mixin(dest: any, source: any, copyFunc: Function): any {
21 10 : var name, s, i, empty = Object.create(null);
22 10 : for (name in source) {
23 : // the (!(name in empty) || empty[name] !== s) condition avoids copying properties in "source"
24 : // inherited from Object.prototype. For example, if dest has a custom toString() method,
25 : // don't overwrite it with the toString() method that source inherited from Object.prototype
26 15 : s = source[name];
27 0 : if (
28 0 : !(name in dest) ||
29 0 : (dest[name] !== s && (!(name in empty) || empty[name] !== s))
30 0 : ) {
31 0 : dest[name] = copyFunc ? copyFunc(s) : s;
32 15 : }
33 10 : }
34 10 : return dest;
35 7 : }
36 7 : if (
37 7 : !src || typeof src != "object" ||
38 7 : Object.prototype.toString.call(src) === "[object Function]"
39 7 : ) {
40 : // null, undefined, any non-object, or function
41 10 : return src; // anything
42 10 : }
43 0 : if (src.nodeType && "cloneNode" in src) {
44 : // DOM Node
45 0 : return src.cloneNode(true); // Node
46 0 : }
47 0 : if (src instanceof Date) {
48 : // Date
49 0 : return new Date(src.getTime()); // Date
50 0 : }
51 0 : if (src instanceof RegExp) {
52 : // RegExp
53 0 : return new RegExp(src); // RegExp
54 0 : }
55 10 : var r, i, l;
56 0 : if (src instanceof Array) {
57 : // array
58 0 : r = [];
59 0 : for (i = 0, l = src.length; i < l; ++i) {
60 0 : if (i in src) {
61 0 : r.push(deep_clone(src[i]));
62 0 : }
63 0 : }
64 0 : } else {
65 : // generic objects
66 0 : r = src.constructor ? new src.constructor() : {};
67 7 : }
68 10 : return mixin(r, src, deep_clone);
69 1 : } /** Specialized clone function that only clones scalars and nested arrays where
70 : * each array either consists of all arrays or all numbers. This function
71 : * is meant as a fast way of cloning parameter draws within the mcmc sampling
72 : * loop.
73 : */
74 :
75 10001 : function clone_param_draw(x: Array<any>): Array<any> {
76 0 : if (Array.isArray(x)) {
77 0 : if (Array.isArray(x[0])) {
78 : // x is an array of arrays so we need to clone it recursively
79 0 : var x_copy = [];
80 0 : for (var i = 0, length = x.length; i < length; i++) {
81 0 : x_copy.push(clone_param_draw(x[i]));
82 0 : }
83 0 : return x_copy;
84 0 : } else { // We'll assume x is a arrays of scalars
85 0 : return x.slice(0);
86 0 : }
87 0 : } else { // We'll assume x is a scalar
88 10001 : return x;
89 10001 : }
90 1 : } /** Returns true if object is a number.
91 : */
92 :
93 3 : function is_number(object: any): boolean {
94 3 : return typeof object == "number" ||
95 3 : (typeof object == "object" && object.constructor === Number);
96 1 : } /**
97 : * Creates and initializes a (possibly multidimensional/nested) array.
98 : * @param dim - An array giving the dimension of the array. For example,
99 : * [5] would yield a 5 element array, and [3,3] would yield a 3 by 3 matrix.
100 : * @param init - A value or a function used to fill in the each element in
101 : * the array. If it is a function it should take no arguments, it will be
102 : * evaluated once for each element, and it's return value will be used to
103 : * fill in each element.
104 : * @example
105 : * // The following would return [[1,1],[1,1],[1,1]]
106 : * create_array([2,3], 1)
107 : */
108 :
109 0 : function create_array(dim: any, init: any): Array<any> {
110 0 : var arr = new Array(dim[0]);
111 0 : var i;
112 0 : if (dim.length == 1) { // Fill it up with init
113 0 : if (typeof init === "function") {
114 0 : for (i = 0; i < dim[0]; i++) {
115 0 : arr[i] = init();
116 0 : }
117 0 : } else {
118 0 : for (i = 0; i < dim[0]; i++) {
119 0 : arr[i] = init;
120 0 : }
121 0 : }
122 0 : } else if (dim.length > 1) {
123 0 : for (i = 0; i < dim[0]; i++) {
124 0 : arr[i] = create_array(dim.slice(1), init);
125 0 : }
126 0 : } else {
127 0 : throw "create_array can't create a dimensionless array";
128 0 : }
129 0 : return arr;
130 0 : } /**
131 : * Return the dimensions of a possibly nested array as an array. For
132 : * example, array_dim( [[1, 2], [1, 2]] ) should return [2, 2]
133 : * Assumes that all arrays inside another array are of the same length.
134 : * @example
135 : * // Should return [4, 2, 1]
136 : * array_dim(create_array([4, 2, 1], 0))
137 : */
138 :
139 0 : function array_dim(a: Array<any>): Array<number> {
140 0 : if (Array.isArray(a[0])) {
141 0 : return [a.length].concat(array_dim(a[0]));
142 0 : } else {
143 0 : return [a.length];
144 0 : }
145 0 : } /**
146 : * Checks if two arrays are equal in the sense that they contain the same elements
147 : * as judged by the "==" operator. Returns true or false.
148 : * Adapted from http://stackoverflow.com/a/14853974/1001848
149 : */
150 :
151 7 : function array_equal(a1: Array<any>, a2: Array<any>): boolean {
152 7 : if (a1.length != a2.length) return false;
153 7 : for (var i = 0; i < a1.length; i++) {
154 : // Check if we have nested arrays
155 0 : if (Array.isArray(a1[i]) && Array.isArray(a2[i])) {
156 : // recurse into the nested arrays
157 0 : if (!array_equal(a1[i], a2[i])) return false;
158 0 : } else if (a1[i] != a2[i]) {
159 : // Warning - two different object instances will never be equal: {x:20} != {x:20}
160 0 : return false;
161 0 : }
162 7 : }
163 7 : return true;
164 1 : } /**
165 : * Traverses a possibly nested array a and applies fun to all "leaf nodes",
166 : * that is, values that are not arrays. Returns an array of the same size as
167 : * a.
168 : */
169 :
170 0 : function nested_array_apply(a: Array<any>, fun: Function) {
171 0 : if (Array.isArray(a)) {
172 0 : var result = new Array(a.length);
173 0 : for (var i = 0; i < a.length; i++) {
174 0 : result[i] = nested_array_apply(a[i], fun);
175 0 : }
176 0 : return result;
177 0 : } else {
178 0 : return fun(a);
179 0 : }
180 0 : } /** Randomizing the array element order in-place. Using Durstenfeld
181 : * shuffle algorithm. Adapted from here:
182 : * http://stackoverflow.com/a/12646864/1001848
183 : */
184 :
185 12001 : function shuffle_array(array: Array<any>): Array<any> {
186 12001 : for (var i = array.length - 1; i > 0; i--) {
187 18001 : var j = Math.floor(Math.random() * (i + 1));
188 18001 : var temp = array[i];
189 18001 : array[i] = array[j];
190 18001 : array[j] = temp;
191 12001 : }
192 12001 : return array;
193 1 : }
194 :
195 : /**
196 : * Does the same thing as nested_array_apply, that is, traverses a possibly
197 : * nested array a and applies fun to all "leaf nodes" and returns an array
198 : * of the same size as a. The difference is that nested_array_random_apply
199 : * branches randomly.
200 : */
201 0 : function nested_array_random_apply(a: Array<any>, fun: Function) {
202 0 : if (Array.isArray(a)) {
203 0 : var len = a.length;
204 0 : var i;
205 0 : var array_is = [];
206 0 : for (i = 0; i < len; i++) {
207 0 : array_is[i] = i;
208 0 : }
209 0 : shuffle_array(array_is);
210 0 : var result = [];
211 :
212 0 : for (i = 0; i < len; i++) {
213 0 : var array_i = array_is[i];
214 0 : result[array_i] = nested_array_apply(a[array_i], fun);
215 0 : }
216 0 : return result;
217 0 : } else {
218 0 : return fun(a);
219 0 : }
220 0 : } /**
221 : * Allows a pretty way of setting default options where the defults can be
222 : * overridden by an options object.
223 : * @param option_name - the name of the option as a string
224 : * @param my_options - an option object that could have option_name
225 : * as a member.
226 : * @param defaul_value - defult value that is returned if option_name
227 : * is not defined in my_options.
228 : * @example
229 : * var my_options = {pi: 3.14159}
230 : * var pi = get_option("pi", my_options, 3.14)
231 : */
232 :
233 : // Pretty way of setting default options where the defaults can be overridden
234 : // by an options object. For example:
235 : // var pi = get_option("pi", my_options, 3.14)
236 :
237 16 : function get_option(option_name: string, options: any, defaul_value: any) {
238 16 : options = options || {};
239 16 : return options.hasOwnProperty(option_name) &&
240 0 : options[option_name] !== undefined &&
241 0 : options[option_name] !== null
242 0 : ? options[option_name]
243 16 : : defaul_value;
244 1 : } /** Version of get_option where the option should be a one or multi-dimensional
245 : * array and where the default can be overridden either by a scalar or by an array.
246 : * If it's a scalar the that scalar is used to initialize an array with
247 : * dim dimensions.
248 : *
249 : */
250 :
251 0 : var get_multidim_option = function (
252 0 : option_name: string,
253 0 : options: any,
254 0 : dim: Array<number>,
255 0 : defaul_value: any,
256 : ) {
257 0 : var value = get_option(option_name, options, defaul_value);
258 0 : if (!Array.isArray(value)) {
259 0 : value = create_array(dim, value);
260 0 : }
261 0 : if (!array_equal(array_dim(value), dim)) {
262 0 : throw "The option " + option_name + " is of dimension [" +
263 0 : array_dim(value) + "] but should be [" + dim + "].";
264 0 : }
265 0 : return value;
266 0 : };
267 :
268 : ////////// Functions for handling parameter objects //////////
269 : //////////////////////////////////////////////////////////////
270 :
271 : /**
272 : * Returns a fixed (same every time) number that could be used to initialize
273 : * a parameter of a certain type, possibly with lower and upper bounds.
274 : * The possile types are "real", "int", and "binary".
275 : */
276 1 : var param_init_fixed = function (
277 1 : type: ParameterType,
278 1 : lower: number,
279 1 : upper: number,
280 : ) {
281 0 : if (lower > upper) {
282 0 : throw "Can not initialize parameter where lower bound > upper bound";
283 0 : }
284 3 : if (type == ParameterType.Real) {
285 3 : if (lower === -Infinity && upper === Infinity) {
286 4 : return 0.5;
287 4 : } else if (lower === -Infinity) {
288 0 : return upper - 0.5;
289 0 : } else if (upper === Infinity) {
290 4 : return lower + 0.5;
291 0 : } else if (lower <= upper) {
292 0 : return (lower + upper) / 2;
293 0 : }
294 0 : } else if (type == ParameterType.Int) {
295 0 : if (lower === -Infinity && upper === Infinity) {
296 0 : return 1;
297 0 : } else if (lower === -Infinity) {
298 0 : return upper - 1;
299 0 : } else if (upper === Infinity) {
300 0 : return lower + 1;
301 0 : } else if (lower <= upper) {
302 0 : return Math.round((lower + upper) / 2);
303 0 : }
304 0 : } else if (type == ParameterType.Binary) {
305 0 : return 1;
306 0 : }
307 0 : throw "Could not initialize parameter of type " + type + "[" + lower + ", " +
308 0 : upper + "]";
309 1 : };
310 :
311 : /**
312 : * Completes params_to_complete, an object containing parameter descriptions,
313 : * and initializes non-initialized parameters. This modified version of
314 : * params_to_complete is returned as a deep copy and not modified in place.
315 : * Initialization is done by supplying a param_init function with signature
316 : * function(type, lower, upper) that should return a single number
317 : * (like param_init_fixed, for example).
318 : * @example
319 : * var params = { "mu": {"type": "real"} }
320 : * params = complete_params(params);
321 : * // params should now be:
322 : * // {"mu": { "type": "real", "dim": [1], "upper": Infinity,
323 : * // "lower": -Infinity, "init": 0.5 }}
324 : */
325 1 : var complete_params = function (params_to_complete: any, param_init: any) {
326 2 : var params = deep_clone(params_to_complete);
327 2 : for (var param_name in params) {
328 4 : if (!params.hasOwnProperty(param_name)) continue;
329 4 : var param = params[param_name];
330 0 : if (!param.hasOwnProperty("type")) {
331 0 : param.type = ParameterType.Real;
332 0 : }
333 4 : if (!param.hasOwnProperty("dim")) {
334 4 : param.dim = [1];
335 4 : }
336 0 : if (is_number(param.dim)) {
337 0 : param.dim = [param.dim];
338 0 : }
339 0 : if (param.type == ParameterType.Binary) {
340 0 : param.upper = 1;
341 0 : param.lower = 0;
342 0 : }
343 4 : if (!param.hasOwnProperty("upper")) {
344 4 : param.upper = Infinity;
345 4 : }
346 4 : if (!param.hasOwnProperty("lower")) {
347 5 : param.lower = -Infinity;
348 4 : }
349 :
350 0 : if (param.hasOwnProperty("init")) {
351 : // If this is just a number or a nested array we leave it alone, but if...
352 0 : if (array_equal(param.dim, [1]) && typeof param.init === "function") {
353 : // param.init is a function, use that to initialize the parameter.
354 0 : param.init = param.init();
355 0 : } else if (!array_equal(param.dim, [1]) && !Array.isArray(param.init)) {
356 : // We have a multidimensional parameter where the param.init exist but
357 : // is not an array. Then assume it is a number or a function and use
358 : // it to initialize the parameter.
359 0 : param.init = create_array(param.dim, param.init);
360 0 : }
361 0 : } else { // We use the default initialization function.
362 4 : if (array_equal(param.dim, [1])) {
363 4 : param.init = param_init(param.type, param.lower, param.upper);
364 0 : } else {
365 0 : param.init = create_array(param.dim, function () {
366 0 : return param_init(param.type, param.lower, param.upper);
367 0 : });
368 0 : }
369 4 : }
370 2 : }
371 2 : return params;
372 1 : };
373 :
374 : ////////// Stepper Functions ///////////
375 : ////////////////////////////////////////
376 :
377 : /**
378 : * @interface
379 : * A Stepper is an object responsible for pushing around one
380 : * or more parameter values in a state according to the distribution
381 : * defined by the log posterior. This defines the Stepper "interface",
382 : * where "interface" means that Stepper defines a class that is never
383 : * meant to be instantiated, but just to be subclassed by specialized
384 : * stepper functions.
385 : * @interface
386 : * @param params - An object with parameter definitions, for example:
387 : * {"mu": { "type": "real", "dim": [1], "upper": Infinity,
388 : * "lower": -Infinity, "init": 0.5 }}
389 : * The parameter definitions are expected to be "complete", that is,
390 : * specifying all relevant attributes such as dim, lower and upper.
391 : * @param state - an object containing the state of all parameters in params
392 : * (and possibly more). The parameter names are given as keys and the states
393 : * as scalars or, possibly nested, arrays. For example:
394 : * {mu: 10, sigma: 5, beta: [1, 2.5]}
395 : * @param log_post - A function *taking no parameters* that returns the
396 : * log density that depends on the state. That is, the value of log_post
397 : * should change if the the values in state are changed.
398 :
399 : */
400 1 : abstract class Stepper {
401 : params: any;
402 : state: any;
403 : log_post: Function;
404 1 : constructor(params: any, state: any, log_post: Function) {
405 4 : this.params = params;
406 4 : this.state = state;
407 4 : this.log_post = log_post;
408 1 : }
409 : /**
410 : * Takes a step in the parameter space. Should return the new state,
411 : * but is mainly called for it's side effect of making a change in the
412 : * state object.
413 : */
414 : abstract step(): any;
415 :
416 : /**
417 : * If implemented, makes the stepper adapt while stepping.
418 : */
419 0 : start_adaptation() {
420 0 : }
421 : /**
422 : * If implemented, makes the stepper cease adapting while stepping.
423 : */
424 0 : stop_adaptation() {
425 : // Optional, some steppers might not be adaptive. */
426 0 : }
427 : /**
428 : * Returns an object containg info regarding the stepper.
429 : */
430 : abstract info(): any;
431 1 : }
432 :
433 : /**
434 : * @class
435 : * @implements {Stepper}
436 : * Constructor for an object that implements the metropolis step in
437 : * the Adaptive Metropolis-Within-Gibbs algorithm in "Examples of Adaptive MCMC"
438 : * by Roberts and Rosenthal (2008).
439 : * @param params - An object with a single parameter definition.
440 : * @param state - an object containing the state of all parameters.
441 : * @param log_post - A function that returns the log density that depends on the state.
442 : * @param options - an object with options to the stepper.
443 : * @param generate_proposal - a function returning a proposal (as a number)
444 : * with signature function(param_state, log_scale) where param_state is a
445 : * number and log_scale defines the scale of the proposal somehow.
446 : */
447 1 : export class OnedimMetropolisStepper extends Stepper {
448 : param_name: any;
449 : lower: number;
450 : upper: number;
451 : prop_log_scale: any;
452 : batch_size: any;
453 : max_adaptation: any;
454 : initial_adaptation: any;
455 : target_accept_rate: any;
456 : is_adapting: any;
457 : generate_proposal: any;
458 : acceptance_count: any;
459 : batch_count: any;
460 : iterations_since_adaption: any;
461 1 : constructor(
462 1 : params: any,
463 1 : state: any,
464 1 : log_post: Function,
465 1 : options: any,
466 1 : generate_proposal: Function,
467 : ) {
468 3 : super(params, state, log_post);
469 :
470 3 : var param_names = Object.keys(this.params);
471 0 : if (param_names.length != 1) {
472 0 : throw "OnedimMetropolisStepper can only handle one parameter.";
473 0 : }
474 3 : this.param_name = param_names[0];
475 3 : var param = this.params[this.param_name];
476 0 : if (!array_equal(param.dim, [1])) {
477 0 : throw "OnedimMetropolisStepper can only handle one one-dimensional parameter.";
478 0 : }
479 3 : this.lower = param.lower;
480 3 : this.upper = param.upper;
481 :
482 3 : this.prop_log_scale = get_option("prop_log_scale", options, 0);
483 3 : this.batch_size = get_option("batch_size", options, 50);
484 3 : this.max_adaptation = get_option("max_adaptation", options, 0.33);
485 3 : this.initial_adaptation = get_option("initial_adaptation", options, 1.0);
486 3 : this.target_accept_rate = get_option("target_accept_rate", options, 0.44);
487 3 : this.is_adapting = get_option("is_adapting", options, true);
488 :
489 3 : this.generate_proposal = generate_proposal;
490 :
491 3 : this.acceptance_count = 0;
492 3 : this.batch_count = 0;
493 3 : this.iterations_since_adaption = 0;
494 1 : }
495 :
496 1 : step() {
497 12001 : var param_state = this.state[this.param_name];
498 12001 : var param_proposal = this.generate_proposal(
499 12001 : param_state,
500 12001 : this.prop_log_scale,
501 12001 : );
502 12001 : if (param_proposal < this.lower || param_proposal > this.upper) {
503 : // Outside of limits of the parameter, reject the proposal
504 : // and stay at the current state.
505 12001 : } else { // make a Metropolis step
506 23741 : var curr_log_dens = this.log_post();
507 23741 : this.state[this.param_name] = param_proposal;
508 23741 : var prop_log_dens = this.log_post();
509 23741 : var accept_prob = Math.exp(prop_log_dens - curr_log_dens);
510 23741 : if (accept_prob > Math.random()) {
511 : // We do nothing as the state of param has already been changed to the proposal
512 29273 : if (this.is_adapting) this.acceptance_count++;
513 23741 : } else {
514 : // revert state back to the old state of param
515 29949 : this.state[this.param_name] = param_state;
516 23741 : }
517 12001 : }
518 12001 : if (this.is_adapting) {
519 12001 : this.iterations_since_adaption++;
520 12001 : if (this.iterations_since_adaption >= this.batch_size) { // then adapt
521 12241 : this.batch_count++;
522 12241 : var log_sd_adjustment = Math.min(
523 12241 : this.max_adaptation,
524 12241 : this.initial_adaptation / Math.sqrt(this.batch_count),
525 12241 : );
526 12241 : if (this.acceptance_count / this.batch_size > this.target_accept_rate) {
527 12363 : this.prop_log_scale += log_sd_adjustment;
528 12241 : } else {
529 12359 : this.prop_log_scale -= log_sd_adjustment;
530 12241 : }
531 12241 : this.acceptance_count = 0;
532 12241 : this.iterations_since_adaption = 0;
533 12001 : }
534 12001 : }
535 12001 : return this.state[this.param_name];
536 1 : }
537 :
538 0 : start_adaptation() {
539 0 : this.is_adapting = true;
540 0 : }
541 :
542 0 : stop_adaptation() {
543 0 : this.is_adapting = false;
544 0 : }
545 :
546 0 : info() {
547 0 : return {
548 0 : prop_log_scale: this.prop_log_scale,
549 0 : is_adapting: this.is_adapting,
550 0 : acceptance_count: this.acceptance_count,
551 0 : iterations_since_adaption: this.iterations_since_adaption,
552 0 : batch_count: this.batch_count,
553 0 : };
554 0 : }
555 1 : }
556 :
557 : /**
558 : * Function returning a Normal proposal.
559 : */
560 1 : const normal_proposal = function (param_state: any, prop_log_scale: any) {
561 12001 : return rnorm(param_state, Math.exp(prop_log_scale));
562 1 : };
563 :
564 : /**
565 : * @class
566 : * @augments {OnedimMetropolisStepper}
567 : * A "subclass" of OnedimMetropolisStepper making continous Normal proposals.
568 : */
569 1 : class RealMetropolisStepper extends OnedimMetropolisStepper {
570 1 : constructor(params: any, state: any, log_post: Function, options: any) {
571 3 : super(params, state, log_post, options, normal_proposal);
572 1 : }
573 1 : }
574 :
575 : /**
576 : * Function returning a discretized Normal proposal.
577 : */
578 0 : const discrete_normal_proposal = function (
579 0 : param_state: any,
580 0 : prop_log_scale: any,
581 : ) {
582 0 : return Math.round(rnorm(param_state, Math.exp(prop_log_scale)));
583 0 : };
584 :
585 : /**
586 : * @class
587 : * @augments {OnedimMetropolisStepper}
588 : * A "subclass" of OnedimMetropolisStepper making discretized Normal proposals.
589 : */
590 1 : class IntMetropolisStepper extends OnedimMetropolisStepper {
591 0 : constructor(params: any, state: any, log_post: Function, options: any) {
592 0 : super(params, state, log_post, options, discrete_normal_proposal);
593 0 : }
594 1 : }
595 :
596 : /**
597 : * @class
598 : * @implements {Stepper}
599 : * Constructor for an object that implements the metropolis step in
600 : * the Adaptive Metropolis-Within-Gibbs algorithm in "Examples of Adaptive MCMC"
601 : * by Roberts and Rosenthal (2008) for possibly multidimensional arrays. That
602 : * is, instead of just taking a step for a one-dimensional parameter like
603 : * OnedimMetropolisStepper, this Stepper is responsible for taking steps
604 : * for a multidimensional array. It's still pretty dumb and just takes
605 : * one-dimensional steps for each parameter component, though.
606 : * @param params - An object with a single parameter definition for a
607 : * multidimensional parameter.
608 : * @param state - an object containing the state of all parameters.
609 : * @param log_post - A function that returns the log density that depends on the state.
610 : * @param options - an object with options to the stepper.
611 : * @param SubStepper - a constructor for the type of one dimensional Stepper to apply on
612 : * all the components of the multidimensional parameter.
613 : */
614 :
615 1 : export class MultidimComponentMetropolisStepper extends Stepper {
616 : param_name: any;
617 : lower: number;
618 : upper: number;
619 : dim: any;
620 : prop_log_scale: any;
621 : batch_size: any;
622 : max_adaptation: any;
623 : initial_adaptation: any;
624 : target_accept_rate: any;
625 : is_adapting: any;
626 : substeppers: any;
627 0 : constructor(
628 0 : params: any,
629 0 : state: any,
630 0 : log_post: Function,
631 0 : options: any,
632 0 : SubStepper: any,
633 : ) {
634 0 : super(params, state, log_post);
635 :
636 0 : var param_names = Object.keys(this.params);
637 0 : if (param_names.length != 1) {
638 0 : throw "MultidimComponentMetropolisStepper can't handle more than one parameter.";
639 0 : }
640 0 : this.param_name = param_names[0];
641 0 : var param = this.params[this.param_name];
642 0 : this.lower = param.lower;
643 0 : this.upper = param.upper;
644 0 : this.dim = param.dim;
645 :
646 0 : this.prop_log_scale = get_multidim_option(
647 0 : "prop_log_scale",
648 0 : options,
649 0 : this.dim,
650 0 : 0,
651 0 : );
652 0 : this.batch_size = get_multidim_option("batch_size", options, this.dim, 50);
653 0 : this.max_adaptation = get_multidim_option(
654 0 : "max_adaptation",
655 0 : options,
656 0 : this.dim,
657 0 : 0.33,
658 0 : );
659 0 : this.initial_adaptation = get_multidim_option(
660 0 : "initial_adaptation",
661 0 : options,
662 0 : this.dim,
663 0 : 1.0,
664 0 : );
665 0 : this.target_accept_rate = get_multidim_option(
666 0 : "target_accept_rate",
667 0 : options,
668 0 : this.dim,
669 0 : 0.44,
670 0 : );
671 0 : this.is_adapting = get_multidim_option(
672 0 : "is_adapting",
673 0 : options,
674 0 : this.dim,
675 0 : true,
676 0 : );
677 :
678 : // This hack below is a recursive function that creates an array of
679 : // one dimensional steppers according to dim.
680 0 : var create_substeppers = function (
681 0 : dim: Array<number>,
682 0 : substate: any,
683 0 : log_post: Function,
684 0 : prop_log_scale: any,
685 0 : batch_size: any,
686 0 : max_adaptation: any,
687 0 : initial_adaptation: any,
688 0 : target_accept_rate: any,
689 0 : is_adapting: any,
690 : ): any {
691 0 : var substeppers = [];
692 0 : if (dim.length === 1) {
693 0 : for (var i = 0; i < dim[0]; i++) {
694 0 : var suboptions = {
695 0 : prop_log_scale: prop_log_scale[i],
696 0 : batch_size: batch_size[i],
697 0 : max_adaptation: max_adaptation[i],
698 0 : initial_adaptation: initial_adaptation[i],
699 0 : target_accept_rate: target_accept_rate[i],
700 0 : is_adapting: is_adapting[i],
701 0 : };
702 0 : var subparam: any = {};
703 0 : subparam[i] = deep_clone(param);
704 0 : subparam[i].dim = [1]; // As this should now be a one-dim parameter
705 0 : delete subparam[i].init; // As it sould not be needed
706 0 : substeppers[i] = new SubStepper(
707 0 : subparam,
708 0 : substate,
709 0 : log_post,
710 0 : suboptions,
711 0 : );
712 0 : }
713 0 : } else {
714 0 : for (var i = 0; i < dim[0]; i++) {
715 0 : substeppers[i] = create_substeppers(
716 0 : dim.slice(1),
717 0 : substate[i],
718 0 : log_post,
719 0 : prop_log_scale[i],
720 0 : batch_size[i],
721 0 : max_adaptation[i],
722 0 : initial_adaptation[i],
723 0 : target_accept_rate[i],
724 0 : is_adapting[i],
725 0 : );
726 0 : }
727 0 : }
728 0 : return substeppers;
729 0 : };
730 :
731 0 : this.substeppers = create_substeppers(
732 0 : this.dim,
733 0 : this.state[this.param_name],
734 0 : this.log_post,
735 0 : this.prop_log_scale,
736 0 : this.batch_size,
737 0 : this.max_adaptation,
738 0 : this.initial_adaptation,
739 0 : this.target_accept_rate,
740 0 : this.is_adapting,
741 0 : );
742 0 : }
743 :
744 0 : step() {
745 : // Go through the substeppers in a random order and call step() on them.
746 0 : return nested_array_random_apply(
747 0 : this.substeppers,
748 0 : function (substepper: any) {
749 0 : return substepper.step();
750 0 : },
751 0 : );
752 0 : }
753 :
754 0 : start_adaptation() {
755 0 : nested_array_apply(this.substeppers, function (substepper: any) {
756 0 : substepper.start_adaptation();
757 0 : });
758 0 : }
759 :
760 0 : stop_adaptation() {
761 0 : nested_array_apply(this.substeppers, function (substepper: any) {
762 0 : substepper.stop_adaptation();
763 0 : });
764 0 : }
765 :
766 0 : info() {
767 0 : return nested_array_apply(this.substeppers, function (substepper: any) {
768 0 : return substepper.info();
769 0 : });
770 0 : }
771 1 : }
772 :
773 : /**
774 : * @class
775 : * @augments {MultidimComponentMetropolisStepper}
776 : * A "subclass" of MultidimComponentMetropolisStepper making continous Normal proposals.
777 : */
778 1 : class MultiRealComponentMetropolisStepper
779 1 : extends MultidimComponentMetropolisStepper {
780 0 : constructor(params: any, state: any, log_post: any, options: any) {
781 0 : super(params, state, log_post, options, RealMetropolisStepper);
782 0 : }
783 1 : }
784 :
785 : /**
786 : * @class
787 : * @augments {MultidimComponentMetropolisStepper}
788 : * A "subclass" of MultidimComponentMetropolisStepper making discretized Normal proposals.
789 : */
790 1 : class MultiIntComponentMetropolisStepper
791 1 : extends MultidimComponentMetropolisStepper {
792 0 : constructor(params: any, state: any, log_post: any, options: any) {
793 0 : super(params, state, log_post, options, IntMetropolisStepper);
794 0 : }
795 1 : }
796 :
797 : /**
798 : * @class
799 : * @implements {Stepper}
800 : * Constructor for an object that implements a step for a binary parameter.
801 : * This is done by evaluating the log posterior for both states of the
802 : * parameter and then selecting a state randomly with probability relative
803 : * to the posterior of each state.
804 : * @param params - An object with a single parameter definition.
805 : * @param state - an object containing the state of all parameters.
806 : * @param log_post - A function that returns the log density that depends on the state.
807 : * @param options - an object with options to the stepper.
808 : */
809 1 : export class BinaryStepper extends Stepper {
810 : param_name: any;
811 0 : constructor(params: any, state: any, log_post: Function, options?: any) {
812 0 : super(params, state, log_post);
813 0 : var param_names = Object.keys(this.params);
814 0 : if (param_names.length == 1) {
815 0 : this.param_name = param_names[0];
816 0 : } else {
817 0 : throw "BinaryStepper can't handle more than one parameter.";
818 0 : }
819 0 : }
820 0 : step() {
821 0 : this.state[this.param_name] = 0;
822 0 : var zero_log_dens = this.log_post();
823 0 : this.state[this.param_name] = 1;
824 0 : var one_log_dens = this.log_post();
825 0 : var max_log_dens = Math.max(zero_log_dens, one_log_dens);
826 0 : zero_log_dens -= max_log_dens;
827 0 : one_log_dens -= max_log_dens;
828 0 : var zero_prob = Math.exp(
829 0 : zero_log_dens -
830 0 : Math.log(Math.exp(zero_log_dens) + Math.exp(one_log_dens)),
831 0 : );
832 0 : if (Math.random() < zero_prob) {
833 0 : this.state[this.param_name] = 0;
834 0 : return 0;
835 0 : } // else keep the param at 1 .
836 0 : return 1;
837 0 : }
838 0 : info() {
839 0 : return {};
840 0 : }
841 1 : }
842 :
843 : /**
844 : * @class
845 : * @implements {Stepper}
846 : * Just like MultidimComponentMetropolisStepper this Stepper takes a steps for
847 : * a multidimensional parameter by updating each component in turn. The difference
848 : * is that this stepper works on binary parameters.
849 : * @param params - An object with a single parameter definition for a
850 : * multidimensional parameter.
851 : * @param state - an object containing the state of all parameters.
852 : * @param log_post - A function that returns the log density that depends on the state.
853 : * @param options - an object with options to the stepper.
854 : */
855 1 : export class BinaryComponentStepper extends Stepper {
856 : substeppers: any;
857 : param_name: any;
858 : dim: any;
859 0 : constructor(params: any, state: any, log_post: any, options: any) {
860 0 : super(params, state, log_post);
861 : // Stepper.call(this, );
862 :
863 0 : var param_names = Object.keys(this.params);
864 0 : if (param_names.length == 1) {
865 0 : this.param_name = param_names[0];
866 0 : var param = this.params[this.param_name];
867 0 : this.dim = param.dim;
868 0 : } else {
869 0 : throw "BinaryComponentStepper can't handle more than one parameter.";
870 0 : }
871 0 : var create_substeppers = function (
872 0 : dim: Array<number>,
873 0 : substate: any,
874 0 : log_post: Function,
875 : ): any {
876 0 : var substeppers = [];
877 0 : var i;
878 0 : if (dim.length === 1) {
879 0 : for (i = 0; i < dim[0]; i++) {
880 0 : var subparams: any = {};
881 0 : subparams[i] = param;
882 0 : substeppers[i] = new BinaryStepper(subparams, substate, log_post);
883 0 : }
884 0 : } else {
885 0 : for (i = 0; i < dim[0]; i++) {
886 0 : substeppers[i] = create_substeppers(
887 0 : dim.slice(1),
888 0 : substate[i],
889 0 : log_post,
890 0 : );
891 0 : }
892 0 : }
893 0 : return substeppers;
894 0 : };
895 :
896 0 : this.substeppers = create_substeppers(
897 0 : this.dim,
898 0 : this.state[this.param_name],
899 0 : this.log_post,
900 0 : );
901 0 : }
902 :
903 0 : step() {
904 : // Go through the substeppers in a random order and call step() on them.
905 0 : return nested_array_random_apply(
906 0 : this.substeppers,
907 0 : function (substepper: any) {
908 0 : return substepper.step();
909 0 : },
910 0 : );
911 0 : }
912 :
913 0 : info() {
914 0 : return {};
915 0 : }
916 1 : } /**
917 : * @class
918 : * @implements {Stepper}
919 : * This stepper can be responsible for taking a step for one or more parameters.
920 : * For real and int parameters it takes Metropolis within Gibbs steps, and for
921 : * binary parameters it does evaluates the posterior for both paramter values and
922 : * randomly changes to a certain value proportionally to that value's posterior
923 : * (this is also done for each parameter, so also a * within Gibbs approach).
924 : * This stepper is also adaptive and can be efficient when the number of parameters
925 : * are not too high and the correlations between parameters are low.
926 : * @param params - An object with a one or more parameter definitions
927 : * @param state - an object containing the state of all parameters.
928 : * @param log_post - A function that returns the log density that depends on the state.
929 : * @param options - an object with options to the stepper.
930 : */
931 :
932 1 : export enum ParameterType {
933 2 : Real,
934 2 : Binary,
935 2 : Int,
936 1 : }
937 :
938 1 : export class AmwgStepper extends Stepper {
939 : param_names: any;
940 : substeppers: any;
941 1 : constructor(params: any, state: any, log_post: Function, options: any) {
942 2 : super(params, state, log_post);
943 2 : this.param_names = Object.keys(this.params);
944 2 : this.substeppers = [];
945 2 : for (var i = 0; i < this.param_names.length; i++) {
946 4 : var param = params[this.param_names[i]];
947 4 : var SelectStepper;
948 4 : switch (param.type) {
949 4 : case ParameterType.Real:
950 4 : if (array_equal(param.dim, [1])) {
951 4 : SelectStepper = RealMetropolisStepper;
952 0 : } else {
953 0 : SelectStepper = MultiRealComponentMetropolisStepper;
954 0 : }
955 4 : break;
956 0 : case ParameterType.Int:
957 0 : if (array_equal(param.dim, [1])) {
958 0 : SelectStepper = IntMetropolisStepper;
959 0 : } else {
960 0 : SelectStepper = MultiIntComponentMetropolisStepper;
961 0 : }
962 0 : break;
963 0 : case ParameterType.Binary:
964 0 : if (array_equal(param.dim, [1])) {
965 0 : SelectStepper = BinaryStepper;
966 0 : } else {
967 0 : SelectStepper = BinaryComponentStepper;
968 0 : }
969 0 : break;
970 0 : default:
971 0 : throw "AmwgStepper can't handle parameter " + this.param_names[i] +
972 0 : " with type " + param.type;
973 4 : }
974 4 : var param_object_wrap: any = {};
975 4 : param_object_wrap[this.param_names[i]] = param;
976 4 : options = options || {};
977 0 : var param_options =
978 0 : options.params && options.params[this.param_names[i]] || {};
979 4 : param_options.prop_log_scale = param_options.prop_log_scale ||
980 4 : options.prop_log_scale;
981 4 : param_options.batch_size = param_options.batch_size || options.batch_size;
982 4 : param_options.max_adaptation = param_options.max_adaptation ||
983 4 : options.max_adaptation;
984 4 : param_options.initial_adaptation = param_options.initial_adaptation ||
985 4 : options.initial_adaptation;
986 4 : param_options.target_accept_rate = param_options.target_accept_rate ||
987 4 : options.target_accept_rate;
988 4 : param_options.is_adapting = param_options.is_adapting ||
989 4 : options.is_adapting;
990 4 : this.substeppers[i] = new SelectStepper(
991 4 : param_object_wrap,
992 4 : state,
993 4 : log_post,
994 4 : param_options,
995 4 : );
996 2 : }
997 1 : }
998 1 : step() {
999 6001 : shuffle_array(this.substeppers);
1000 6001 : for (var i = 0; i < this.substeppers.length; i++) {
1001 18001 : this.substeppers[i].step();
1002 6001 : }
1003 6001 : return this.state;
1004 1 : }
1005 :
1006 0 : start_adaptation() {
1007 0 : for (var i = 0; i < this.substeppers.length; i++) {
1008 0 : this.substeppers[i].start_adaptation();
1009 0 : }
1010 0 : }
1011 :
1012 0 : stop_adaptation() {
1013 0 : for (var i = 0; i < this.substeppers.length; i++) {
1014 0 : this.substeppers[i].stop_adaptation();
1015 0 : }
1016 0 : }
1017 :
1018 0 : info() {
1019 0 : var info: any = {};
1020 0 : for (var i = 0; i < this.substeppers.length; i++) {
1021 0 : info[this.param_names[i]] = this.substeppers[i].info();
1022 0 : }
1023 0 : return info;
1024 0 : }
1025 1 : }
1026 :
1027 : /////////// Sampler Functions //////////
1028 : ////////////////////////////////////////
1029 :
1030 : /**
1031 : * @interface
1032 : * While you could fit a model by pasting together Steppers, a
1033 : // Sampler is here is a convenience class where an instance of Sampler
1034 : // sets up the Steppers, checks the parameter definition,
1035 : // and manages the sampling. This here defines the Sampler "interface".
1036 : * @interface
1037 : * @param params - An object with parameter definitions, for example:
1038 : * {"mu": {"type": "real"}, "sigma": {"type": "real", "lower" = 0}}
1039 : * The parameter definitions doesn't have to be "complete" and properties
1040 : * left out (like lower and upper) will be filled in by defaults.
1041 : * @param log_post - A function with signature function(state, data). Here
1042 : * state will be an object representing the state with each parameter as a
1043 : * key and the parameter values as numbers or arrays. For example:
1044 : * {"mu": 3, "sigma": 1.5}. The data argument will be the same object as
1045 : * the data argument given below.
1046 : * @param data - an object that will be passed on to the log_post function
1047 : * when sampling.
1048 : * @param options - an object with options to the sampler.
1049 : */
1050 1 : abstract class Sampler {
1051 : params: any;
1052 : data: any;
1053 : param_names: any;
1054 : param_init_fun: any;
1055 : options: any;
1056 : log_post: Function;
1057 : state: any;
1058 : steppers: any;
1059 : monitored_params: any;
1060 : thinning_interval: any;
1061 1 : constructor(params: any, log_post: Function, data: any, options: any) {
1062 2 : this.params = params;
1063 2 : this.data = data;
1064 2 : this.param_names = Object.keys(this.params);
1065 :
1066 : // Setting default options if not passed through the options object
1067 2 : this.param_init_fun = get_option(
1068 2 : "param_init_fun",
1069 2 : options,
1070 2 : param_init_fixed,
1071 2 : );
1072 2 : var thinning_interval = get_option("thin", options, 1);
1073 2 : var params_to_monitor = get_option("monitor", options, null);
1074 2 : this.thin(thinning_interval);
1075 2 : this.monitor(params_to_monitor);
1076 2 : this.options = options;
1077 : // Completing the params and initializing the state.
1078 2 : this.params = complete_params(this.params, this.param_init_fun);
1079 2 : var state: any = {};
1080 2 : for (var i = 0; i < this.param_names.length; i++) {
1081 4 : state[this.param_names[i]] = this.params[this.param_names[i]].init;
1082 2 : }
1083 2 : this.log_post = function () {
1084 23483 : return log_post(state, data);
1085 2 : };
1086 : // Running the log_post function once in case it further modifies the state
1087 : // for example adding derived quantities.
1088 2 : this.log_post();
1089 2 : this.state = state;
1090 2 : this.steppers = this.create_stepper_ensamble(
1091 2 : this.params,
1092 2 : this.state,
1093 2 : this.log_post,
1094 2 : this.options,
1095 2 : );
1096 1 : }
1097 : /** Should return a vector of steppers that when called
1098 : * should take a step in the parameter space.
1099 : */
1100 : abstract create_stepper_ensamble(
1101 : params: any,
1102 : state: any,
1103 : log_post: Function,
1104 : options: any,
1105 : ): any;
1106 :
1107 : /** Returns an object with info about the state of the Sampler.
1108 : */
1109 0 : info() {
1110 0 : return {
1111 0 : state: this.state,
1112 0 : thin: this.thin,
1113 0 : monitor: this.monitor,
1114 0 : steppers: this.steppers,
1115 0 : };
1116 0 : }
1117 :
1118 : /** Takes a step in the parameter space. Returns the new space
1119 : * but also modifies the state in place.
1120 : */
1121 1 : step() {
1122 6001 : shuffle_array(this.steppers);
1123 6001 : for (var i = 0; i < this.steppers.length; i++) {
1124 6001 : this.steppers[i].step();
1125 6001 : }
1126 0 : if (Object.keys(this.state).length > Object.keys(this.params).length) {
1127 : // The state contains devived quantities (not only parameters) and we
1128 : // need to run the log_post once more in order to set the derived quantities
1129 : // for the final parameter state
1130 0 : this.log_post();
1131 0 : }
1132 6001 : return this.state;
1133 1 : }
1134 :
1135 : /**
1136 : * Takes n_iterations steps in the parameter space and returns them
1137 : * as an object of arrays with one array per parameter. For example:
1138 : * {mu: [1, -1, 2, 3, ...], sigma: [1, 2, 2, 1, ...]}.
1139 : * If thin is > 1 then n_iterations / thin samples are returned.
1140 : */
1141 1 : sample(n_iterations: number) {
1142 : // Initializing curr_sample where the sample is going to be saved
1143 : // as an object containing one array per parameter to be monitored.
1144 2 : var i, j, monitored_params;
1145 2 : if (this.monitored_params === null) {
1146 2 : monitored_params = Object.keys(this.state);
1147 0 : } else {
1148 0 : monitored_params = this.monitored_params;
1149 0 : }
1150 :
1151 2 : var curr_sample: any = {};
1152 2 : for (j = 0; j < monitored_params.length; j++) {
1153 4 : curr_sample[monitored_params[j]] = [];
1154 2 : }
1155 :
1156 2 : for (i = 0; i < n_iterations; i++) {
1157 5002 : if (i % this.thinning_interval === 0) {
1158 5002 : for (j = 0; j < monitored_params.length; j++) {
1159 15002 : var param = monitored_params[j];
1160 15002 : curr_sample[param].push(clone_param_draw(this.state[param]));
1161 5002 : }
1162 5002 : }
1163 5002 : this.step();
1164 2 : }
1165 2 : return curr_sample;
1166 1 : }
1167 :
1168 : /**
1169 : * Takes n_iteration steps in parameter space but returns nothing.
1170 : */
1171 1 : burn(n_iterations: number) {
1172 2 : for (var i = 0; i < n_iterations; i++) {
1173 1002 : this.step();
1174 2 : }
1175 1 : }
1176 :
1177 : /**
1178 : * Sets what parameters should be monitored and returned when calling
1179 : * sample.
1180 : */
1181 1 : monitor(params_to_monitor: any) {
1182 2 : this.monitored_params = params_to_monitor;
1183 1 : }
1184 :
1185 : /**
1186 : * Sets the thinning. For example thin == 10 means that every 10th posterior
1187 : * draw will be kept.
1188 : */
1189 1 : thin(thinning_interval: any) {
1190 2 : this.thinning_interval = thinning_interval;
1191 1 : }
1192 :
1193 : /**
1194 : * Sets adaptation on, if applicable, in all steppers.
1195 : */
1196 0 : start_adaptation() {
1197 0 : for (var i = 0; i < this.steppers.length; i++) {
1198 0 : this.steppers[i].start_adaptation();
1199 0 : }
1200 0 : }
1201 :
1202 : /**
1203 : * Sets adaptation off, if applicable, in all steppers.
1204 : */
1205 0 : stop_adaptation() {
1206 0 : for (var i = 0; i < this.steppers.length; i++) {
1207 0 : this.steppers[i].stop_adaptation();
1208 0 : }
1209 0 : }
1210 1 : }
1211 :
1212 : /**
1213 : * @class
1214 : * @implements {Sampler}
1215 : * This sampler uses the AmwgStepper as the stepper function which implements the
1216 : * Adaptive Metropolis-Within-Gibbs algorithm in "Examples of Adaptive MCMC"
1217 : * by Roberts and Rosenthal (2008). An adition is that it handles int parameters
1218 : * by making discrete Normal proposals and binary parameters by taking on a new
1219 : * value proportional to the posterior of the two possible states of the
1220 : * parameter. This sampler can be efficient when the number of parameters
1221 : * are not too high and the correlations between parameters are low.
1222 : * @param params - An object with a one or more parameter definitions
1223 : * @param state - an object containing the state of all parameters.
1224 : * @param log_post - A function that returns the log density that depends on the state.
1225 : * @param options - an object with options to the stepper.
1226 : */
1227 1 : export class AmwgSampler extends Sampler {
1228 1 : constructor(params: any, log_post: Function, data: any, options?: any) {
1229 2 : super(params, log_post, data, options);
1230 1 : }
1231 1 : create_stepper_ensamble(
1232 1 : params: any,
1233 1 : state: any,
1234 1 : log_post: Function,
1235 1 : options: any,
1236 : ) {
1237 2 : return [new AmwgStepper(params, state, log_post, options)];
1238 1 : }
1239 1 : }
|
