Chabrier 1/m term

[segasai]

Hi,

Are you sure that the current Chabrier implementation is correct ? I think it is currently lacking the 1/m term. As Chabrier distribution is strictly a log-normal distribution. I.e dn/dm ~ 1/m exp(-(log(m/m0))^2/s^2), or dn/dlogm ~ exp(-(log(m/m0)^2/s^2) Your current implementation is supposed to return dn/dm, but according to your lognormal() definition doesn't have the 1/m term.

Thanks, Sergey

[keflavich]

Good question. I haven't used the Chabrier IMF in production at all, so it may have errors, and there are no correctness tests yet. It would be very helpful to add those tests.

I think you're right, though, the Chabrier IMF function looks like it's returning dn/dlogm not dn/dm. A fix and/or tests would be welcome as a PR

[segasai]

Okay. I'll see if I have time. I am also not quite sure where the width 0.57 for a simple Chabrier comes from. The problem also is that there are many Chabrier MFs even in 2003 paper. Also, in your code, do you ever care about the normalization constant for dn/dm ? Because if you don't, for the lognormal at least it is easier to just plugin scipy.stats.lognorm().pdf() (that'd be good for all other MFs as well, as you can get rid of all the maths, and just use scipy.stats.powerlaw.pdf, scipy.stats.powerlaw.cdf etc. )

[keflavich]

well, the Chabrier MF that I want to use (which I think is why I refer to their 2005 paper... have to look again) is the lognormal + power-law. I've mostly been using the mass functions here to study the high-mass end of the IMF, which means a pure lognormal is generally not valid.

Yes, I care about the normalization. It would be great to make sure I got it right... I think the tests cover normalization for Kroupa, but I don't remember.

Also... not sure scipy.stats.powerlaw is any good for broken power laws?

[segasai]

I see. I am mostly interested in low mass stuff (and I don't care about normalization). I take back my .powerlaw() comment, as I see it is only powerlaw on [0,1] domain.

[segasai]

I guess what I meant is using normalized distributions everywhere in the code makes things way easier in terms of math. Here is just a toy implementation of a broken power law with arbitrary number of segments and trivial maths

import numpy as np
class LogNorm:
    def __init__(self, mu, s):
        self.d = scipy.stats.lognorm(s=s, scale=mu)

    def pdf(self):
        return self.d.pdf(x)

    def cdf(self, x):
        return self.d.cdf(x)

class PowerL:
    def __init__(self, slope, m1, m2):
        self.slope = slope
        self.m1 = m1
        self.m2 = m2

    def pdf(self, x):
        return x**self.slope * (self.slope + 1) / (
            self.m2**(self.slope + 1) -
            self.m1**(self.slope + 1)) * (x >= self.m1) * (x <= self.m2)

    def cdf(self, x):
        return (np.clip(x, self.m1, self.m2)**(self.slope + 1) -
                (self.m1**(self.slope + 1))) / (self.m2**(self.slope + 1) -
                                                self.m1**(self.slope + 1))

class BrokenPowerL:
    def __init__(self, slopes, breaks):
        assert (len(slopes) == len(breaks) - 1)
        nsegm = len(slopes)
        pows = []
        for i in range(nsegm):
            pows.append(PowerL(slopes[i], breaks[i], breaks[i + 1]))
        weights = [1]
        for i in range(1, nsegm):
            rat = pows[i].pdf(breaks[i]) / pows[i - 1].pdf(breaks[i])
            weights.append(weights[-1] / rat)
        weights = np.array(weights)
        self.slopes = slopes
        self.breaks = breaks
        self.pows = pows
        self.weights = weights / np.sum(weights)
        self.nsegm = nsegm

    def pdf(self, x):
        x1 = np.asarray(x)
        ret = x1 * 0.
        for i in range(self.nsegm):
            xind = (x1 < self.breaks[i + 1]) & (x1 > self.breaks[i])
            if xind.sum()>0:
                ret[xind] = self.weights[i] * self.pows[i].pdf(x1[xind])
        return ret
###
xgrid = np.linspace(0.10001,3.999,100000)
plt.plot(xgrid,BrokenPowerL([-1.1,-2,-3],[0.1,1,2,4]).pdf(xgrid))

These distributions all integrate to 1 by default. If you care about absolute normalization, that's just one extra multiplication. And you can use the same principal to combine the log-normal and power-law -- it is trivial.

[keflavich]

Yeah, I think that's a nicer version of the general case power law. I think you'd need another special case for an upper mass limit, right? You could just assign the last segment zero weight...

Is your definition any different from BrokenPowerLaw in imf.py as is now? I'm not sure, haven't dug through in enough detail

[segasai]

For the upper mass limit, putting infinity should work fine. The main difference to the existing code is the integration is much more trivial, and all the maths are locked inside the PowerLaw function() In fact I extended the code a bit, so now you can combine an arbitrary number of distributions in intervals. https://gist.github.com/segasai/2a5bbd6f5bedbac9e7d2c841d0796de9 This basically makes Chabrier2003, Chabrier2005, Kroupa -- one-liners.

[keflavich]

👍 want to submit these changes/additions as a PR? It might be good to keep some of the old, complex code around to perform consistency tests with

[keflavich]

I'm closing this now that your code is incorporated. Thanks!