search

JohannesBuchner / PyMultiNest

Pythonic Bayesian inference and visualization for the MultiNest Nested Sampling Algorithm and PyCuba's cubature algorithms.
http://johannesbuchner.github.io/PyMultiNest/
Other
194 stars 88 forks source link

Problem converging at the edge of the prior #252

Closed sumitaghosh closed 5 months ago

sumitaghosh commented 5 months ago

When I try to run on an asimov dataset with one parameter set to 0, PyMultiNest never seems to consider 0 as an option.

This code:

from pymultinest.solve import Solver
import numpy as np
from scipy.special import ndtri

class EdgeTest(Solver):
    def __init__(self, **kwargs):
        # now add my stuff
        self.prior_dictionary = {
            'help': ('type', 'mean', ('low sigma', 'high sigma'), ('low limit', 'high limit')),
            "Rate_Xe136_0nu_XeLS": ("uniform", 0.00, (0, 0), (0, 100)),
            "Rate_Xe136_2nu_XeLS": ("uniform", 107743, (0, 0), (0, 5000000)),
            "Rate_U238_S2_XeLS": ("gaussian", 31.119518, (1.599997, 1.599997), (0, 1000))
        }
        self.keys = ['Rate_Xe136_0nu_XeLS',
                     'Rate_Xe136_2nu_XeLS',
                     'Rate_U238_S2_XeLS']
        self.energies = [1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.8, 1.85, 1.9, 1.95, 2.0, 2.05, 2.1, 2.15, 2.2, 2.25, 2.3,
                         2.35, 2.4, 2.45, 2.5, 2.55, 2.6, 2.65, 2.7, 2.75, 2.8, 2.85, 2.9, 2.95, 3.0]
        self.num_bins = len(self.energies) - 1
        self.singles_pdfs = [
            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0202388627378247e-05, 9.322164424936633e-05,
                      0.0007690453438840364, 0.00540569002010639, 0.028218190513198682, 0.12680151725234257,
                      0.42734238135610686, 1.1160401165758493, 2.3096656576084875, 3.7430419096654917,
                      4.775387753471517, 4.877187329942036, 3.9618515568174786, 2.594490026200128, 1.3725374362439022,
                      0.5930706356033242, 0.20785100505456394, 0.06497153968436205, 0.016736819363366334,
                      0.004067050837079363, 0.0012914943198711777, 0.00026701527555763585]),
            np.array([0.391718430066749, 0.3296560935827085, 0.27291943475258806, 0.2213611740895031,
                      0.17677717566964501, 0.1382513575889907, 0.10423149762768026, 0.07747041685101658,
                      0.05615467371189704, 0.03922946899057661, 0.026363843916357824, 0.0166612669007192,
                      0.010197667428678139, 0.005879412520757319, 0.003191151426513824, 0.001603981554648969,
                      0.000742813576244242, 0.0003184041485716475, 0.00012381364974478314, 4.2718660269486025e-05,
                      1.2725186619549708e-05, 3.4520783175797974e-06, 8.534196569655529e-07, 2.6808887101253576e-07,
                      3.891115771731739e-08, 0.0, 0.0, 0.0, 0.0, 0.0]),
            np.array([0.00013970272910318227, 0.0001695444980874317, 0.0002130338691449695, 0.00026077700622307936,
                      0.00030592190753134046, 0.00034896740739312113, 0.0003776835718418515, 0.00040675468222436516,
                      0.00042855590716811315, 0.0004481596968498875, 0.0004665415263580842, 0.00047659673461909356,
                      0.00048798526518804245, 0.0004914968855994857, 0.0004953852533590487, 0.0004920550234329798,
                      0.0004830504806493089, 0.0004765854698489232, 0.00045760723400835746, 0.000430421700289296,
                      0.00040365107203471553, 0.00036916303155576987, 0.0003315153570877531, 0.0002973201488910318,
                      0.0002569909587913854, 0.00021615705815814851, 0.00017729908906460946, 0.00014155612982320797,
                      0.0001068346093281903, 7.87202087528126e-05])
        ]
        self.ll_pdfs = [
            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.3248256513969964e-06, 7.998673213754771e-06,
                      6.917575811100022e-05, 0.0004765254709399369, 0.0024969538874807544, 0.01122931768685254,
                      0.03783047451080829, 0.09897734061932392, 0.20520961308004582, 0.3330744525409708,
                      0.425322617169032, 0.4343943102006812, 0.35304251455284624, 0.23107437355200136,
                      0.12209713096825664, 0.052632003341252245, 0.018346841724733803, 0.005670245939716155,
                      0.0014213919597989948, 0.0003370737792653955, 0.0001024308743483086, 2.0021083743262294e-05]),
            np.array([0.034794924836000427, 0.029275462285747123, 0.024249784015954644, 0.019677734478505283,
                      0.015710674797310885, 0.012281235026670882, 0.009256651482260273, 0.006884535594277013,
                      0.004997695587889409, 0.0034854488486913554, 0.0023421154378630736, 0.001478700892105994,
                      0.0009054456368268835, 0.0005218672440850021, 0.00028295070270060617, 0.0001422432252774825,
                      6.561424443550626e-05, 2.81740685973128e-05, 1.0925448010457665e-05, 3.745235473707396e-06,
                      1.1205462636825924e-06, 3.0395082023730854e-07, 7.375334367608048e-08, 2.0968791718140986e-08,
                      2.9123653288131696e-09, 0.0, 0.0, 0.0, 0.0, 0.0]),
            np.array([1.2465657214400514e-05, 1.5183311162005786e-05, 1.9100606198328514e-05, 2.3331568170041787e-05,
                      2.7357261703310833e-05, 3.110906227836709e-05, 3.36310823265831e-05, 3.620776257794921e-05,
                      3.810878115396978e-05, 3.977754016072002e-05, 4.138710438701383e-05, 4.2330490970748957e-05,
                      4.331792712311154e-05, 4.361558488588878e-05, 4.3981707457409194e-05, 4.362484220507875e-05,
                      4.2877048216007714e-05, 4.221149735133398e-05, 4.053596414657667e-05, 3.809229630022501e-05,
                      3.563567037608486e-05, 3.258457306428801e-05, 2.9241409530376086e-05, 2.616571944712311e-05,
                      2.2566739909996788e-05, 1.893405868338155e-05, 1.550385405335905e-05, 1.233331294310511e-05,
                      9.2731532407586e-06, 6.795393018000643e-06])
        ]
        self.true_values = [0, 107743, 31.1195180]
        self.singles_data = np.zeros(self.num_bins)
        self.long_lived_data = np.zeros(self.num_bins)
        for t, s_pdf, l_pdf in zip(self.true_values, self.singles_pdfs, self.ll_pdfs):
            self.singles_data += t * s_pdf
            self.long_lived_data += t * l_pdf
        super().__init__(**kwargs)

    @staticmethod
    def return_gaussian(uniform, mu, low_sigma, high_sigma):
        assert low_sigma != 0 or high_sigma != 0  # at least one of them should be nonzero
        if uniform > 0.5:  # ndtri(uniform) is positive between 0 and 1, so this will be positive
            return high_sigma * ndtri(uniform) + mu
        else:  # ndtri(uniform) would be negative (between -1 and 0)
            theta = low_sigma * ndtri(uniform) + mu
        return max(theta, 0.)  # ensure we're positive semi-definite

    def Prior(self, uniform):
        theta = np.zeros(len(uniform))
        for i, k in enumerate(self.keys):
            distribution_type, mu, sigmas, limits = self.prior_dictionary[k]
            assert distribution_type in ('gaussian', 'uniform', 'fixed')
            if distribution_type == 'gaussian':
                theta[i] = EdgeTest.return_gaussian(uniform[i], mu, sigmas[0], sigmas[1])
            elif distribution_type == 'uniform':
                theta[i] = uniform[i] * limits[1]
            elif distribution_type == 'fixed':
                theta[i] = mu
        assert np.all(np.isfinite(theta))
        assert np.all(theta >= 0)  # ensure theta is non-negative
        return theta

    def LogLikelihood(self, theta):
        arrival_rate_shape = self.singles_data.shape
        singles_arrival = np.zeros(arrival_rate_shape)
        long_lived_rate = np.zeros(arrival_rate_shape)
        for par, s_pdf, l_pdf in zip(theta, self.singles_pdfs, self.ll_pdfs):
            singles_arrival += par * s_pdf
            long_lived_rate += par * l_pdf
        to_sum = ((self.singles_data * np.log(singles_arrival, where=singles_arrival > 0) - singles_arrival) +
                  (self.long_lived_data * np.log(long_lived_rate, where=long_lived_rate > 0) - long_lived_rate))
        return np.sum(to_sum)

solution = EdgeTest(n_dims=3)

print(solution)
print('Recall that the true values are', solution.true_values)

results in this message:

 MultiNest Warning: no resume file found, starting from scratch
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  400
 dimensionality =    3
 *****************************************************

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            1  is converging towards the edge of the prior.
 Parameter            2  of mode            1  is converging towards the edge of the prior.
 ln(ev)=   1930588.3089250897      +/-  0.18600451493200451     
 Total Likelihood Evaluations:        10309
 Sampling finished. Exiting MultiNest
  analysing data from /var/folders/ck/05t1k3v16z916q_jc7ht4qtc0000gn/T/tmph_u590z6pymultinest/.txt
Model in "(temporary directory)" (3 dimensions)
Evidence ln Z = 1930588.3 +- 0.2
Parameter values:
   Parameter  1 : 0.148 +- 0.125
   Parameter  2 : 107738.549 +- 231.254
   Parameter  3 : 31.036 +- 1.504
Recall that the true values are [0, 107743, 31.119518]

Process finished with exit code 0
JohannesBuchner commented 5 months ago

Closing for now, not enough information provided. Please give a minimal code example that reproduces the issue.

sumitaghosh commented 5 months ago

I have done so! Please reopen this issue @JohannesBuchner and thank you for the instructions!

JohannesBuchner commented 5 months ago

Can you plot the posterior probability distribution? The results look like expected to me: Your prior on the first variable is forced to be non-negative. Keep in mind when working with probability distributions, the probability of any specific value (such as zero) is zero.

sumitaghosh commented 5 months ago

testmarg

I expected the error to be equal to the mode so that the value was consistent with zero, so thanks for explaining why that's not correct!

sumitaghosh commented 4 months ago

I have an unrelated issue I haven't been able to solve with google. I ran another asimov dataset with more parameters, and it seems to be stuck in a loop (running 13+ hours so far). I think the problem is that the evidence keeps registering as something really large, so the following is repeating over and over:

Importance Nested Sampling ln(Z):            NaN +/-       NaN
Acceptance Rate:                        0.497320
Replacements:                             869800
Total Samples:                           1748975

I think this is happening because my second parameter has gotten stuck on the value -0.437431652236425247E+11 (when I check the .txt file, this is what I see in the second column), but I don't know how this is possible when my prior defines this parameter as theta[1] = uniform[1] * 5000000.

Do you know why this could be happening, and what I can do to fix it? What would cause a number that can't go below 0 or above 5e6 become -4e10? Thank you for your help!

JohannesBuchner commented 4 months ago

maybe put some prints in your likelihood to see what it is doing. Maybe some asserts too.