Gaussian init had problems

Often the range was outside of what cerf can handle causing NaNs.
Using the gsl module fixes this.

src/fitting/ensemble-samplers/emcee.sl

@@ -3,6 +3,7 @@
 require("rand");
 require("fork");
 require("socket");
+require("gsl", "gsl");
 %require("select");
 
 % Implementation of the emcee hammer (Foreman-Mackey 2013) with
@@ -1025,10 +1026,14 @@ EmceeInitRegister["uniform"] =  &emceeInitUniform;
 %{{{ Gauss initialization function
 private define rand_gauss_cut (sigma, v, bmin, bmax) %{{{
 {
-  variable upper = Real(cerf((bmax-v)/sqrt(2.)/sigma));
-  variable lower = Real(cerf((bmin-v)/sqrt(2.)/sigma));
+  variable mu = (v-bmin)/(bmax-bmin);
+  variable alpha = -mu/sigma/sqrt(2);
+  variable beta = (1.0-mu)/sigma/sqrt(2);
+  variable phia = 0.5*(1+gsl->erf(alpha));
+  variable phib = 0.5*(1+gsl->erf(beta));
+  variable k = phia+rand_uniform(length(v))*(phib-phia);
 
-  return _min(bmax,_max(bmin,sqrt(2)*erfinv(rand_uniform(length(v))*(upper-lower)+lower)*sigma+v));
+  return _min(bmax,_max(bmin,gsl->cdf_gaussian_Pinv(k, sigma)*sigma*(bmax-bmin)+v));
 }
 %}}}
 
@@ -1039,7 +1044,6 @@ private define emceeInitGaussPick (init, engine) %{{{
   variable file = engine.leader.inFile;
   variable par = __parameters(engine.fit.object);
   variable numParameter = length(par.value);
-  variable sigma;
 
   if (file.mode & EMCEE_FILE_RANGE)
     (feed, ) = file.read(engine, engine.totalNumberWalkers);
@@ -1051,8 +1055,7 @@ private define emceeInitGaussPick (init, engine) %{{{
   k = 0;
   _for i (0, length(c)-1) {
     _for j (0, c[i]-1) {
-      sigma = (p[i].max-p[i].min)*init.sigma;
-      engine.walkers[j+k] = rand_gauss_cut(sigma, p[i].value, p[i].min, p[i].max);
+      engine.walkers[j+k] = rand_gauss_cut(init.sigma, p[i].value, p[i].min, p[i].max);
     }
     k += c[i];
   }