Commit 943edcc8 authored by Philippe Veber's avatar Philippe Veber
Browse files

prc: revamped + added interface

parent c0ca7c3c
open Core_kernel
open Gsl
type classification_data = Classification_data of (float * bool) list
type evaluation_data =
Evaluation_data of (float * bool) list
let maybe_min x = Array.fold x ~init:None ~f:(fun acc x_i ->
match acc, x_i with
| Some u, Some v -> Some (Float.min u v)
| None, Some _ -> x_i
| Some _, None -> acc
| None, None -> None
)
let evaluation_data x y =
let open Result.Monad_infix in
(
if Array.length x = Array.length y then Ok ()
else Error (`Msg "expected arrays of equal lengths")
) >>= fun () ->
Result.of_option (maybe_min x) ~error:(`Msg "empty score table") >>| fun min_x ->
Evaluation_data (
List.init (Array.length x) ~f:(fun i ->
match x.(i), y.(i) with
| Some x_i, y_i -> (x_i, y_i)
| None, y_i -> (min_x, y_i)
)
)
type curve = {
x : float array ;
y : float array ;
col : string option ;
lwd : int option ;
label : string option ;
}
let empirical_curve (Classification_data xs) =
let xy xs =
let curve = Array.of_list xs in
let x = Array.map curve ~f:fst in
let y = Array.map curve ~f:snd in
x, y
let empirical_curve_points (Evaluation_data xs) =
let npos = List.count xs ~f:snd in
List.sort xs ~compare:(Tuple.T2.compare ~cmp1:Float.descending ~cmp2:Bool.descending)
|> List.fold ~init:(0, 0, [ 0., 1. ]) ~f:(fun (pos, neg, points) (_, b) ->
let pos, neg = if b then pos + 1, neg else pos, neg +1 in
let recall = float pos /. float npos in
let precision = float pos /. float (pos + neg) in
(pos, neg, (recall, precision) :: points)
)
|> trd3
|> List.rev
let compare (x, _) (y, _) = Float.compare y x in
let recall pos = float pos /. float npos in
let precision pos neg = float pos /. float (pos + neg) in
let data = List.sort xs ~compare in
let rec loop pos neg previous_threshold = function
| [] -> [ recall pos, precision pos neg ]
| (x, b) :: t ->
let pos, neg = if b then pos + 1, neg else pos, neg + 1 in
let tail = loop pos neg (Some x) t in
let group_end =
Option.value_map previous_threshold ~default:false ~f:Float.(( <> ) x)
in
if group_end then (recall pos, precision pos neg) :: tail
else tail
in
loop 0 0 None data
let empirical_curve ?col ?label ?lwd data =
let x, y = xy (empirical_curve_points data) in
{ x ; y ; col ; label ; lwd }
(** Binormal model *)
(**
The binormal assumption on precision-recall curves
Kay H. Brodersen ∗† , Cheng Soon Ong ∗ , Klaas E. Stephan † and Joachim M. Buhmann ∗
*)
type binormal_params = {
mu_pos : float ;
sigma_pos : float ;
......@@ -30,7 +77,6 @@ type binormal_params = {
let binormal_params ?(mu_pos = 1.) ?(sigma_pos = 1.) ?(mu_neg = 0.) ?(sigma_neg = 1.) alpha =
{ mu_pos ; mu_neg ; sigma_pos ; sigma_neg ; alpha }
(** Binormal simulator: Z ~ alpha N(mu+,sigma+) + (1 - alpha) N(mu-, sigma-) *)
let binormal_simulation rng ~n { mu_pos ; mu_neg ; sigma_pos ; sigma_neg ; alpha } =
let sample = List.init n ~f:(fun _ ->
let mu, sigma, label = match Randist.bernoulli rng ~p:alpha with
......@@ -41,30 +87,32 @@ let binormal_simulation rng ~n { mu_pos ; mu_neg ; sigma_pos ; sigma_neg ; alpha
mu +. Randist.gaussian rng ~sigma, label
)
in
Classification_data sample
Evaluation_data sample
let phi ~mu ~sigma x = Cdf.gaussian_P ~x:(x -. mu) ~sigma
let inv_phi ~mu ~sigma p = Cdf.gaussian_Pinv ~p ~sigma +. mu
let binormal_curve ~n p =
let binormal_precision p recall =
let phi_neg = phi ~mu:p.mu_neg ~sigma:p.sigma_neg in
let inv_phi_pos = inv_phi ~mu:p.mu_pos ~sigma:p.sigma_pos in
let alpha = p.alpha in
let xs = List.init n ~f:(fun i -> float i /. float n) in
List.map xs ~f:(fun x_i ->
let y_i =
alpha *. x_i
/.
(alpha *. x_i
+.
(1. -. alpha) *. (1. -. phi_neg (inv_phi_pos (1. -. x_i))))
in
x_i, y_i
)
let binormal_estimate (Classification_data xs) =
alpha *. recall
/.
(alpha *. recall
+.
(1. -. alpha) *. (1. -. phi_neg (inv_phi_pos (1. -. recall))))
let binormal_curve ?(n = 100) ?col ?lwd ?label p =
let max_i = float (n - 1) in
let x = Array.init n ~f:(fun i -> float i /. max_i) in
let y = Array.map x ~f:(binormal_precision p) in
{ x ; y ; col ; lwd ; label }
let binormal_estimate (Evaluation_data xs) =
let select label =
List.filter_map xs ~f:(fun (x, b) -> if Bool.(b = label) then Some x else None)
List.filter_map xs ~f:(fun (x, b) ->
if Bool.(b = label) then Some x else None
)
|> Array.of_list
in
let x = select false in
......@@ -78,45 +126,11 @@ let binormal_estimate (Classification_data xs) =
alpha = len y /. (len x +. len y) ;
}
let xy xs =
let curve = Array.of_list xs in
let x = Array.map curve ~f:fst in
let y = Array.map curve ~f:snd in
x, y
let plot_binormal_fit ?(seed = 42) ?(sigma = 1.) ~n ~alpha () =
let p = binormal_params ~sigma_pos:sigma ~sigma_neg:sigma alpha in
let plot ?lwd ?col points =
let x, y = xy points in
OCamlR_graphics.lines ?lwd ?col ~x ~y ()
in
let rng = Rng.(make (default ())) in
Rng.set rng (Nativeint.of_int seed) ;
OCamlR_graphics.plot
~xlab:"Recall" ~ylab:"Precision"
~xlim:(0., 1.) ~ylim:(0., 1.) ~x:[| |] ~y:[| |] () ;
for _ = 1 to 10 do
binormal_simulation ~n rng p
|> empirical_curve
|> plot ~col:"lightgrey"
done ;
let sim = binormal_simulation rng ~n p in
let estimation = binormal_estimate sim in
plot ~col:"red" (empirical_curve sim) ;
plot ~lwd:2 ~col:"black" (binormal_curve ~n:100 p) ;
plot ~lwd:2 ~col:"green" (binormal_curve ~n:100 estimation) ;
OCamlR_graphics.legend
~col:[| "lightgrey" ; "black" |]
~lwd:[| 1. ; 2. |]
`topright
[|"Simulations" ; "Theoretical curve"|]
(**
Area Under the Precision-Recall Curve:
Point Estimates and Confidence Intervals
Kendrick Boyd 1 , Kevin H. Eng 2 , and C. David Page
*)
let auc_binormal p =
let n = 1_000 in
let x = Array.init n ~f:(fun i -> float (i + 1) /. float n) in
let y = Array.map x ~f:(binormal_precision p) in
Stats.mean y
let sum a b ~f =
let rec loop i acc =
......@@ -126,7 +140,7 @@ let sum a b ~f =
loop a 0.
let auc_trapezoidal_lt classification =
let curve = empirical_curve classification in
let curve = empirical_curve_points classification in
let n, r, pmin, pmax =
let data =
List.group curve ~break:(fun (r1,_) (r2, _) -> Float.(r1 <> r2))
......@@ -143,7 +157,7 @@ let auc_trapezoidal_lt classification =
in
sum 0 (n - 2) ~f:(fun i -> (pmin.(i) +. pmax.(i + 1)) /. 2. *. (r.(i + 1) -. r.(i)))
let auc_average_precision (Classification_data xs) =
let auc_average_precision (Evaluation_data xs) =
let pos, _, sum =
List.sort xs ~compare:(Tuple.T2.compare ~cmp1:Float.descending ~cmp2:Bool.descending)
|> List.fold ~init:(0, 0, 0.) ~f:(fun (pos, neg, sum) (_, b) ->
......@@ -159,32 +173,13 @@ let auc_average_precision (Classification_data xs) =
in
sum /. float pos
let plot_auc_boxplots ?(seed = 42) ?(sigma = 1.) ~n ~alpha () =
let open OCamlR_base in
let p = binormal_params ~sigma_pos:sigma ~sigma_neg:sigma alpha in
let rng = Rng.(make (default ())) in
Rng.set rng (Nativeint.of_int seed) ;
let samples = List.init 1_000 ~f:(fun _ -> binormal_simulation ~n rng p) in
let compute f =
List.map samples ~f
|> Numeric.of_list
|> Numeric.to_sexp
in
let trapezoidal_lt = compute auc_trapezoidal_lt in
let average_precision = compute auc_average_precision in
let l = List_.create [
Some "trapezoidal", trapezoidal_lt ;
Some "average_precision", average_precision ;
]
in
OCamlR_graphics.list_boxplot l
let logit p = Float.log (p /. (1. -. p))
let sigmoid x =
let exp_x = Float.exp x in
exp_x /. (1. +. exp_x)
let logit_confidence_interval ~alpha ~theta_hat ~n =
let eta_hat = logit theta_hat in
let tau_hat = (float n *. theta_hat *. (1. -. theta_hat)) ** (-0.5) in
......@@ -192,19 +187,12 @@ let logit_confidence_interval ~alpha ~theta_hat ~n =
sigmoid (eta_hat -. delta),
sigmoid (eta_hat +. delta)
let logit_confidence_interval_coverage_test ?(seed = 42) ?(sigma = 1.) ~n ~alpha () =
let p = binormal_params ~sigma_pos:sigma ~sigma_neg:sigma alpha in
let rng = Rng.(make (default ())) in
Rng.set rng (Nativeint.of_int seed) ;
let nsims = 1_000 in
let samples = List.init nsims ~f:(fun _ -> binormal_simulation ~n rng p) in
let npos = List.map samples ~f:(fun (Classification_data xs) -> List.count xs ~f:snd) in
let estimates = List.map samples ~f:auc_trapezoidal_lt in
let confidence_intervals = List.map2_exn estimates npos ~f:(fun theta_hat n ->
logit_confidence_interval ~alpha:0.05 ~theta_hat ~n
)
in
let best_estimate = Array.of_list estimates |> Stats.mean in
let nsuccessess = List.count confidence_intervals ~f:(fun (lb, ub) -> Float.(lb <= best_estimate && best_estimate <= ub)) in
float nsuccessess /. float nsims
let plot_curve { x ; y ; col ; lwd ; _ } =
OCamlR_graphics.lines ?lwd ?col ~x ~y ()
let recall_precision_plot ?main curves =
OCamlR_graphics.plot
~plot_type:`Lines ?main
~xlab:"Recall" ~ylab:"Precision"
~xlim:(0., 1.) ~ylim:(0., 1.) ~x:[||] ~y:[||] () ;
List.(iter (rev curves)) ~f:plot_curve
(**
References:
- The binormal assumption on precision-recall curves,
Kay H. Brodersen, Cheng Soon Ong, Klaas E. Stephan and Joachim M. Buhmann
- Area Under the Precision-Recall Curve: Point Estimates and Confidence Intervals,
Kendrick Boyd, Kevin H. Eng and C. David Page
*)
type evaluation_data = Evaluation_data of (float * bool) list
val evaluation_data :
float option array ->
bool array ->
(evaluation_data, [> `Msg of string]) result
type curve = {
x : float array ;
y : float array ;
col : string option ;
lwd : int option ;
label : string option ;
}
val empirical_curve :
?col:string ->
?label:string ->
?lwd:int ->
evaluation_data ->
curve
(** Binormal model *)
type binormal_params = {
mu_pos : float ;
sigma_pos : float ;
mu_neg : float ;
sigma_neg : float ;
alpha : float ;
}
val binormal_params :
?mu_pos:float ->
?sigma_pos:float ->
?mu_neg:float ->
?sigma_neg:float ->
float ->
binormal_params
val binormal_simulation :
Gsl.Rng.t ->
n:int ->
binormal_params ->
evaluation_data
val binormal_curve :
?n:int ->
?col:string ->
?lwd:int ->
?label:string ->
binormal_params ->
curve
val binormal_estimate : evaluation_data -> binormal_params
val auc_binormal : binormal_params -> float
val auc_trapezoidal_lt : evaluation_data -> float
val auc_average_precision : evaluation_data -> float
val logit_confidence_interval :
alpha:float ->
theta_hat:float ->
n:int ->
float * float
(** [logit_confidence_interval ~alpha ~theta_hat ~n] computes an
asymptotically valid confidence interval at level 1 - [alpha], when
the estimate [theta_hat] was obtained from a sample of [n]
observations. *)
val recall_precision_plot :
?main:string ->
curve list ->
unit
......@@ -190,11 +190,15 @@ module Mutsel = struct
let ancestral_counts = collect_counts `Ancestral in
let convergent_counts = collect_counts `Convergent in
let make_classification_data x y =
Prc.Classification_data (
let min_x = Array.fold x ~init:Float.neg_infinity ~f:(fun acc x_i ->
Option.value_map ~default:acc x_i ~f:(Float.min acc)
) in
Prc.Evaluation_data (
List.init (Array.length x) ~f:(fun i ->
match x.(i), y.(i) with
| Some x_i, Some y_i -> Some (x_i, y_i)
| None, _ | _, None -> None
| None, Some y_i -> Some (min_x, y_i)
| _, None -> None
)
|> List.filter_opt
)
......@@ -202,7 +206,7 @@ module Mutsel = struct
let average_precision =
let oracle = Array.init nsites ~f:(fun i -> if i < n_h0 then Some false else Some true) in
List.map results ~f:(fun (_, scores) ->
let Prc.Classification_data xs as data = make_classification_data scores oracle in
let Prc.Evaluation_data xs as data = make_classification_data scores oracle in
let n = List.count xs ~f:snd in
let theta_hat = Prc.auc_trapezoidal_lt data in
let lb, ub = Prc.logit_confidence_interval ~alpha:0.05 ~theta_hat ~n in
......@@ -215,7 +219,7 @@ module Mutsel = struct
average_precision ; profiles ; n_h0 ; n_ha
}
in
Workflow.plugin ~descr:"simulation_pipeline.benchmark" f
Workflow.plugin ~version:3 ~descr:"simulation_pipeline.benchmark" f
let rds_of_benchmark b =
let f = fun%workflow dest ->
......
......@@ -104,18 +104,22 @@ let average_precision_plot ~labels ~oracle ~results : pdf file =
let result_paths = [%eval Bistro.Workflow.path_list (List.map ~f:fst results)] in
let results = List.map2_exn result_paths (List.map ~f:snd results) ~f:load_results in
let make_classification_data x y =
Prc.Classification_data (
let min_x = Array.fold x ~init:Float.neg_infinity ~f:(fun acc x_i ->
Option.value_map ~default:acc x_i ~f:(Float.min acc)
) in
Prc.Evaluation_data (
List.init (Array.length x) ~f:(fun i ->
match x.(i), y.(i) with
| Some x_i, Some y_i -> Some (x_i, y_i)
| None, _ | _, None -> None
| None, Some y_i -> Some (min_x, y_i)
| _, None -> None
)
|> List.filter_opt
)
in
let estimates, lower_bounds, upper_bounds =
List.map results ~f:(fun scores ->
let Prc.Classification_data xs as data = make_classification_data scores oracle in
let Prc.Evaluation_data xs as data = make_classification_data scores oracle in
let n = List.count xs ~f:snd in
let theta_hat = Prc.auc_trapezoidal_lt data in
let lb, ub = Prc.logit_confidence_interval ~alpha:0.05 ~theta_hat ~n in
......
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