LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2025-01-03 01:55:40 +01:00
Code Releases Activity

Optimized Incomplete Gamma

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-09 01:32:07 +01:00
parent 58a5ec1abc
commit 6e2c5130ca
1 changed files with 13 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Utils/Util.ml
View File

@ -20,8 +20,12 @@ let factmax = 150
sha, Tokyo, 1991) p.227 [in Japanese] *) sha, Tokyo, 1991) p.227 [in Japanese] *)
let incomplete_gamma ~alpha x = let incomplete_gamma ~alpha x =
let rec p_gamma a x loggamma_a = let a = alpha in
if x >= 1. +. a then 1. -. q_gamma a x loggamma_a let a_inv = 1./. a in
let gf = gamma_float alpha in
let loggamma_a = log gf in
let rec p_gamma x =
if x >= 1. +. a then 1. -. q_gamma x
else if x = 0. then 0. else if x = 0. then 0.
else else
let rec pg_loop prev res term k = let rec pg_loop prev res term k =
@ -31,21 +35,22 @@ let incomplete_gamma ~alpha x =
let term = term *. x /. (a +. k) in let term = term *. x /. (a +. k) in
pg_loop res (res +. term) term (k +. 1.) pg_loop res (res +. term) term (k +. 1.)
in in
let r0 = exp (a *. log x -. x -. loggamma_a) /. a in let r0 = exp (a *. log x -. x -. loggamma_a) *. a_inv in
pg_loop min_float r0 r0 1. pg_loop min_float r0 r0 1.
and q_gamma a x loggamma_a = and q_gamma x =
if x < 1. +. a then 1. -. p_gamma a x loggamma_a if x < 1. +. a then 1. -. p_gamma x
else else
let rec qg_loop prev res la lb w k = let rec qg_loop prev res la lb w k =
if k > 1000. then failwith "q_gamma did not converge." if k > 1000. then failwith "q_gamma did not converge."
else if prev = res then res else if prev = res then res
else else
let k_inv = 1. /. k in let k_inv = 1. /. k in
let kma = (k -. 1. -. a) *. k_inv in
let la, lb = let la, lb =
lb, ((k -. 1. -. a) *. (lb -. la) +. (k +. x) *. lb) *. k_inv lb, kma *. (lb -. la) +. (k +. x) *. lb *. k_inv
in in
let w = w *. (k -. 1. -. a) *. k_inv in let w = w *. kma in
let prev, res = res, res +. w /. (la *. lb) in let prev, res = res, res +. w /. (la *. lb) in
qg_loop prev res la lb w (k +. 1.) qg_loop prev res la lb w (k +. 1.)
in in
@ -53,8 +58,7 @@ let incomplete_gamma ~alpha x =
let lb = (1. +. x -. a) in let lb = (1. +. x -. a) in
qg_loop min_float (w /. lb) 1. lb w 2.0 qg_loop min_float (w /. lb) 1. lb w 2.0
in in
let gf = gamma_float alpha in gf *. p_gamma x
gf *. p_gamma alpha x (log gf)