Issue with priors

[DavidMoiseNataf]

Hi Tim,

I get a crash when I try and use a power-law prior with alpha=-1, e.g. model1.set_prior(AV=PowerLawPrior(alpha=-1., bounds=(0.01,BoundExtinction)))

I assume that this is because the integral of 1/x is a logarithm rather than a power-law, and so the definition of "C" in the priors.py needs to have an "if" statement for those special cases (see below). Can you make the fix in a consistent manner? Thank you.

def __init__(self, alpha, bounds=None):
    self.alpha = alpha
    super(PowerLawPrior, self).__init__(bounds=bounds)

def _pdf(self, x):
    lo, hi = self.bounds
    C = (1 + self.alpha) / (hi ** (1 + self.alpha) - lo ** (1 + self.alpha))
    # C = 1/(1/(self.alpha+1)*(1 - lo**(self.alpha+1)))
    return C * x ** self.alpha

def _lnpdf(self, x):
    lo, hi = self.bounds
    C = (1 + self.alpha) / (hi ** (1 + self.alpha) - lo ** (1 + self.alpha))
    return np.log(C) + self.alpha * np.log(x)

def sample(self, n):
    """

    cdf = C * int(x**a)|x=lo..x

        = C * [x**(a+1) / (a+1)] | x=lo..x
        = C * ([x**(a+1) / (a+1)] - [lo**(a+1) / (a+1)])
     u  =

     u/C + [lo**(a+1) / (a+1)] = x**(a+1) / (a+1)
     (a+1) * (u/C + [lo**(a+1) / (a+1)]) = x**(a+1)
     [(a+1) * (u/C + [lo**(a+1) / (a+1)])] ** (1/(a+1)) = x
    """
    lo, hi = self.bounds
    C = (1 + self.alpha) / (hi ** (1 + self.alpha) - lo ** (1 + self.alpha))
    u = np.random.random(n)
    a = self.alpha
    return ((a + 1) * (u / C + (lo ** (a + 1) / (a + 1)))) ** (1 / (a + 1))

[timothydmorton]

Yes, you are correct-- please feel free to make the fix yourself and submit a PR! Another quick-fix solution would be just to set alpha=-1.00001...