

In [12]:

    
import numpy as np
import scipy.stats as ss
from scipy.stats import norm
from scipy.special import erf,erfc
import matplotlib.pyplot as plt
%matplotlib inline

From Reed $$f(x) = \frac{\alpha\beta}{\alpha+\beta}\left[A(\alpha,\nu,\tau^2)x^{-\alpha-1}\Phi\left(\frac{\log x-\nu-\alpha \tau^2}{\tau} \right)+A(-\beta,\nu,\tau^2)x^{\beta-1}\Phi^c\left(\frac{\log x-\nu+\beta \tau^2}{\tau} \right)\right]$$ where $A(\theta,\nu,\tau^2) = \exp\left(\theta \nu + \theta^2\tau^2/2\right)$. Parameter shapes are something like

$\nu$ is like the mean of the lognormal
$\alpha,\beta$ are the PL slopes of the tails
$\tau^2$ is the precision of the normal



In [13]:

    
def dpln(x,alpha=1,beta=0.5,nu=-1,tau2=1):
    A1 = np.exp(alpha*nu+alpha**2*tau2/2)
    A2 = np.exp(-beta*nu+beta**2*tau2/2)
    term1 = A1*x**(-alpha-1)*norm.cdf((np.log(x)-nu-alpha*tau2)/tau2**0.5)
    term2 = A2*x**(beta-1)*norm.sf((np.log(x)-nu+beta*tau2)/tau2**0.5)
    fofx = alpha*beta/(alpha+beta)*(term2+term1)
    return fofx

This also suggests the Power-log normal, which is just one side of the equation (and also likely needs a proper normalization). $$f(x) = \alpha x^{-\alpha-1}\Phi\left(\frac{\log x-\nu-\alpha \tau^2}{\tau} \right)$$



In [27]:

    
def pln(x,alpha=1,tau2=1,nu=1):
    A1 = np.exp(alpha*nu+alpha**2*tau2/2)
    fofx = alpha*x**(-alpha-1)*norm.cdf((np.log(x)-nu-alpha*tau2)/tau2**0.5)
    return fofx



In [29]:

    
trialx = np.logspace(-2,3,100)
pdf = dpln(trialx,alpha=1,nu=-1, beta = 3.0, tau2=0.50)
pdf2 = pln(trialx, alpha=1,nu=-1,tau2=0.5)
plt.loglog(trialx,pdf,label='Double power lognorm')
plt.loglog(trialx,pdf2,label='Power lognorm')
plt.legend(loc=4)









    Out[29]:





<matplotlib.legend.Legend at 0x1157f16d0>

Using the scipy rv_continuous class, we can define dPlN as a subclass.



In [50]:

    
from dpln_distrib import dpln

Let's load in some data to test it...



In [51]:

    
from astropy.io.fits import getdata
img = getdata("../testingdata/hd22.13co.intintensity.fits")
data = img[np.isfinite(img)]

scipy's general definitions are a little weird, and dpln will likely need to be adjusted to some degree (like doing away with nu and tau2 in favour of the built-ins). For now, I've found it necessary to freeze the loc parameter to zero. Also parameter constraints should assist with the fitting.



In [65]:

    
fit = dpln.fit(data, floc=0.0, fscale=1.0)
print(fit)









    



(1.2231409354398668, 1.4283787570663817, 1.3551816091802553, 7.9469599648260317e-10, 0.0, 1.0)

How did it do?



In [66]:

    
plt.hist(data, bins=100, normed=1)
print(dpln.ppf(0.05, *fit), dpln.ppf(0.95, *fit))
x = np.linspace(dpln.ppf(0.05, *fit), dpln.ppf(0.95, *fit), 1000)
plt.plot(x, dpln.pdf(x, *fit))









    



(0.0, 27.07614232728028)






    Out[66]:





[<matplotlib.lines.Line2D at 0x11ad83a90>]

Clearly, there's still some issues with this fit... Now try the single powerlaw form.



In [136]:

    
from pln_distrib import pln



In [143]:

    
fit2 = pln.fit(data, fscale=1.0, floc=0.0)
print(fit2)









    



(1.9312049622065837, 1.0029882024062222, 0.051537141731991747, 0.0, 1.0)



In [144]:

    
plt.hist(data-data.min(), bins=100, normed=1)
print(pln.ppf(0.05, *fit2))
x2 = np.linspace(pln.ppf(0.05, *fit2), pln.ppf(0.95, *fit2), 1000)
plt.plot(x2, pln.pdf(x2, *fit2))









    



0.0






    Out[144]:





[<matplotlib.lines.Line2D at 0x128c41090>]



In [145]:

    
percents = np.linspace(0.01, 1.0, 100)
plt.plot(percents, pln.ppf(percents, *fit2))









    Out[145]:





[<matplotlib.lines.Line2D at 0x129a09c10>]



In [147]:

    
rv = pln(3.0, 1.0, 1.0)
print(rv.ppf(0.01), rv.ppf(0.95))
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
plt.plot(x, rv.pdf(x))









    



(0.0, 21.791615130134865)






    Out[147]:





[<matplotlib.lines.Line2D at 0x129b62350>]



In [148]:

    
plt.plot(x, rv.cdf(x))









    Out[148]:





[<matplotlib.lines.Line2D at 0x129c1a0d0>]



In [149]:

    
plt.plot(x, rv.logpdf(x))









    Out[149]:





[<matplotlib.lines.Line2D at 0x129cb3850>]



In [150]:

    
plt.plot(np.linspace(0.01,1.0,100), rv.ppf(np.linspace(0.01,1.0,100)))









    Out[150]:





[<matplotlib.lines.Line2D at 0x129e6a590>]



In [ ]: