Permalink
Browse files

Adding ln versions of gaussian and gamma, fixing some function signat…

…ures for compatibility
  • Loading branch information...
1 parent 935f262 commit d5788548b2fb7ff94e4ccbfecbcf206629548a73 @jerrytalton jerrytalton committed Sep 28, 2011
Showing with 24 additions and 9 deletions.
  1. +19 −6 scheme-compilers/javascript/math-functions.js
  2. +5 −3 scheme-compilers/scheme2js-template.sc
View
25 scheme-compilers/javascript/math-functions.js
@@ -58,7 +58,7 @@ function sample_gamma(a,b)
while(true)
{
- do{x = sample_gaussian(1); v = 1+c*x;} while(v <= 0);
+ do{x = sample_gaussian(0,1); v = 1+c*x;} while(v <= 0);
v=v*v*v;
u=random();
@@ -77,6 +77,12 @@ function gamma_pdf(x,a,b)
return Math.exp((a - 1)*Math.log(x/b) - x/b - log_gamma(a))/b;
}
+// Evaluate log gammma pdf
+function gamma_lnpdf(x,a,b)
+{
+ return (1 - a)*Math.log(x) - x/b - log_gamma(a) - a*Math.log(b);
+}
+
// Draw a sample from a Binomial distribution
// Knuth TAOCP 2 (could be improved, i.e. via Kachitvichyanukul & Schmeiser)
function sample_binomial(p,n)
@@ -142,7 +148,7 @@ function sample_beta(a, b)
// Draw a sample from a Gaussian distribution
// Leva '92 (could be improved, i.e. via Ziggurat method)
-function sample_gaussian(sigma)
+function sample_gaussian(mu,sigma)
{
var u, v, x, y, q;
@@ -156,17 +162,24 @@ function sample_gaussian(sigma)
}
while(q >= 0.27597 && (q > 0.27846 || v*v > -4 * u * u * Math.log(u)))
- return sigma*v/u;
+ return mu + sigma*v/u;
}
// Evaluate the gaussian distribution
-function gaussian_pdf(x,sigma)
+function gaussian_pdf(x,mu,sigma)
{
+ x-=mu;
var asigma = Math.abs(sigma);
var u = x/asigma;
return (1/ Math.sqrt(2*Math.PI) * asigma) * Math.exp(-u*u/2);
}
+// Evaluate the log gaussian distribution
+function gaussian_lnpdf(x,mu,sigma)
+{
+ return -.5*(1.8378770664093453 + Math.log(sigma) + (x - mu)*(x - mu)/sigma);
+}
+
// Draw a sample from a Dirichlet distribution
// Law & Kelton (roughly optimal)
// TODO: may need to match function signature for Ikarus compatibility
@@ -198,12 +211,12 @@ function dirichlet_lnpdf(alpha, theta)
// Marsaglia '80
function sample_tdist(nu)
{
- if(nu <= 2) return sample_gaussian(1) / sqrt( 2 * sample_gamma(nu/2, 1) / nu);
+ if(nu <= 2) return sample_gaussian(0,1) / sqrt( 2 * sample_gamma(nu/2, 1) / nu);
var a,b,c,t;
do
{
- a = sample_gaussian(1);
+ a = sample_gaussian(0,1);
b = -1 / (nu/2 - 1) * log1p(-random());
c = a*a/(nu - 2);
}
View
8 scheme-compilers/scheme2js-template.sc
@@ -1,4 +1,4 @@
-;; Broken placeholder constants
+3;; Broken placeholder constants
(define infinity 999999999999)
(define minus-infinity (- 999999999999))
(define nan (/ 1 0))
@@ -18,12 +18,14 @@
;; These functions are defined in math-functions.js, but need aliases
(define sample-gamma sample_gamma)
-(define gamma-pdf gamma_pdf) ;; Why was this an ln version? Score computations are done in ln domain, so need ln-score...
+(define gamma-pdf gamma_pdf)
+(define gamma-lnpdf gamma_lnpdf)
(define sample-poisson sample_poisson)
(define sample-binomial sample_binomial)
(define sample-beta sample_beta)
(define sample-gaussian sample_gaussian)
-(define gaussian-pdf gaussian_pdf) ;; Why was this an ln version? Score computations are done in ln domain, so need ln-score...
+(define gaussian-pdf gaussian_pdf)
+(define gaussian-lnpdf gaussian_lnpdf)
(define sample-dirichlet sample_dirichlet)
(define dirichlet-lnpdf dirichlet_lnpdf)
(define sample-tdist sample_tdist)

0 comments on commit d578854

Please sign in to comment.