superbobry / pareto

GSL powered OCaml statistics library
http://superbobry.github.io/pareto/0.2
MIT License
40 stars 5 forks source link

BCa bootstrap is wrong #39

Open edwintorok opened 4 years ago

edwintorok commented 4 years ago
let z_2s = Array.map jack ~f:(fun v -> sqr (jack_mean -. v)) in
let sum_cubes   = Array.fold_left z_2s ~init:0.
          ~f:(fun acc z_2 -> acc +. sqr z_2)
and sum_squares = Array.fold_left z_2s ~f:(+.) ~init:0. in

Instead of cubes (x^3) this is calculating x^4 here. This causes it to give obviously wrong results, like in the following test:

open Pareto.Resampling.Bootstrap
let () =
  let s1 = [|1.0;2.0;3.0;4.0;5.0; 1000.0;2000.0|] in
  let b = bca ~estimator:Pareto.Sample.mean ~n_iter:1000 s1 in
  Printf.printf "%f, %f, %f\n" b.lower_bound b.point b.upper_bound

It prints the same value for low/upper-bound and a point that is outside of it: 430.571429, 430.714286, 430.571429

Fixing the sum_cubes to actually sum cubes I get something more reasonably looking: 3.142857, 430.714286, 1286.428571

superbobry commented 4 years ago

Thanks for the bug report! I'm not actively looking after pareto these days and will happily accept a PR fixing that.