Commit d8691ea5 authored by Philippe Veber's avatar Philippe Veber
Browse files

tk/Multinomial_test: added chi-squared test implementation

parent efa16120
......@@ -60,6 +60,42 @@ let simulation_test ?(sample_size = 10_000) ~likelihood_log_ratio d =
let pvalue = float !c /. float sample_size in
{ _T_ = t_obs ; pvalue }
module Chisq = struct
let statistic d =
let n1 = isum d.x1
and n2 = isum d.x2 in
let n = float (n1 + n2) in
let p_hat = Array.init d.k ~f:(fun j -> float (d.x1.(j) + d.x2.(j)) /. n) in
sum 0 d.k (fun j ->
let e_hat_1j = p_hat.(j) *. float n1
and e_hat_2j = p_hat.(j) *. float n2 in
((float d.x1.(j) -. e_hat_1j) ** 2.) /. e_hat_1j +.
((float d.x2.(j) -. e_hat_2j) ** 2.) /. e_hat_2j
)
let test d =
let _T_ = statistic d in
let nu = float (d.k - 1) in
let pvalue = Gsl.Cdf.chisq_Q ~x:_T_ ~nu in
{ _T_ ; pvalue }
let%test _ =
let x1 = [|3; 0; 5|] and x2 = [|2; 3; 10|] in
let d = data ~x1 ~x2 in
let r = test d in
Float.(abs (r._T_ - 3.0156) < 0.0001) &&
Float.(abs (r.pvalue - 0.2214) < 0.0001)
(*
> chisq.test(x)
Pearson's Chi-squared test
data: x
X-squared = 3.0156, df = 2, p-value = 0.2214
*)
end
module LRT = struct
let likelihood_log_ratio ({ x1 ; x2 ; _ } as d) =
let logpdf = Owl.Stats.multinomial_logpdf in
......
(** Two-sample multinomial tests *)
(** Two-sample multinomial tests
References:
- Exact inference for categorical data: recent advances and continuing controversies, Alan Agresti
*)
type data
......@@ -12,6 +16,11 @@ type result = {
pvalue : float ;
}
module Chisq : sig
val statistic : data -> float
val test : data -> result
end
module LRT : sig
val likelihood_log_ratio : data -> float
val asymptotic_test : ?df_correction:bool -> data -> result
......
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