Permalink
Browse files

"Working" implementations of all the samplers and distributions excep…

…t for multinomial and discrete. These all need to be tested rigorously: there are undoubtedly bugs. Moreover, most of the samplers in some way depend on fast implementations of the Binomial or Gaussian samplers, and the implemented versions are inferior to the state of the art algorithms (but much simpler).

In addition, also added a set of replacements for Math.random().
  • Loading branch information...
1 parent 457ec87 commit 296f991101816e602bca879bcd79ee0444b5ebba @jerrytalton jerrytalton committed with longouyang Sep 22, 2011
View
18 bher
@@ -98,14 +98,14 @@ Examples:
for (i, (param_key, param_value)) in enumerate(params):
if param_key:
- vprint("* %s: %s" % (param_key, param_value), True)
+ vprint("* %s: %s" % (param_key, param_value), True)
settings["out_path"] = abspath(file) + "." + str(i) + ".ss"
else:
settings["out_path"] = abspath(file) + compilers[compiler]["extension"]
vprint("removing old compilation files (if existent)", verbose)
call("rm -f '%(out_path)s'" % settings, verbose, allow_fail=True)
- call("rm -f '%(out_path)s.tmp'" % settings, verbose, allow_fail=True)
+ call("rm -f '%(out_path)s.tmp'" % settings, verbose, allow_fail=True)
call("rm -f '%(in_path)s.tmp'" % settings, verbose, allow_fail=True)
vprint("processing input file", verbose)
@@ -123,7 +123,7 @@ Examples:
settings["churchprogram"] = open("%(out_path)s.tmp" % settings).read()
call("rm -f '%(in_path)s.tmp'" % settings, verbose, allow_fail=True)
call("rm -f '%(out_path)s.tmp'" % settings, verbose, allow_fail=True)
-
+
## vprint("adding compiler-specific header & footer", verbose)
## call(("cat '%(header_path)s' '%(out_path)s.tmp' '%(trailer_path)s' " +
## "> '%(out_path)s'") % settings, verbose)
@@ -138,27 +138,29 @@ Examples:
pre = datetime.now()
call("vicare %(debug)s --r6rs-script '%(out_path)s'" % settings, verbose)
post = datetime.now()
-
+
if time:
delta = post-pre
seconds = delta.seconds + delta.microseconds/1000000.0
print("runtime: %fs" % seconds)
-
+
if not keep:
vprint("removing compiled file", verbose)
call("rm -f '%(out_path)s'" % settings, verbose)
-
+
if compiler=="scheme2js":
vprint("compiling scheme to js via scheme2js ", verbose)
settings["js_out_path"] = settings["in_path"] + ".js"
- settings["js_math_path"] = abspath(os.path.join(settings["bher_path"], "scheme-compilers/math-functions.js"))
+ settings["js_mash_path"] = abspath(os.path.join(settings["bher_path"], "scheme-compilers/javascript/random/Mash.js"))
+ settings["js_random_path"] = abspath(os.path.join(settings["bher_path"], "scheme-compilers/javascript/random/MRG32k3a.js"))
+ settings["js_math_path"] = abspath(os.path.join(settings["bher_path"], "scheme-compilers/javascript/math-functions.js"))
settings["scheme2js_path"] = abspath(os.path.join(settings["bher_path"], "external/scheme2js/scheme2js.jar"))
settings["scheme2js_runtime_path"] = abspath(os.path.join(settings["bher_path"], "external/scheme2js/runtime.js"))
call("java -jar '%(scheme2js_path)s' '%(out_path)s' -o '%(js_out_path)s' -O1" % settings, verbose)
if not keep:
call("rm -f '%(out_path)s'" % settings, verbose)
vprint("constructing html file ", verbose)
- scheme2js_template = open("scheme-compilers/template-scheme2js.html").read()
+ scheme2js_template = open("scheme-compilers/template-scheme2js.html").read()
scheme2js_html = scheme2js_template % settings
f = open(settings["in_path"] + ".html", "w")
f.write(scheme2js_html)
@@ -1,3 +1,23 @@
+// Set the random number generator
+var random = new MRG32k3a();
+var intRandom = random.uint32;
+
+function random_integer(n)
+{
+ return intRandom() % n;
+}
+
+function random_real()
+{
+ return random();
+}
+
+function seed_rng(seed)
+{
+ random = new MRG32k3a(seed);
+ intRandom = random.uint32;
+}
+
// Draw sample from Poisson distribution
// Knuth TAOCP 2 (roughly optimal)
function sample_poisson(mu)
@@ -15,16 +35,22 @@ function sample_poisson(mu)
var emu = Math.exp(-mu);
var p = 1;
- do{ p *= Math.random(); k++; } while(p > emu);
+ do{ p *= random(); k++; } while(p > emu);
return k-1;
}
+// Poisson probability distribution function via iterative expansion
+function poisson_pdf(k, mu)
+{
+ return Math.exp(k * Math.log(mu) - mu - lnfact(k));
+}
+
// Draw sample from a Gamma distribution
// Marsagli and Tsang '00 (roughly optimal)
function sample_gamma(a,b)
{
- if(a < 1) return sample_gamma(1+a,b) * Math.pow(Math.random(), 1/a);
+ if(a < 1) return sample_gamma(1+a,b) * Math.pow(random(), 1/a);
var x,v,u;
var d = a-1/3;
@@ -35,12 +61,22 @@ function sample_gamma(a,b)
do{x = sample_gaussian(1); v = 1+c*x;} while(v <= 0);
v=v*v*v;
- u=Math.random();
+ u=random();
if((u < 1 - .331*x*x*x*x) || (Math.log(u) < .5*x*x + d*(1 - v + Math.log(v)))) return b*d*v;
}
}
+// Evaluate gamma pdf
+function gamma_pdf(x,a,b)
+{
+ if(x<0) return 0;
+ if(x==0) return a==1 ? 1/b : 0;
+ if(a==1) return Math.exp(-x/b)/b;
+
+ return Math.exp((a - 1)*Math.log(x/b) - x/b - log_gamma(a))/b;
+}
+
// Draw a sample from a Binomial distribution
// Knuth TAOCP 2 (could be improved, i.e. via Kachitvichyanukul & Schmeiser)
function sample_binomial(p,n)
@@ -63,13 +99,39 @@ function sample_binomial(p,n)
var u;
for(i=0; i<n; i++)
{
- u = Math.random();
+ u = random();
if(u<p) k++;
}
return k;
}
+// Binomial probability distribution function via Normal approximation
+// Peizer & Pratt 1968, JASA 63: 1416-1456 (may not be optimal...)
+function binomial_pdf(k, p, n)
+{
+ var inv2 = 1/2, inv3 = 1/3, inv6 = 1/6;
+
+ if (k >= n) return 1;
+
+ var q = 1 - p;
+ var s = k + inv2;
+ var t = n - k - inv2;
+ var d1 = s + inv6 - (n + inv3) * p;
+ var d2 = q /(s+inv2) - p/(t+inv2) + (q-inv2)/(n+1);
+
+ d2 = d1 + 0.02 * d2;
+
+ var num = 1 + q * g(s/(n*p)) + p * g(t/(n*q));
+ var den = (n + inv6) * p * q;
+ var z = num / den;
+
+ z = d2 * Math.sqrt(z);
+ z = normal_cdf(z);
+
+ return z;
+}
+
// Draw a sample from a Beta distribution
// Knuth TAOCP 2 (roughly optimal)
function sample_beta(a, b)
@@ -86,8 +148,8 @@ function sample_gaussian(sigma)
do
{
- u = 1 - Math.random();
- v = 1.7156 * (Math.random() - .5);
+ u = 1 - random();
+ v = 1.7156 * (random() - .5);
x = u - 0.449871;
y = Math.abs(v) + 0.386595;
q = x*x + y*(0.196*y - 0.25472*x);
@@ -97,6 +159,14 @@ function sample_gaussian(sigma)
return sigma*v/u;
}
+// Evaluate the gaussian distribution
+function gaussian_pdf(x,sigma)
+{
+ var asigma = Math.abs(sigma);
+ var u = x/asigma;
+ return (1/ Math.sqrt(2*Math.PI) * asigma) * Math.exp(-u*u/2);
+}
+
// Draw a sample from a Dirichlet distribution
// Law & Kelton (roughly optimal)
// TODO: may need to match function signature for Ikarus compatibility
@@ -112,6 +182,18 @@ function sample_dirichlet(alpha)
return theta;
}
+// Evaluate the logarithm of the Dirichlet distribution
+function dirichlet_lnpdf(alpha, theta)
+{
+ var logp = 0;
+
+ for(i=0; i<alpha.length; i++) logp = (alpha[i] - 1)*Math.log(theta[i]);
+ logp += log_gamma(sum(alpha));
+ for(i=0; i<alpha.length; i++) logp -= log_gamma(alpha[i]);
+
+ return logp;
+}
+
// Draw a sample from a Student's t-distribution
// Marsaglia '80
function sample_tdist(nu)
@@ -122,14 +204,84 @@ function sample_tdist(nu)
do
{
a = sample_gaussian(1);
- b = -1 / (nu/2 - 1) * log1p(-Math.random());
+ b = -1 / (nu/2 - 1) * log1p(-random());
c = a*a/(nu - 2);
}
while(1-c < 0 || Math.exp(-b-c) > (1-c));
return a / Math.sqrt((1-c/nu) * (1-c));
}
+// Evaluate t-distribution
+function tdist_pdf(x,nu)
+{
+ var a = log_gamma(nu/2);
+ var b = log_gamma((nu+1)/2);
+
+ return Math.exp(b-a)/Math.sqrt(Math.PI*nu) * Math.pow(1 + x*x/nu, -(nu+1)/2);
+}
+
+// Draw a sample from a generalized t-distribution
+function sample_generalized_tdist(nu,mu,sigma_squared)
+{
+ return sample_tdist(nu)*Math.sqrt(sigma_squared) + mu;
+}
+
+// Return the log of a sum of exponentials, to minimize under/overflow
+function logsumexp(v)
+{
+ var t=0;
+ var v;
+
+ for(i=0;i<v.length;i++)
+ {
+ abs=Math.abs(v[i]);
+ if(abs>t){ t=abs; val=v[i]; }
+ }
+
+ var sum=0;
+ for(i=0;i<v.length;i++) sum += Math.exp(v[i]-val);
+
+ return Math.log(sum) + v;
+}
+
+// Evaluate the log of gamma(x)
+// Lancsoz approximation from Numerical Recipes in C
+function log_gamma(xx)
+{
+ var cof = [76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5];
+
+ var x = xx - 1.0;
+ var tmp = x + 5.5; tmp -= (x + 0.5)*Math.log(tmp);
+ var ser=1.000000000190015;
+ for (j=0;j<=5;j++){ x++; ser += cof[j]/x; }
+ return -tmp+Math.log(2.5066282746310005*ser);
+}
+
+// Calculate the sum of elements in a vector
+function sum(v)
+{
+ var sum=0;
+ for(i=0;i<v.length;i++) sum += v[i];
+ return sum;
+}
+
+// Calculate the mean of elements in a vector
+function mean(v)
+{
+ return sum(v)/v.length;
+}
+
+// Normalize a vector
+function normalize(v)
+{
+ var s=0;
+ for(i=0;i<v.length;i++) s += v[i]*v[i];
+ s = Math.sqrt(s);
+ for(i=0;i<v.length;i++) v[i] /= s;
+ return v;
+}
+
// Returns log(1 + x) in a numerically stable way
function log1p(x)
{
@@ -228,35 +380,3 @@ function g(x)
}
return z;
}
-
-// Binomial probability distribution function via Normal approximation
-// Peizer & Pratt 1968, JASA 63: 1416-1456
-function binomial_pdf(k, p, n)
-{
- var inv2 = 1/2, inv3 = 1/3, inv6 = 1/6;
-
- if (k >= n) return 1;
-
- var q = 1 - p;
- var s = k + inv2;
- var t = n - k - inv2;
- var d1 = s + inv6 - (n + inv3) * p;
- var d2 = q /(s+inv2) - p/(t+inv2) + (q-inv2)/(n+1);
-
- d2 = d1 + 0.02 * d2;
-
- var num = 1 + q * g(s/(n*p)) + p * g(t/(n*q));
- var den = (n + inv6) * p * q;
- var z = num / den;
-
- z = d2 * Math.sqrt(z);
- z = normal_cdf(z);
-
- return z;
-}
-
-// Poisson probability distribution function via iterative expansion
-function poisson_pdf(k, mu)
-{
- return Math.exp(k * Math.log(mu) - mu - lnfact(k));
-}
@@ -0,0 +1,52 @@
+// From http://baagoe.com/en/RandomMusings/javascript/
+function Alea() {
+ return (function(args) {
+ // Johannes Baagøe <baagoe@baagoe.com>, 2010
+ var s0 = 0;
+ var s1 = 0;
+ var s2 = 0;
+ var c = 1;
+
+ if (args.length == 0) {
+ args = [+new Date];
+ }
+ var mash = Mash();M
+ s0 = mash(' ');
+ s1 = mash(' ');
+ s2 = mash(' ');
+
+ for (var i = 0; i < args.length; i++) {
+ s0 -= mash(args[i]);
+ if (s0 < 0) {
+ s0 += 1;
+ }
+ s1 -= mash(args[i]);
+ if (s1 < 0) {
+ s1 += 1;
+ }
+ s2 -= mash(args[i]);
+ if (s2 < 0) {
+ s2 += 1;
+ }
+ }
+ mash = null;
+
+ var random = function() {
+ var t = 2091639 * s0 + c * 2.3283064365386963e-10; // 2^-32
+ s0 = s1;
+ s1 = s2;
+ return s2 = t - (c = t | 0);
+ };
+ random.uint32 = function() {
+ return random() * 0x100000000; // 2^32
+ };
+ random.fract53 = function() {
+ return random() +
+ (random() * 0x200000 | 0) * 1.1102230246251565e-16; // 2^-53
+ };
+ random.version = 'Alea 0.9';
+ random.args = args;
+ return random;
+
+ } (Array.prototype.slice.call(arguments)));
+};
Oops, something went wrong.

0 comments on commit 296f991

Please sign in to comment.