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
# 1 "src/owl/stats/owl_stats_sampler.ml"
(*
 * OWL - OCaml Scientific and Engineering Computing
 * Copyright (c) 2016-2020 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