package owl

  1. Overview
  2. Docs
package owl
Legend:
Page
Library
Module
Module type
Parameter
Class
Class type
Source

Source file owl_stats_sampler.ml

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# 1 "src/owl/stats/owl_stats_sampler.ml"
(*
 * OWL - OCaml Scientific and Engineering Computing
 * Copyright (c) 2016-2019 Liang Wang <liang.wang@cl.cam.ac.uk>
 *)

type 'a t = {
  samples : 'a array;
  accept  : float;
}


let rejection ~m ~proprvs ~proppdf ~pdf nsamples =
  assert (m > 0.);
  let log_m = log m in
  let accept = ref 0 in
  let total = ref 0 in
  let samples = Array.make nsamples 0. in

  while !accept < nsamples do
    let x = proprvs () in
    let a = log (Owl_stats_dist.std_uniform_rvs ()) in
    let b = (log (pdf x)) -. (log (proppdf x)) -. log_m in
    assert (b < 0.);
    if a < b then (
      samples.(!accept) <- x;
      accept := !accept + 1;
    );
    total := !total + 1
  done;

  let accept = (float_of_int !accept) /. (float_of_int !total) in
  { samples; accept }


let ars = None


let arms = None


let metropolis ?burnin ?thin ~initial ~proprvs ~proppdf ~pdf nsamples =
  let burnin = match burnin with Some a -> a | None -> 1000 in
  let thin = match thin with Some a -> a | None -> 10 in

  let accept = ref 0 in
  let niter = burnin + thin * nsamples in
  let samples = Array.make niter [||] in
  samples.(0) <- Array.copy initial;

  for i = 1 to niter - 1 do
    let x  = samples.(i - 1) in
    let x' = proprvs x in
    let p  = log (pdf x) in
    let p' = log (pdf x') in
    let g  = log (proppdf x x') in
    let g' = log (proppdf x' x) in
    let a  = min 0. (p' -. p +. g -. g') in
    let b =
      if a >= 0. then true
      else if log (Owl_stats_dist.std_uniform_rvs ()) < a then true
      else false
    in
    let x' = if b then x' else Array.copy x in
    if b then accept := !accept + 1;
    samples.(i) <- x'
  done;

  let samples = Owl_utils.Array.filteri (fun i _ ->
    (i >= burnin) && (i mod thin = 0)
  ) samples
  in
  let accept = (float_of_int !accept) /. (float_of_int niter) in
  { samples; accept }


let gibbs = None


let slice = None


let adaptive_rejection = None


let to_string t =
  Printf.sprintf "[ samples: %i; accept_ratio: %g ]"
  (Array.length t.samples) t.accept


let pp_t formatter t =
  Format.open_box 0;
  Format.fprintf formatter "%s" (to_string t);
  Format.close_box ()
OCaml

Innovation. Community. Security.

GitHub Discord Twitter Peertube RSS
About
Changelog Releases Industrial Users Academic Users Why OCaml
Ecosystem
Packages Community Events OCaml Planet Jobs
Resources
Install OCaml Get Started Platform Tools Language Manual Standard Library API Books Exercises Papers OCaml Playground Logo
Policies
Carbon Footprint Governance Privacy Code of Conduct