Merge pull request #4 from mattpitkin/inverse_transform_sampling

Convert to using inverse transform sampling from unit hypercube to prior space

src/nested_sampler.m

@@ -151,24 +151,8 @@
     extraparvals = [];
 end
 
-% draw the set of initial live points from the prior
-livepoints = zeros(Nlive, D);
-
-for i=1:D
-    priortype = prior{i,2};
-    p3 = prior{i,3};
-    p4 = prior{i,4};
-
-    % currently only handles uniform or Gaussian priors
-    if strcmp(priortype, 'uniform')
-        livepoints(:,i) = p3 + (p4-p3)*rand(Nlive,1);
-    elseif strcmp(priortype, 'gaussian')
-        livepoints(:,i) = p3 + p4*randn(Nlive,1);
-    elseif strcmp(priortype, 'jeffreys')
-        % uniform in log space
-        livepoints(:,i) = 10.^(log10(p3) + (log10(p4)-log10(p3))*rand(Nlive,1));
-    end
-end
+% draw the set of initial live points from the unit hypercube
+livepoints = rand(Nlive, D);
 
 % check whether likelihood is a function handle, or a string that is a
 % function name
@@ -184,16 +168,11 @@
 logL = zeros(Nlive,1);
 
 for i=1:Nlive
-    parvals = cat(1, num2cell(livepoints(i,:)'), extraparvals);
+    % rescale points for calculating likelihood
+    parvals = cat(1, num2cell(rescale_parameters(prior, livepoints(i,:))), extraparvals);
     logL(i) = feval(flike, data, model, parnames, parvals);
 end
 
-% now scale the parameters, so that uniform parameters range from 0->1,
-% and Gaussian parameters have a mean of zero and unit standard deviation
-for i=1:Nlive
-    livepoints(i,:) = scale_parameters(prior, livepoints(i,:));
-end
-
 % initial tolerance
 tol = inf;
 

src/rescale_parameters.m

@@ -2,7 +2,9 @@
 
 % scaled = rescale_parameters(prior, params)
 %
-% This function will do the reverse of scale_parameters.
+% This function will convert from a unit hypercube to the true
+% parameter values. For the Gaussian prior inverse transform
+% sampling is used to convert from the unit hypercube to a Gaussian.
 
 lp = length(params);
 
@@ -17,7 +19,8 @@
     if strcmp(priortype, 'uniform')
         scaled(i) = params(i)*(p4 - p3) + p3;
     elseif strcmp(priortype, 'gaussian')
-        scaled(i) = params(i)*p4 + p3;
+        % use inverse transform sampling
+        scaled(i) = p3 + p4*sqrt(2)*erfinv(2.*params(i) - 1.);
     elseif strcmp(priortype, 'jeffreys')
         scaled(i) = 10^(params(i)*(log10(p4) - log10(p3)) + log10(p3));
     end