* Fwd: log probability density functions (PDFs)
[not found] <CAPEJTf6T8vk+KgnoDBaZgJzjRUvJfLK+w-S-xwoZ_=f_wgNH3g@mail.gmail.com>
@ 2012-07-09 17:30 ` Sam Mason
2012-07-09 17:39 ` Patrick Alken
[not found] ` <d75321b0fc2d4b7fb010d63b6252dd89@IDMEXCHT02.ads.warwick.ac.uk>
0 siblings, 2 replies; 3+ messages in thread
From: Sam Mason @ 2012-07-09 17:30 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 1201 bytes --]
Hi,
I sent the following message to the help-gsl list a while ago, but
haven't received a reply...
Is this sort of contribution useful to the project?
Thanks,
Sam
---------- Forwarded message ----------
From: Sam Mason <s.a.mason@warwick.ac.uk>
Date: 15 June 2012 17:59
To: help-gsl@gnu.org
Hi,
I sent a message to this a couple of months ago enquiring about the
availability of log-density functions. I've been using my own
versions of these since then, and thought it would be good to send
them back to the community.
I've gone through most of the distributions and added in code
calculating their log-densities. I've copied the naming and layout
convention of the dirichlet distribution; which already had the
gsl_ran_dirichlet_pdf() and gsl_ran_dirichlet_lnpdf() functions. I've
not updated the docs yet, as I wanted to gage interest first.
The code compiles and the tests run successfully on my computer. I've
also used the test code as a mini-benchmark, and the changed code
takes basically the same amount of time to execute (median = 1.0038
times GSL 1.15, 95%CI=1.0009-1.0064).
Let me know what would make this more useful.
Sam
[-- Attachment #2: log-pdfs-v1.patch --]
[-- Type: application/octet-stream, Size: 32548 bytes --]
diff -u gsl-1.15/randist/beta.c gsl-sam/randist/beta.c
--- gsl-1.15/randist/beta.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/beta.c 2012-06-15 15:55:57.000000000 +0100
@@ -42,9 +42,15 @@
double
gsl_ran_beta_pdf (const double x, const double a, const double b)
{
+ return exp(gsl_ran_beta_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_beta_lnpdf (const double x, const double a, const double b)
+{
if (x < 0 || x > 1)
{
- return 0 ;
+ return GSL_NEGINF;
}
else
{
@@ -54,20 +60,13 @@
double ga = gsl_sf_lngamma (a);
double gb = gsl_sf_lngamma (b);
- if (x == 0.0 || x == 1.0)
- {
- if (a > 1.0 && b > 1.0)
- {
- p = 0.0;
- }
- else
- {
- p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
- }
- }
+ if ((x == 0.0 || x == 1.0) && a > 1.0 && b > 1.0)
+ {
+ p = GSL_NEGINF;
+ }
else
{
- p = exp (gab - ga - gb + log(x) * (a - 1) + log1p(-x) * (b - 1));
+ p = gab - ga - gb + log (x) * (a - 1) + log1p (-x) * (b - 1);
}
return p;
diff -u gsl-1.15/randist/bigauss.c gsl-sam/randist/bigauss.c
--- gsl-1.15/randist/bigauss.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/bigauss.c 2012-06-15 13:38:21.000000000 +0100
@@ -61,10 +61,18 @@
const double sigma_x, const double sigma_y,
const double rho)
{
+ return exp (gsl_ran_bivariate_gaussian_lnpdf (x, y, sigma_x, sigma_y, rho));
+}
+
+double
+gsl_ran_bivariate_gaussian_lnpdf (const double x, const double y,
+ const double sigma_x, const double sigma_y,
+ const double rho)
+{
double u = x / sigma_x ;
double v = y / sigma_y ;
double c = 1 - rho*rho ;
- double p = (1 / (2 * M_PI * sigma_x * sigma_y * sqrt(c)))
- * exp (-(u * u - 2 * rho * u * v + v * v) / (2 * c));
+ double p = -log(2 * M_PI * sigma_x * sigma_y * sqrt(c))
+ - (u * u - 2 * rho * u * v + v * v) / (2 * c);
return p;
}
diff -u gsl-1.15/randist/binomial.c gsl-sam/randist/binomial.c
--- gsl-1.15/randist/binomial.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/binomial.c 2012-06-15 16:21:17.000000000 +0100
@@ -23,6 +23,7 @@
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_math.h>
/* The binomial distribution has the form,
@@ -72,9 +73,16 @@
gsl_ran_binomial_pdf (const unsigned int k, const double p,
const unsigned int n)
{
+ return exp (gsl_ran_binomial_lnpdf (k, p, n));
+}
+
+double
+gsl_ran_binomial_lnpdf (const unsigned int k, const double p,
+ const unsigned int n)
+{
if (k > n)
{
- return 0;
+ return GSL_NEGINF;
}
else
{
@@ -82,17 +90,16 @@
if (p == 0)
{
- P = (k == 0) ? 1 : 0;
+ P = (k == 0) ? 0 : GSL_NEGINF;
}
else if (p == 1)
{
- P = (k == n) ? 1 : 0;
+ P = (k == n) ? 0 : GSL_NEGINF;
}
else
{
double ln_Cnk = gsl_sf_lnchoose (n, k);
P = ln_Cnk + k * log (p) + (n - k) * log1p (-p);
- P = exp (P);
}
return P;
diff -u gsl-1.15/randist/cauchy.c gsl-sam/randist/cauchy.c
--- gsl-1.15/randist/cauchy.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/cauchy.c 2012-06-15 15:00:07.000000000 +0100
@@ -45,7 +45,13 @@
double
gsl_ran_cauchy_pdf (const double x, const double a)
{
+ return exp (gsl_ran_cauchy_lnpdf (x, a));
+}
+
+double
+gsl_ran_cauchy_lnpdf (const double x, const double a)
+{
double u = x / a;
- double p = (1 / (M_PI * a)) / (1 + u * u);
+ double p = -log(M_PI * a) - log1p(u * u);
return p;
}
diff -u gsl-1.15/randist/chisq.c gsl-sam/randist/chisq.c
--- gsl-1.15/randist/chisq.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/chisq.c 2012-06-15 16:21:23.000000000 +0100
@@ -22,6 +22,7 @@
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The chisq distribution has the form
@@ -39,22 +40,28 @@
double
gsl_ran_chisq_pdf (const double x, const double nu)
{
+ return exp (gsl_ran_chisq_lnpdf (x, nu));
+}
+
+double
+gsl_ran_chisq_lnpdf (const double x, const double nu)
+{
if (x < 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
if(nu == 2.0)
{
- return exp(-x/2.0) / 2.0;
+ return (-x/2.0) - log (2.0);
}
else
{
double p;
double lngamma = gsl_sf_lngamma (nu / 2);
- p = exp ((nu / 2 - 1) * log (x/2) - x/2 - lngamma) / 2;
+ p = (nu / 2 - 1) * log (x/2) - x/2 - lngamma - log (2.0);
return p;
}
}
diff -u gsl-1.15/randist/erlang.c gsl-sam/randist/erlang.c
--- gsl-1.15/randist/erlang.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/erlang.c 2012-06-15 16:21:53.000000000 +0100
@@ -22,6 +22,7 @@
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The sum of N samples from an exponential distribution gives an
Erlang distribution
@@ -39,16 +40,22 @@
double
gsl_ran_erlang_pdf (const double x, const double a, const double n)
{
+ return exp (gsl_ran_erlang_lnpdf (x, a, n));
+}
+
+double
+gsl_ran_erlang_lnpdf (const double x, const double a, const double n)
+{
if (x <= 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
double p;
double lngamma = gsl_sf_lngamma (n);
- p = exp ((n - 1) * log (x/a) - x/a - lngamma) / a;
+ p = (n - 1) * log (x/a) - x/a - lngamma - log (a);
return p;
}
}
diff -u gsl-1.15/randist/exponential.c gsl-sam/randist/exponential.c
--- gsl-1.15/randist/exponential.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/exponential.c 2012-06-15 16:02:41.000000000 +0100
@@ -40,13 +40,19 @@
double
gsl_ran_exponential_pdf (const double x, const double mu)
{
+ return exp (gsl_ran_exponential_lnpdf (x, mu));
+}
+
+double
+gsl_ran_exponential_lnpdf (const double x, const double mu)
+{
if (x < 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
- double p = exp (-x/mu)/mu;
+ double p = (-x/mu) - log (mu);
return p;
}
diff -u gsl-1.15/randist/exppow.c gsl-sam/randist/exppow.c
--- gsl-1.15/randist/exppow.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/exppow.c 2012-06-15 14:59:19.000000000 +0100
@@ -115,8 +115,14 @@
double
gsl_ran_exppow_pdf (const double x, const double a, const double b)
{
+ return exp (gsl_ran_exppow_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_exppow_lnpdf (const double x, const double a, const double b)
+{
double p;
double lngamma = gsl_sf_lngamma (1 + 1 / b);
- p = (1 / (2 * a)) * exp (-pow (fabs (x / a), b) - lngamma);
+ p = -log(2 * a) + (-pow (fabs (x / a), b) - lngamma);
return p;
}
diff -u gsl-1.15/randist/fdist.c gsl-sam/randist/fdist.c
--- gsl-1.15/randist/fdist.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/fdist.c 2012-06-15 16:22:01.000000000 +0100
@@ -22,6 +22,7 @@
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The F distribution has the form
@@ -46,9 +47,15 @@
double
gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2)
{
+ return exp (gsl_ran_fdist_lnpdf (x, nu1, nu2));
+}
+
+double
+gsl_ran_fdist_lnpdf (const double x, const double nu1, const double nu2)
+{
if (x < 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
@@ -59,9 +66,8 @@
double lg1 = gsl_sf_lngamma (nu1 / 2);
double lg2 = gsl_sf_lngamma (nu2 / 2);
- p =
- exp (lglg + lg12 - lg1 - lg2 + (nu1 / 2 - 1) * log (x) -
- ((nu1 + nu2) / 2) * log (nu2 + nu1 * x));
+ p = (lglg + lg12 - lg1 - lg2 + (nu1 / 2 - 1) * log (x) -
+ ((nu1 + nu2) / 2) * log (nu2 + nu1 * x));
return p;
}
diff -u gsl-1.15/randist/flat.c gsl-sam/randist/flat.c
--- gsl-1.15/randist/flat.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/flat.c 2012-06-15 16:22:08.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* This is the uniform distribution in the range [a, b)
@@ -42,12 +43,18 @@
double
gsl_ran_flat_pdf (double x, const double a, const double b)
{
+ return exp (gsl_ran_flat_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_flat_lnpdf (double x, const double a, const double b)
+{
if (x < b && x >= a)
{
- return 1 / (b - a);
+ return -log (b - a);
}
else
{
- return 0;
+ return GSL_NEGINF;
}
}
diff -u gsl-1.15/randist/gamma.c gsl-sam/randist/gamma.c
--- gsl-1.15/randist/gamma.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gamma.c 2012-06-15 12:34:34.000000000 +0100
@@ -150,26 +150,32 @@
double
gsl_ran_gamma_pdf (const double x, const double a, const double b)
{
+ return exp (gsl_ran_gamma_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_gamma_lnpdf (const double x, const double a, const double b)
+{
if (x < 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else if (x == 0)
{
if (a == 1)
- return 1/b ;
+ return -log(b) ;
else
- return 0 ;
+ return GSL_NEGINF ;
}
else if (a == 1)
{
- return exp(-x/b)/b ;
+ return (-x/b) - log(b) ;
}
else
{
double p;
double lngamma = gsl_sf_lngamma (a);
- p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
+ p = (a - 1) * log (x/b) - x/b - lngamma - log (b);
return p;
}
}
diff -u gsl-1.15/randist/gauss.c gsl-sam/randist/gauss.c
--- gsl-1.15/randist/gauss.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gauss.c 2012-06-15 13:31:48.000000000 +0100
@@ -118,8 +118,14 @@
double
gsl_ran_gaussian_pdf (const double x, const double sigma)
{
+ return exp (gsl_ran_gaussian_lnpdf (x, sigma));
+}
+
+double
+gsl_ran_gaussian_lnpdf (const double x, const double sigma)
+{
double u = x / fabs (sigma);
- double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
+ double p = -log(sqrt (2 * M_PI) * fabs (sigma)) + (-u * u / 2);
return p;
}
@@ -136,8 +142,13 @@
}
double
-gsl_ran_ugaussian_pdf (const double x)
+gsl_ran_ugaussian_lnpdf (const double x)
{
- return gsl_ran_gaussian_pdf (x, 1.0);
+ return gsl_ran_gaussian_lnpdf (x, 1.0);
}
+double
+gsl_ran_ugaussian_pdf (const double x)
+{
+ return exp (gsl_ran_gaussian_lnpdf (x, 1.0));
+}
diff -u gsl-1.15/randist/gausstail.c gsl-sam/randist/gausstail.c
--- gsl-1.15/randist/gausstail.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gausstail.c 2012-06-15 13:34:07.000000000 +0100
@@ -75,9 +75,15 @@
double
gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma)
{
+ return exp (gsl_ran_gaussian_tail_lnpdf (x, a, sigma));
+}
+
+double
+gsl_ran_gaussian_tail_lnpdf (const double x, const double a, const double sigma)
+{
if (x < a)
{
- return 0;
+ return GSL_NEGINF;
}
else
{
@@ -88,7 +94,7 @@
N = 0.5 * f;
- p = (1 / (N * sqrt (2 * M_PI) * sigma)) * exp (-u * u / 2);
+ p = -log(N * sqrt (2 * M_PI) * sigma) + (-u * u / 2);
return p;
}
@@ -103,5 +109,11 @@
double
gsl_ran_ugaussian_tail_pdf (const double x, const double a)
{
- return gsl_ran_gaussian_tail_pdf (x, a, 1.0) ;
+ return exp (gsl_ran_gaussian_tail_lnpdf (x, a, 1.0)) ;
+}
+
+double
+gsl_ran_ugaussian_tail_lnpdf (const double x, const double a)
+{
+ return gsl_ran_gaussian_tail_lnpdf (x, a, 1.0) ;
}
diff -u gsl-1.15/randist/geometric.c gsl-sam/randist/geometric.c
--- gsl-1.15/randist/geometric.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/geometric.c 2012-06-15 16:22:18.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* Geometric distribution (bernoulli trial with probability p)
@@ -51,17 +52,23 @@
double
gsl_ran_geometric_pdf (const unsigned int k, const double p)
{
+ return exp (gsl_ran_geometric_lnpdf (k, p));
+}
+
+double
+gsl_ran_geometric_lnpdf (const unsigned int k, const double p)
+{
if (k == 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else if (k == 1)
{
- return p ;
+ return log (p) ;
}
else
{
- double P = p * pow (1 - p, k - 1.0);
+ double P = log (p) + log1p (-p) * (k - 1.0);
return P;
}
}
diff -u gsl-1.15/randist/gsl_randist.h gsl-sam/randist/gsl_randist.h
--- gsl-1.15/randist/gsl_randist.h 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gsl_randist.h 2012-06-15 15:49:21.000000000 +0100
@@ -38,23 +38,29 @@
double gsl_ran_beta (const gsl_rng * r, const double a, const double b);
double gsl_ran_beta_pdf (const double x, const double a, const double b);
+double gsl_ran_beta_lnpdf (const double x, const double a, const double b);
unsigned int gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n);
unsigned int gsl_ran_binomial_knuth (const gsl_rng * r, double p, unsigned int n);
unsigned int gsl_ran_binomial_tpe (const gsl_rng * r, double p, unsigned int n);
double gsl_ran_binomial_pdf (const unsigned int k, const double p, const unsigned int n);
+double gsl_ran_binomial_lnpdf (const unsigned int k, const double p, const unsigned int n);
double gsl_ran_exponential (const gsl_rng * r, const double mu);
double gsl_ran_exponential_pdf (const double x, const double mu);
+double gsl_ran_exponential_lnpdf (const double x, const double mu);
double gsl_ran_exppow (const gsl_rng * r, const double a, const double b);
double gsl_ran_exppow_pdf (const double x, const double a, const double b);
+double gsl_ran_exppow_lnpdf (const double x, const double a, const double b);
double gsl_ran_cauchy (const gsl_rng * r, const double a);
double gsl_ran_cauchy_pdf (const double x, const double a);
+double gsl_ran_cauchy_lnpdf (const double x, const double a);
double gsl_ran_chisq (const gsl_rng * r, const double nu);
double gsl_ran_chisq_pdf (const double x, const double nu);
+double gsl_ran_chisq_lnpdf (const double x, const double nu);
void gsl_ran_dirichlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
double gsl_ran_dirichlet_pdf (const size_t K, const double alpha[], const double theta[]);
@@ -62,16 +68,20 @@
double gsl_ran_erlang (const gsl_rng * r, const double a, const double n);
double gsl_ran_erlang_pdf (const double x, const double a, const double n);
+double gsl_ran_erlang_lnpdf (const double x, const double a, const double n);
double gsl_ran_fdist (const gsl_rng * r, const double nu1, const double nu2);
double gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2);
+double gsl_ran_fdist_lnpdf (const double x, const double nu1, const double nu2);
double gsl_ran_flat (const gsl_rng * r, const double a, const double b);
double gsl_ran_flat_pdf (double x, const double a, const double b);
+double gsl_ran_flat_lnpdf (double x, const double a, const double b);
double gsl_ran_gamma (const gsl_rng * r, const double a, const double b);
double gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a);
double gsl_ran_gamma_pdf (const double x, const double a, const double b);
+double gsl_ran_gamma_lnpdf (const double x, const double a, const double b);
double gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b);
double gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b);
@@ -79,43 +89,55 @@
double gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma);
double gsl_ran_gaussian_ziggurat (const gsl_rng * r, const double sigma);
double gsl_ran_gaussian_pdf (const double x, const double sigma);
+double gsl_ran_gaussian_lnpdf (const double x, const double sigma);
double gsl_ran_ugaussian (const gsl_rng * r);
double gsl_ran_ugaussian_ratio_method (const gsl_rng * r);
double gsl_ran_ugaussian_pdf (const double x);
+double gsl_ran_ugaussian_lnpdf (const double x);
double gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma);
double gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma);
+double gsl_ran_gaussian_tail_lnpdf (const double x, const double a, const double sigma);
double gsl_ran_ugaussian_tail (const gsl_rng * r, const double a);
double gsl_ran_ugaussian_tail_pdf (const double x, const double a);
+double gsl_ran_ugaussian_tail_lnpdf (const double x, const double a);
void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y);
double gsl_ran_bivariate_gaussian_pdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);
+double gsl_ran_bivariate_gaussian_lnpdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);
double gsl_ran_landau (const gsl_rng * r);
double gsl_ran_landau_pdf (const double x);
unsigned int gsl_ran_geometric (const gsl_rng * r, const double p);
double gsl_ran_geometric_pdf (const unsigned int k, const double p);
+double gsl_ran_geometric_lnpdf (const unsigned int k, const double p);
unsigned int gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2, unsigned int t);
double gsl_ran_hypergeometric_pdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);
+double gsl_ran_hypergeometric_lnpdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);
double gsl_ran_gumbel1 (const gsl_rng * r, const double a, const double b);
double gsl_ran_gumbel1_pdf (const double x, const double a, const double b);
+double gsl_ran_gumbel1_lnpdf (const double x, const double a, const double b);
double gsl_ran_gumbel2 (const gsl_rng * r, const double a, const double b);
double gsl_ran_gumbel2_pdf (const double x, const double a, const double b);
+double gsl_ran_gumbel2_lnpdf (const double x, const double a, const double b);
double gsl_ran_logistic (const gsl_rng * r, const double a);
double gsl_ran_logistic_pdf (const double x, const double a);
+double gsl_ran_logistic_lnpdf (const double x, const double a);
double gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma);
double gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma);
+double gsl_ran_lognormal_lnpdf (const double x, const double zeta, const double sigma);
unsigned int gsl_ran_logarithmic (const gsl_rng * r, const double p);
double gsl_ran_logarithmic_pdf (const unsigned int k, const double p);
+double gsl_ran_logarithmic_lnpdf (const unsigned int k, const double p);
void gsl_ran_multinomial (const gsl_rng * r, const size_t K,
const unsigned int N, const double p[],
@@ -128,35 +150,44 @@
unsigned int gsl_ran_negative_binomial (const gsl_rng * r, double p, double n);
double gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n);
+double gsl_ran_negative_binomial_lnpdf (const unsigned int k, const double p, double n);
unsigned int gsl_ran_pascal (const gsl_rng * r, double p, unsigned int n);
double gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n);
+double gsl_ran_pascal_lnpdf (const unsigned int k, const double p, unsigned int n);
double gsl_ran_pareto (const gsl_rng * r, double a, const double b);
double gsl_ran_pareto_pdf (const double x, const double a, const double b);
+double gsl_ran_pareto_lnpdf (const double x, const double a, const double b);
unsigned int gsl_ran_poisson (const gsl_rng * r, double mu);
void gsl_ran_poisson_array (const gsl_rng * r, size_t n, unsigned int array[],
double mu);
double gsl_ran_poisson_pdf (const unsigned int k, const double mu);
+double gsl_ran_poisson_lnpdf (const unsigned int k, const double mu);
double gsl_ran_rayleigh (const gsl_rng * r, const double sigma);
double gsl_ran_rayleigh_pdf (const double x, const double sigma);
+double gsl_ran_rayleigh_lnpdf (const double x, const double sigma);
double gsl_ran_rayleigh_tail (const gsl_rng * r, const double a, const double sigma);
double gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma);
+double gsl_ran_rayleigh_tail_lnpdf (const double x, const double a, const double sigma);
double gsl_ran_tdist (const gsl_rng * r, const double nu);
double gsl_ran_tdist_pdf (const double x, const double nu);
+double gsl_ran_tdist_lnpdf (const double x, const double nu);
double gsl_ran_laplace (const gsl_rng * r, const double a);
double gsl_ran_laplace_pdf (const double x, const double a);
+double gsl_ran_laplace_lnpdf (const double x, const double a);
double gsl_ran_levy (const gsl_rng * r, const double c, const double alpha);
double gsl_ran_levy_skew (const gsl_rng * r, const double c, const double alpha, const double beta);
double gsl_ran_weibull (const gsl_rng * r, const double a, const double b);
double gsl_ran_weibull_pdf (const double x, const double a, const double b);
+double gsl_ran_weibull_lnpdf (const double x, const double a, const double b);
void gsl_ran_dir_2d (const gsl_rng * r, double * x, double * y);
void gsl_ran_dir_2d_trig_method (const gsl_rng * r, double * x, double * y);
diff -u gsl-1.15/randist/gumbel.c gsl-sam/randist/gumbel.c
--- gsl-1.15/randist/gumbel.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gumbel.c 2012-06-15 16:22:29.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The Type I Gumbel distribution has the form,
@@ -45,7 +46,13 @@
double
gsl_ran_gumbel1_pdf (const double x, const double a, const double b)
{
- double p = a * b * exp (-(b * exp(-a * x) + a * x));
+ return exp (gsl_ran_gumbel1_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_gumbel1_lnpdf (const double x, const double a, const double b)
+{
+ double p = log(a) + log(b) - (b * exp(-a * x) + a * x);
return p;
}
@@ -62,13 +69,19 @@
double
gsl_ran_gumbel2_pdf (const double x, const double a, const double b)
{
+ return exp (gsl_ran_gumbel2_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_gumbel2_lnpdf (const double x, const double a, const double b)
+{
if (x <= 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
- double p = b * a * pow(x,-(a+1)) * exp (-b * pow(x, -a));
+ double p = log(b) + log(a) - log(x) * (a+1) - (b * pow(x, -a));
return p;
}
}
diff -u gsl-1.15/randist/hyperg.c gsl-sam/randist/hyperg.c
--- gsl-1.15/randist/hyperg.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/hyperg.c 2012-06-15 16:22:39.000000000 +0100
@@ -22,6 +22,7 @@
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_math.h>
/* The hypergeometric distribution has the form,
@@ -96,6 +97,15 @@
const unsigned int n2,
unsigned int t)
{
+ return exp (gsl_ran_hypergeometric_lnpdf (k, n1, n2, t));
+}
+
+double
+gsl_ran_hypergeometric_lnpdf (const unsigned int k,
+ const unsigned int n1,
+ const unsigned int n2,
+ unsigned int t)
+{
if (t > n1 + n2)
{
t = n1 + n2 ;
@@ -103,11 +113,11 @@
if (k > n1 || k > t)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else if (t > n2 && k + n2 < t )
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
@@ -117,7 +127,7 @@
double c2 = gsl_sf_lnchoose(n2,t-k);
double c3 = gsl_sf_lnchoose(n1+n2,t);
- p = exp(c1 + c2 - c3) ;
+ p = c1 + c2 - c3 ;
return p;
}
diff -u gsl-1.15/randist/laplace.c gsl-sam/randist/laplace.c
--- gsl-1.15/randist/laplace.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/laplace.c 2012-06-15 16:22:59.000000000 +0100
@@ -51,7 +51,14 @@
double
gsl_ran_laplace_pdf (const double x, const double a)
{
- double p = (1/(2*a)) * exp (-fabs (x)/a);
+ return exp (gsl_ran_laplace_lnpdf (x, a));
+}
+
+double
+gsl_ran_laplace_lnpdf (const double x, const double a)
+{
+ double p = -log(2*a) - (fabs (x)/a);
+
return p;
}
diff -u gsl-1.15/randist/logarithmic.c gsl-sam/randist/logarithmic.c
--- gsl-1.15/randist/logarithmic.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/logarithmic.c 2012-06-15 16:23:18.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* Logarithmic distribution
@@ -64,13 +65,19 @@
double
gsl_ran_logarithmic_pdf (const unsigned int k, const double p)
{
+ return exp (gsl_ran_logarithmic_lnpdf (k, p));
+}
+
+double
+gsl_ran_logarithmic_lnpdf (const unsigned int k, const double p)
+{
if (k == 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
- double P = pow(p, (double)k) / (double) k / log(1/(1-p)) ;
+ double P = log(p) * k - log(k) - log(-log1p(-p)) ;
return P;
}
}
diff -u gsl-1.15/randist/logistic.c gsl-sam/randist/logistic.c
--- gsl-1.15/randist/logistic.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/logistic.c 2012-06-15 14:20:03.000000000 +0100
@@ -47,7 +47,13 @@
double
gsl_ran_logistic_pdf (const double x, const double a)
{
- double u = exp (-fabs(x)/a);
- double p = u / (fabs(a) * (1 + u) * (1 + u));
+ return exp (gsl_ran_logistic_lnpdf (x, a));
+}
+
+double
+gsl_ran_logistic_lnpdf (const double x, const double a)
+{
+ double u = -fabs(x)/a;
+ double p = u - (log(fabs(a)) + log1p (exp (u)) * 2);
return p;
}
diff -u gsl-1.15/randist/lognormal.c gsl-sam/randist/lognormal.c
--- gsl-1.15/randist/lognormal.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/lognormal.c 2012-06-15 16:14:53.000000000 +0100
@@ -57,14 +57,20 @@
double
gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma)
{
+ return exp (gsl_ran_lognormal_lnpdf (x, zeta, sigma));
+}
+
+double
+gsl_ran_lognormal_lnpdf (const double x, const double zeta, const double sigma)
+{
if (x <= 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
double u = (log (x) - zeta)/sigma;
- double p = 1 / (x * fabs(sigma) * sqrt (2 * M_PI)) * exp (-(u * u) /2);
+ double p = -log(x * fabs(sigma) * sqrt (2 * M_PI)) - (u * u) / 2;
return p;
}
}
diff -u gsl-1.15/randist/nbinomial.c gsl-sam/randist/nbinomial.c
--- gsl-1.15/randist/nbinomial.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/nbinomial.c 2012-06-15 14:46:32.000000000 +0100
@@ -42,13 +42,19 @@
double
gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n)
{
+ return exp( gsl_ran_negative_binomial_lnpdf (k, p, n));
+}
+
+double
+gsl_ran_negative_binomial_lnpdf (const unsigned int k, const double p, double n)
+{
double P;
double f = gsl_sf_lngamma (k + n) ;
double a = gsl_sf_lngamma (n) ;
double b = gsl_sf_lngamma (k + 1.0) ;
- P = exp(f-a-b) * pow (p, n) * pow (1 - p, (double)k);
+ P = (f-a-b) + log (p) * n + log1p (-p) * k;
return P;
}
diff -u gsl-1.15/randist/pareto.c gsl-sam/randist/pareto.c
--- gsl-1.15/randist/pareto.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/pareto.c 2012-06-15 16:23:26.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The Pareto distribution has the form,
@@ -41,14 +42,21 @@
double
gsl_ran_pareto_pdf (const double x, const double a, const double b)
{
+ return exp (gsl_ran_pareto_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_pareto_lnpdf (const double x, const double a, const double b)
+{
if (x >= b)
{
- double p = (a/b) / pow (x/b, a + 1);
+ double lb = log(b);
+ double p = (log(a) - lb) - (log (x) - lb) * (a + 1);
return p;
}
else
{
- return 0;
+ return GSL_NEGINF;
}
}
diff -u gsl-1.15/randist/pascal.c gsl-sam/randist/pascal.c
--- gsl-1.15/randist/pascal.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/pascal.c 2012-06-15 14:49:27.000000000 +0100
@@ -45,6 +45,11 @@
double
gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n)
{
- double P = gsl_ran_negative_binomial_pdf (k, p, (double) n);
- return P;
+ return exp (gsl_ran_negative_binomial_lnpdf (k, p, (double) n));
+}
+
+double
+gsl_ran_pascal_lnpdf (const unsigned int k, const double p, unsigned int n)
+{
+ return gsl_ran_negative_binomial_lnpdf (k, p, (double) n);
}
diff -u gsl-1.15/randist/poisson.c gsl-sam/randist/poisson.c
--- gsl-1.15/randist/poisson.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/poisson.c 2012-06-15 14:57:26.000000000 +0100
@@ -85,9 +85,15 @@
double
gsl_ran_poisson_pdf (const unsigned int k, const double mu)
{
+ return exp(gsl_ran_poisson_lnpdf (k, mu));
+}
+
+double
+gsl_ran_poisson_lnpdf (const unsigned int k, const double mu)
+{
double p;
double lf = gsl_sf_lnfact (k);
- p = exp (log (mu) * k - lf - mu);
+ p = log (mu) * k - lf - mu;
return p;
}
diff -u gsl-1.15/randist/rayleigh.c gsl-sam/randist/rayleigh.c
--- gsl-1.15/randist/rayleigh.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/rayleigh.c 2012-06-15 16:23:34.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The Rayleigh distribution has the form
@@ -39,14 +40,20 @@
double
gsl_ran_rayleigh_pdf (const double x, const double sigma)
{
+ return exp (gsl_ran_rayleigh_lnpdf (x, sigma));
+}
+
+double
+gsl_ran_rayleigh_lnpdf (const double x, const double sigma)
+{
if (x < 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
double u = x / sigma ;
- double p = (u / sigma) * exp(-u * u / 2.0) ;
+ double p = log(u / sigma) - (u * u / 2.0) ;
return p;
}
@@ -69,17 +76,23 @@
double
gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma)
{
+ return exp (gsl_ran_rayleigh_tail_lnpdf (x, a, sigma));
+}
+
+double
+gsl_ran_rayleigh_tail_lnpdf (const double x, const double a, const double sigma)
+{
if (x < a)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else
{
double u = x / sigma ;
double v = a / sigma ;
- double p = (u / sigma) * exp((v + u) * (v - u) / 2.0) ;
-
+ double p = log(u / sigma) + ((v + u) * (v - u) / 2.0) ;
+
return p;
}
}
diff -u gsl-1.15/randist/tdist.c gsl-sam/randist/tdist.c
--- gsl-1.15/randist/tdist.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/tdist.c 2012-06-15 17:18:00.000000000 +0100
@@ -67,13 +67,20 @@
double
gsl_ran_tdist_pdf (const double x, const double nu)
{
+ return exp (gsl_ran_tdist_lnpdf (x, nu));
+}
+
+double
+gsl_ran_tdist_lnpdf (const double x, const double nu)
+{
double p;
double lg1 = gsl_sf_lngamma (nu / 2);
double lg2 = gsl_sf_lngamma ((nu + 1) / 2);
+
+ p = (((lg2 - lg1) - log(sqrt (M_PI * nu)))
+ + log1p (x * x / nu) * ((-1 - nu) / 2));
- p = ((exp (lg2 - lg1) / sqrt (M_PI * nu))
- * pow ((1 + x * x / nu), -(nu + 1) / 2));
return p;
}
diff -u gsl-1.15/randist/weibull.c gsl-sam/randist/weibull.c
--- gsl-1.15/randist/weibull.c 2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/weibull.c 2012-06-15 16:23:46.000000000 +0100
@@ -21,6 +21,7 @@
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
/* The Weibull distribution has the form,
@@ -41,24 +42,30 @@
double
gsl_ran_weibull_pdf (const double x, const double a, const double b)
{
+ return exp (gsl_ran_weibull_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_weibull_lnpdf (const double x, const double a, const double b)
+{
if (x < 0)
{
- return 0 ;
+ return GSL_NEGINF ;
}
else if (x == 0)
{
if (b == 1)
- return 1/a ;
+ return -log(a) ;
else
- return 0 ;
+ return GSL_NEGINF ;
}
else if (b == 1)
{
- return exp(-x/a)/a ;
+ return (-x/a) - log(a) ;
}
else
{
- double p = (b/a) * exp (-pow (x/a, b) + (b - 1) * log (x/a));
+ double p = log(b/a) + (-pow (x/a, b) + (b - 1) * log (x/a));
return p;
}
}
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: Fwd: log probability density functions (PDFs)
2012-07-09 17:30 ` Fwd: log probability density functions (PDFs) Sam Mason
@ 2012-07-09 17:39 ` Patrick Alken
[not found] ` <d75321b0fc2d4b7fb010d63b6252dd89@IDMEXCHT02.ads.warwick.ac.uk>
1 sibling, 0 replies; 3+ messages in thread
From: Patrick Alken @ 2012-07-09 17:39 UTC (permalink / raw)
To: gsl-discuss
The best option would be to make a standalone extension (see the main
gsl page for examples), so people can look at it and try it out. We can
add the extension to the main web page and possibly add it to the main
source code later if it looks useful.
On 07/09/2012 11:29 AM, Sam Mason wrote:
> Hi,
>
> I sent the following message to the help-gsl list a while ago, but
> haven't received a reply...
>
> Is this sort of contribution useful to the project?
>
> Thanks,
>
> Sam
>
>
> ---------- Forwarded message ----------
> From: Sam Mason <s.a.mason@warwick.ac.uk>
> Date: 15 June 2012 17:59
> To: help-gsl@gnu.org
>
>
> Hi,
>
> I sent a message to this a couple of months ago enquiring about the
> availability of log-density functions. I've been using my own
> versions of these since then, and thought it would be good to send
> them back to the community.
>
> I've gone through most of the distributions and added in code
> calculating their log-densities. I've copied the naming and layout
> convention of the dirichlet distribution; which already had the
> gsl_ran_dirichlet_pdf() and gsl_ran_dirichlet_lnpdf() functions. I've
> not updated the docs yet, as I wanted to gage interest first.
>
> The code compiles and the tests run successfully on my computer. I've
> also used the test code as a mini-benchmark, and the changed code
> takes basically the same amount of time to execute (median = 1.0038
> times GSL 1.15, 95%CI=1.0009-1.0064).
>
> Let me know what would make this more useful.
>
> Sam
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: Fwd: log probability density functions (PDFs)
[not found] ` <d75321b0fc2d4b7fb010d63b6252dd89@IDMEXCHT02.ads.warwick.ac.uk>
@ 2012-07-09 21:05 ` Sam Mason
0 siblings, 0 replies; 3+ messages in thread
From: Sam Mason @ 2012-07-09 21:05 UTC (permalink / raw)
To: Patrick Alken; +Cc: gsl-discuss
Hi Patrick,
I've put the following together:
http://www.samason.me.uk/hacking-the-gsl
Is that the sort of thing you'd consider linking to? Maybe with an
accompanying sentence of "log-probability density functions"?
Thanks,
Sam
On 9 July 2012 18:39, Patrick Alken <patrick.alken@colorado.edu> wrote:
> The best option would be to make a standalone extension (see the main
> gsl page for examples), so people can look at it and try it out. We can
> add the extension to the main web page and possibly add it to the main
> source code later if it looks useful.
>
> On 07/09/2012 11:29 AM, Sam Mason wrote:
>> Hi,
>>
>> I sent the following message to the help-gsl list a while ago, but
>> haven't received a reply...
>>
>> Is this sort of contribution useful to the project?
>>
>> Thanks,
>>
>> Sam
>>
>>
>> ---------- Forwarded message ----------
>> From: Sam Mason <s.a.mason@warwick.ac.uk>
>> Date: 15 June 2012 17:59
>> To: help-gsl@gnu.org
>>
>>
>> Hi,
>>
>> I sent a message to this a couple of months ago enquiring about the
>> availability of log-density functions. I've been using my own
>> versions of these since then, and thought it would be good to send
>> them back to the community.
>>
>> I've gone through most of the distributions and added in code
>> calculating their log-densities. I've copied the naming and layout
>> convention of the dirichlet distribution; which already had the
>> gsl_ran_dirichlet_pdf() and gsl_ran_dirichlet_lnpdf() functions. I've
>> not updated the docs yet, as I wanted to gage interest first.
>>
>> The code compiles and the tests run successfully on my computer. I've
>> also used the test code as a mini-benchmark, and the changed code
>> takes basically the same amount of time to execute (median = 1.0038
>> times GSL 1.15, 95%CI=1.0009-1.0064).
>>
>> Let me know what would make this more useful.
>>
>> Sam
>
>
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2012-07-09 21:05 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
[not found] <CAPEJTf6T8vk+KgnoDBaZgJzjRUvJfLK+w-S-xwoZ_=f_wgNH3g@mail.gmail.com>
2012-07-09 17:30 ` Fwd: log probability density functions (PDFs) Sam Mason
2012-07-09 17:39 ` Patrick Alken
[not found] ` <d75321b0fc2d4b7fb010d63b6252dd89@IDMEXCHT02.ads.warwick.ac.uk>
2012-07-09 21:05 ` Sam Mason
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).