Matlab MultiNest

10 Replies to “Matlab MultiNest”

1. was wondering if the author had a plan to implement the extra features in the fortran 90 version.

   Thanks!

   Yun

   Reply

2. Hi Yun,

   I don’t have any plan to implement the extra features in the main fortran version of MultiNest. I’ve only done minimal amounts of coding in fortran and don’t have the time to add this.

   Cheers,

   Matt

   Reply

3. Thanks for the prompt answer Matt!

   Best regards,

   Yun

   Reply

4. Hi,

   I just downloaded your MATLAB implementation of multinest from CCForge. Running the examples with a non-zero MCMC steps parameter (that the examples request when first running it) seems to work all right, but if I enter zero (indicating that the program should use the sampling scheme from multinest), it fails with the errors:

   Input number of iterations for MCMC sampling: (enter 0 for multinest sampling)
   0
   ??? Error using ==> mldivide
   Matrix dimensions must agree.

   Error in ==> calc_ellipsoid at 56
   invC = 1/C;

   Error in ==> optimal_ellipsoids at 27
   [B, invB, mu, VE, flag] = calc_ellipsoid(u, VS);

   Error in ==> nested_sampler at 205
   [Bs, mus, VEs, ns] = optimal_ellipsoids(livepoints, VS);

   Error in ==> example_line_with_outliers_and_marginalization at 29
   [logZ, nest_samples, post_samples] = nested_sampler(data, Nlive, Nmcmc, …

   >>

   —–

   I added all three of the directories to my path so that all of the .m files can be seen by each script before I ran them.

   If I go into the source code to this line in calc_ellipsoid, I uncomment the line that uses inv(C) and comment the line that uses 1/C:

   % calculate inverse covariance matrix
   invC = inv(C);
   %invC = 1/C;

   —

   Then the program appears to run OK using the setting of 0 for MCMC steps.

   Just FYI and hope any comments you can provide will be seen by others; I think actually MATLAB encourages use of 1/C for the inverse instead of inv, as you’ve done, but in this case it seems to fail for some reason.

   David

   Reply

   1. Hi David,

      Sorry for the slow response. Thanks very much for spotting this and contacting me. I’d not seen the problem before, but from the looks of it I think I changed from using inv(C) to 1/C at Matlab’s suggestion without properly testing that it works! Retesting with example_sinusiod.m and I see the same error as you.

      Looking at this again I suggest changing the calc_ellipsiod.m function to not output invB at all, removing the variable invC from the code, changing the line:
      f = (u(i,:)-mu) * invC * (u(i,:)-mu)'; to f = ( (u(i,:)-mu) / C ) * (u(i,:)-mu)';
      and removing the line defining invB.

      This also requires editing the optimal_ellipsoids.m function to remove invB from the call to calc_ellipsiod (it’s not actually used elsewhere in the function), and editing the split_ellipsiod.m function to again remove invB from the two calls to calc_ellipsiod and then changing lines defining du1 and du2 to be:
      du1 = ( (u(i,:)-mu1) / B1 ) * (u(i,:)-mu1)';
      and
      du2 = ( (u(i,:)-mu2) / B2 ) * (u(i,:)-mu2)';

      I’ll update the version on CCForge.

      Cheers,

      Matt

      Reply

5. Hi,
   I’ve downloaded the Matlab version of multinest, and I’m using it to analyse some neutron reflectivity data. I got it running very easily and it seems to run very smoothly, so thanks very much for this code!
   I’m essentially using it to analyse some 1d (i.e. x-y) data fitted with a multiparameter model, and I’m getting very sensible posterior distributions out of it. What I would like to do though is to then plot the resulting fits with shaded 95% (or any other) probability regions shown around the best fit curves. Looking at the main multinest page, the package ‘superplot’ looks like it would do this, but that’s for the Python implementation, so I’m a bit stuck….
   Have you got any functions that would calculate these regions, or at least suggest how I can proceed? So for example its says that Superplot will calculate these from “….Multinest output files”. If I could work out what the format of those is, then maybe I could format the Matlab output in a way that I can then feed through the Superplot calculation routines or something…. Anyway, any ideas as to what I might try would be most welcome…
   Many thanks,
   Arwel

   Reply

   1. Hi Arwel,

      Thanks for your comments. Sorry for the extremely slow response. I hope you found out how to do what you were wanting elsewhere, but below’s an example of how I would go about calculating and plotting a confidence region (or credible interval). Note that this example produces a region that is the shortest region to bound the given interval and as such is not necessarily symmetric about the best fit value.

      % make some posterior samples (instead of producing them with the nested sampling routine)
      post = randn(10000, 1);
      
      psamples = sort(post(:,1)); % sort posterior values in ascending order
      
      lowbound = psamples(1);         % initial low bound
      highbound = psamples(end);      % initial high bound
      cregion = highbound - lowbound; % initial region
      
      ci = 0.95; % required fractional confidence interval (e.g. 95% = 0.95)
      
      nsamp = length(psamples); % number of samples
      cllen = floor(nsamp*ci);
      
      % find shortest interval spanning required confidence interval
      for j=1:(nsamp-cllen)
        if psamples(j+cllen)-psamples(j) < cregion
          lowbound = psamples(j);
          highbound = psamples(j+cllen);
          cregion = highbound - lowbound;
        end
      end
      
      % histogram the data
      nbins = 50; % number of histogram bins
      [n, bins] = hist(psamples, nbins);
      
      % bin width
      dx = bins(2)-bins(1)
      
      % normalise histogram
      n = n/(sum(n)*dx);
      
      % plot histogrammed data
      plot(bins, n, 'b');
      
      % create polygon of confidence region to fill
      i1 = ceil((lowbound-bins(1))/dx); % get index of nearest bin to lower confidence bound
      y1 = n(i1) + (lowbound-bins(i1))*(n(i1+1)-n(i1))/dx; % linearly interpolate between nearest bins
      i2 = ceil((highbound-bins(1))/dx); % get index of nearest bin to upper confidence bound
      y2 = n(i2) + (highbound-bins(i2))*(n(i2+1)-n(i2))/dx;
      X = [lowbound, lowbound, bins(i1+1:i2), highbound, highbound];
      Y = [0., y1, n(i1+1:i2), y2, 0.];
      
      hold on;
      fill(X, Y, 'k', 'facealpha', 0.2, 'edgecolor', 'b'); % fill in region with transparent black
      hold off;
      

      That method gets the interval directly from the samples, but you could get it from the cumulative probability distribution (from which you could get an interval symmetric about the best fit).

      Reply

6. Ah thanks Matthew, that’s great.
   Id does sound as though it’s the cumulative probability distribution I need. How can I get the cumulative distribution from the output of multinest??
   Cheers,
   Arwel

   Reply

   1. Here’s a method for getting a confidence interval from the cumulative probability distribution. Given posterior samples from the nested sampling code in the array post you can get the cumulative distribution for the first parameter with e.g.:

       % get posterior samples for first parameter
      postvals = post(:,1);
      
      % histogram the data to get the distribution
      [n, bins] = hist(postvals, length(postvals));
      
      % get the cumulative probability distribution as normalise it
      cs = cumsum(n)/length(postvals);
      
      % have a look at the cdf! (which should go between 0 and 1)
      plot(bins, cs);
      
      % number of standard deviations you want for your confidence interval
      Nsigma = 1;
      
      % the probability covered by Nsigma (or you could just set 
      % e.g. probarea = 0.95, which is close to the 2sigma value)
      probarea = erf(Nsigma/sqrt(2));
      
      % get unique values of the cdf (otherwise the interpolation function
      % complains)
      [csu, iu, ju] = unique(cs);
      
      % get lower bound on the interval
      lowrange = interp1(csu, bins(iu), 0.5-(probarea/2));
      
       % get upper bound on the interval
      highrange = interp1(csu, bins(iu), 0.5+(probarea/2));
      

      That method has the probability interval centred around the point where 50% of the probability is on either side (the 0.5 in lowrange and highrange).
      Hope that helps.

      Reply

7. Hi, I tried the Matlab version of MultiNest. Is there a reason why the egg-box example runs forever and do not finally converge? Also is there a way to calculate the error in the evidence as in the original MultiNest paper? Also how do you calculate number of likelihood evaluations? The number of iterations (displayed in the workspace) and number of rows of the posterior samples are different. Which one is the number of likelihood calls?

   Reply