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

prc: reorganized API

parent 9806a257
(library
(name prc)
(public_name codepi.prc)
(libraries core_kernel gsl ocaml-r.graphics))
(libraries core_kernel gsl ocaml-r.graphics)
(inline_tests)
(preprocess (pps ppx_jane)))
open Core_kernel
open Gsl
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 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
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
let sum a b ~f =
let rec loop i acc =
if i > b then acc
else loop (i + 1) (acc +. f i)
in
loop 0 0 None data
loop a 0.
let empirical_curve ?col ?label ?lwd data =
let x, y = xy (empirical_curve_points data) in
{ x ; y ; col ; label ; lwd }
let logit p = Float.log (p /. (1. -. p))
(** Binormal model *)
let sigmoid x =
let exp_x = Float.exp x in
exp_x /. (1. +. exp_x)
type dataset = Dataset of (float * bool) list
type binormal_params = {
mu_pos : float ;
sigma_pos : float ;
mu_neg : float ;
sigma_neg : float ;
alpha : float ;
let confusion_matrix_fold (Dataset d) ~init ~f =
let xs = Array.of_list d in
Array.sort xs ~compare:(fun (x, _) (y, _) -> Float.compare y x) ;
let n = Array.length xs in
let npos = Array.count xs ~f:snd in
let nneg = n - npos in
let rec loop i acc ~tp ~fp ~fn ~tn =
if i <= n then
let score, b = xs.(i - 1) in
let tp = if b then tp + 1 else tp
and fp = if b then fp else fp + 1
and fn = if b then fn - 1 else fn
and tn = if b then tn else tn - 1 in
let acc =
if i = n || Float.(score <> fst xs.(i)) then
f ~tp ~fp ~fn ~tn ~score acc
else acc
in
loop (i + 1) acc ~tp ~tn ~fp ~fn
else acc
in
let tp = 0 and tn = nneg and fp = 0 and fn = npos in
loop 1 (f ~tp ~fp ~fn ~tn ~score:Float.infinity init) ~tp ~fp ~tn ~fn
type confusion_matrix = {
tp : int ;
tn : int ;
fp : int ;
fn : int ;
}
[@@deriving sexp]
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 }
type confusion_matrix_series = (float * confusion_matrix) array
[@@deriving sexp]
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
| 0 -> mu_neg, sigma_neg, false
| 1 -> mu_pos, sigma_pos, true
| _ -> assert false
in
mu +. Randist.gaussian rng ~sigma, label
let confusion_matrix_series d =
confusion_matrix_fold d ~init:[] ~f:(fun ~tp ~fp ~fn ~tn ~score acc ->
(score, { tp ; tn ; fp ; fn }) :: acc
)
in
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_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
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.rev
|> Array.of_list
let%expect_test "performance curve 1" =
let scores = [2.1 ; 1.2 ; 5.6 ; 0. ; 5.6] in
let labels = [true ; false ; true ; false ; false] in
let dataset = Dataset (List.zip_exn scores labels) in
let series = confusion_matrix_series dataset in
print_endline (Sexp.to_string_hum (sexp_of_confusion_matrix_series series)) ;
[%expect "
((INF ((tp 0) (tn 3) (fp 0) (fn 2))) (5.6 ((tp 1) (tn 2) (fp 1) (fn 1)))
(2.1 ((tp 2) (tn 2) (fp 1) (fn 0))) (1.2 ((tp 2) (tn 1) (fp 2) (fn 0)))
(0 ((tp 2) (tn 0) (fp 3) (fn 0))))"]
let recall ~tp ~fn = float tp /. float (tp + fn)
let precision ~tp ~fp = float tp /. float (tp + fp)
module Precision_recall = struct
let operating_points d =
confusion_matrix_fold d ~init:[] ~f:(fun ~tp ~fp ~fn ~tn:_ ~score acc ->
(score, recall ~tp ~fn, precision ~tp ~fp) :: acc
)
|> Array.of_list
in
let x = select false in
let y = select true in
let len t = float (Array.length t) in
Stats.{
mu_pos = mean y ;
mu_neg = mean x ;
sigma_pos = sd y ;
sigma_neg = sd x ;
alpha = len y /. (len x +. len y) ;
}
|> List.rev
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 auc_trapezoidal_lt d =
let points = operating_points d in
let n, r, pmin, pmax =
let data =
List.group points ~break:(fun (_, r1, _) (_, r2, _) ->
Float.(r1 <> r2)
)
|> List.map ~f:(function
| [] -> assert false
| (_, r, _) :: _ as xs -> r, List.map xs ~f:trd3
)
|> Array.of_list
in
Array.length data,
Array.map data ~f:fst,
Array.map data ~f:(fun (_, ps) -> List.reduce_exn ps ~f:Float.min),
Array.map data ~f:(fun (_, ps) -> List.reduce_exn ps ~f:Float.max)
in
sum 0 (n - 2) ~f:(fun i -> (pmin.(i) +. pmax.(i + 1)) /. 2. *. (r.(i + 1) -. r.(i)))
let sum a b ~f =
let rec loop i acc =
if i > b then acc
else loop (i + 1) (acc +. f i)
in
loop a 0.
let decreasing_ties_groups (Dataset xs) =
let rec loop closed_groups current_group points =
match points, current_group with
| [], None -> []
| [], Some (s, npos, nneg) ->
(s, npos, nneg) :: closed_groups
| (s, b) :: t, None ->
let npos, nneg = if b then 1, 0 else 0, 1 in
loop closed_groups (Some (s, npos, nneg)) t
| (s', b) :: t, Some (s, npos, nneg) ->
if Float.equal s s' then
let npos, nneg = if b then npos + 1, nneg else npos, nneg + 1 in
loop closed_groups (Some (s, npos, nneg)) t
else
let npos', nneg' = if b then 1, 0 else 0, 1 in
loop ((s, npos, nneg) :: closed_groups) (Some (s', npos', nneg')) t
in
let compare = Tuple.T2.compare ~cmp1:Float.descending ~cmp2:Bool.descending in
List.sort xs ~compare
|> loop [] None
let auc_trapezoidal_lt classification =
let curve = empirical_curve_points classification in
let n, r, pmin, pmax =
let data =
List.group curve ~break:(fun (r1,_) (r2, _) -> Float.(r1 <> r2))
|> List.map ~f:(function
| [] -> assert false
| (r, _) :: _ as xs -> r, List.map xs ~f:snd
let%expect_test "decreasing_ties_groups" =
let data = Dataset [0., false; 1., false; 2., true; 0., false; 1., false; 2., true; 0., true; 1., false; 2., true; 0., true] in
let groups = decreasing_ties_groups data in
Sexplib.Sexp.pp Format.std_formatter ([%sexp_of: (float * int * int) list] groups) ;
[%expect "((0 2 2)(1 0 3)(2 3 0))"]
let auc_average_precision (Dataset 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) ->
let pos, neg = if b then pos + 1, neg else pos, neg +1 in
let sum =
if b then
let precision = float pos /. float (pos + neg) in
precision +. sum
else sum
in
(pos, neg, sum)
)
|> Array.of_list
in
Array.length data,
Array.map data ~f:fst,
Array.map data ~f:(fun (_, ps) -> List.reduce_exn ps ~f:Float.min),
Array.map data ~f:(fun (_, ps) -> List.reduce_exn ps ~f:Float.max)
in
sum 0 (n - 2) ~f:(fun i -> (pmin.(i) +. pmax.(i + 1)) /. 2. *. (r.(i + 1) -. r.(i)))
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) ->
let pos, neg = if b then pos + 1, neg else pos, neg +1 in
let sum =
if b then
let precision = float pos /. float (pos + neg) in
precision +. sum
else sum
sum /. float pos
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
let delta = tau_hat *. Cdf.gaussian_Pinv ~sigma:1. ~p:(1. -. alpha /. 2.) in
sigmoid (eta_hat -. delta),
sigmoid (eta_hat +. delta)
end
module Binormal_model = struct
type t = {
mu_pos : float ;
sigma_pos : float ;
mu_neg : float ;
sigma_neg : float ;
alpha : float ;
}
let make ?(mu_pos = 1.) ?(sigma_pos = 1.) ?(mu_neg = 0.) ?(sigma_neg = 1.) alpha =
{ mu_pos ; mu_neg ; sigma_pos ; sigma_neg ; alpha }
let 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
| 0 -> mu_neg, sigma_neg, false
| 1 -> mu_pos, sigma_pos, true
| _ -> assert false
in
(pos, neg, sum)
mu +. Randist.gaussian rng ~sigma, label
)
in
sum /. float pos
in
Dataset sample
let logit p = Float.log (p /. (1. -. p))
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 sigmoid x =
let exp_x = Float.exp x in
exp_x /. (1. +. exp_x)
let 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
alpha *. recall
/.
(alpha *. recall
+.
(1. -. alpha) *. (1. -. phi_neg (inv_phi_pos (1. -. recall))))
let curve ?(n = 100) p =
let max_i = float (n - 1) in
Array.init n ~f:(fun i ->
let r = float i /. max_i in
r, precision p r
)
let estimate (Dataset xs) =
let select label =
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
let y = select true in
let len t = float (Array.length t) in
Stats.{
mu_pos = mean y ;
mu_neg = mean x ;
sigma_pos = sd y ;
sigma_neg = sd x ;
alpha = len y /. (len x +. len y) ;
}
let auc 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:(precision p) in
Stats.mean y
end
module Plot = struct
type t = {
x : float array ;
y : float array ;
col : string option ;
ty : [`Lines of int | `Points of int] ;
label : string option ;
}
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
let delta = tau_hat *. Cdf.gaussian_Pinv ~sigma:1. ~p:(1. -. alpha /. 2.) in
sigmoid (eta_hat -. delta),
sigmoid (eta_hat +. delta)
let lines ?col ?(lwd = 1) ?label points = {
x = Array.map points ~f:fst ;
y = Array.map points ~f:snd ;
col ; ty = `Lines lwd ; label ;
}
let points ?col ?(pch = 19) ?label points = {
x = Array.map points ~f:fst ;
y = Array.map points ~f:snd ;
col ; ty = `Points pch ; label ;
}
let plot_curve { x ; y ; col ; lwd ; _ } =
OCamlR_graphics.lines ?lwd ?col ~x ~y ()
let plot_curve { x ; y ; col ; ty ; _ } =
match ty with
| `Lines lwd ->
OCamlR_graphics.lines ~lwd ?col ~x ~y ()
| `Points pch ->
OCamlR_graphics.points ~pch ?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
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
end
(**
References:
- The binormal assumption on precision-recall curves.
Kay H. Brodersen, Cheng Soon Ong, Klaas E. Stephan and Joachim M. Buhmann
[1] 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
[2] Area Under the Precision-Recall Curve: Point Estimates and Confidence Intervals.
Kendrick Boyd, Kevin H. Eng and C. David Page
- Precision-Recall-Gain Curves: PR Analysis Done Right.
Peter A. Flach and Meelis Kull
[3] Precision-Recall-Gain Curves: PR Analysis Done Right.
Peter A. Flach and Meelis Kull
- The Relationship Between Precision-Recall and ROC Curves.
Jesse Davis and Mark Goadrich
- Realisable Classifiers: Improving Operating Performance on Variable Cost Problems.
M.J.J. Scott, M. Niranjan, R.W. Prager
[4] The Relationship Between Precision-Recall and ROC Curves.
Jesse Davis and Mark Goadrich
[5] Realisable Classifiers: Improving Operating Performance on Variable Cost Problems.
M.J.J. Scott, M. Niranjan, R.W. Prager
*)
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
type dataset = Dataset of (float * bool) list
(** Binary prediction scores with associated labels *)
module Precision_recall : sig
val operating_points :
dataset ->
(float * float * float) list
(** [operating_points d] computes the list of score threshold, recall
and precision triplets, sorted by decreasing threshold. *)
val auc_trapezoidal_lt : dataset -> float
(** AUC lower triangular estimator (see [2] for reference) *)
val auc_average_precision : dataset -> float
(** AUC average precision (see [2] for reference) *)
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. *)
end
(** Binormal model
A Gaussian mixture model for which the precision-recall curve can
be computed explicitly (see [1])
*)
module Binormal_model : sig
type t = {
mu_pos : float ;
sigma_pos : float ;
mu_neg : float ;
sigma_neg : float ;
alpha : float ;
}
val make :
?mu_pos:float ->
?sigma_pos:float ->
?mu_neg:float ->
?sigma_neg:float ->
float ->
t
val simulation :
Gsl.Rng.t ->
n:int ->
t ->
dataset
val curve :
?n:int ->