-
Notifications
You must be signed in to change notification settings - Fork 2
/
dpln_distrib.py
44 lines (30 loc) · 1.3 KB
/
dpln_distrib.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
'''
Implementation in scipy form of the Double Pareto-Lognormal Distribution
'''
import numpy as np
from scipy.stats import rv_continuous, norm
def _dpln_pdf(x, alpha, beta, nu, tau2):
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)/np.sqrt(tau2))
term2 = A2*x**(beta-1) * \
norm.sf((np.log(x)-nu+beta*tau2)/np.sqrt(tau2))
return alpha*beta/(alpha+beta)*(term2+term1)
def _dlpn_cdf(x, alpha, beta, nu, tau2):
A1 = np.exp(alpha * nu + alpha**2 * tau2/2)
A2 = np.exp(-beta*nu + beta**2 * tau2/2)
term1 = norm.cdf((np.log(x) - nu) / np.sqrt(tau2))
term2 = beta * x**-alpha * A1 * \
norm.cdf((np.log(x)-nu-alpha * tau2)/np.sqrt(tau2))
term3 = alpha * x**beta * A2 * \
norm.sf((np.log(x)-nu+beta*tau2)/np.sqrt(tau2))
return term1 - (alpha + beta)**-1 * (term2 + term3)
class dpln_gen(rv_continuous):
def _pdf(self, x, alpha, beta, nu, tau2):
return _dpln_pdf(x, alpha, beta, nu, tau2)
def _logpdf(self, x, alpha, beta, nu, tau2):
return np.log(_dpln_pdf(x, alpha, beta, nu, tau2))
def _cdf(self, x, alpha, beta, nu, tau2):
return _dlpn_cdf(x, alpha, beta, nu, tau2)
dpln = dpln_gen(name="dpln", a=0.0)