Commit 6815030e authored by Philippe Veber's avatar Philippe Veber
Browse files

tk/Multinomial_test: added a permutation test

parent 5b2db384
......@@ -155,6 +155,38 @@ module Sparse = struct
let simulation_test = simulation_test ~likelihood_log_ratio
end
module Permutation = struct
let expand t =
Array.concat_mapi t ~f:(fun i j ->
Array.create ~len:j i
)
let counts t ~k ~pos ~len =
let res = Array.create 0 ~len:k in
assert (pos + len <= Array.length t) ;
for i = pos to pos + len - 1 do
let j = t.(i) in
res.(j) <- res.(j) + 1
done ;
res
let resample d =
let n1 = isum d.x1 in
let n2 = isum d.x2 in
let t = Array.append (expand d.x1) (expand d.x2) in
fun () ->
Array.permute t ; (* FIXME: use gsl rng *)
{ d with x1 = counts t ~k:d.k ~pos:0 ~len:n1 ;
x2 = counts t ~k:d.k ~pos:n1 ~len:n2 }
let test ?(_B_ = 1_000) d ~statistic =
let _T_ = statistic d in
let resample = resample d in
let n = Utils.range_count 1 _B_ ~f:Float.(fun _ -> statistic (resample ()) > _T_) in
let pvalue = float n /. float _B_ in
{ _T_ ; pvalue }
end
let random_discrete_probability k ~alpha =
let x = Array.init k ~f:(fun _ -> Owl.Stats.gamma_rvs ~shape:alpha ~scale:1.) in
let s = Array.reduce_exn x ~f:( +. ) in
......
......@@ -33,6 +33,13 @@ module Sparse : sig
val simulation_test : ?sample_size:int -> data -> result
end
module Permutation : sig
val test :
?_B_:int ->
data -> statistic:(data -> float) ->
result
end
val uniformity_test :
?alpha:float ->
?filter_out_single_category:bool ->
......
......@@ -47,6 +47,13 @@ let range_find a b ~f =
in
loop a
let range_count a b ~f =
let rec loop acc i =
if i < b then loop (if f i then acc + 1 else acc) (i + 1)
else acc
in
loop 0 a
let int_fold a b ~init ~f =
let rec loop acc i = if i >= b then acc else loop (f acc i) (i + 1) in
loop init a
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment