'''gradient/Jacobian of normal and t loglikelihood
use chain rule
normal derivative wrt mu, sigma and beta
new version: loc-scale distributions, derivative wrt loc, scale
also includes "standardized" t distribution (for use in GARCH)
TODO:
* use sympy for derivative of loglike wrt shape parameters
it works for df of t distribution dlog(gamma(a))da = polygamma(0,a) check
polygamma is available in scipy.special
* get loc-scale example to work with mean = X*b
* write some full unit test examples
A: josef-pktd
'''
from __future__ import print_function
import numpy as np
from scipy import special
from scipy.special import gammaln
[docs]def norm_lls(y, params):
'''normal loglikelihood given observations and mean mu and variance sigma2
Parameters
----------
y : array, 1d
normally distributed random variable
params: array, (nobs, 2)
array of mean, variance (mu, sigma2) with observations in rows
Returns
-------
lls : array
contribution to loglikelihood for each observation
'''
mu, sigma2 = params.T
lls = -0.5*(np.log(2*np.pi) + np.log(sigma2) + (y-mu)**2/sigma2)
return lls
[docs]def norm_lls_grad(y, params):
'''Jacobian of normal loglikelihood wrt mean mu and variance sigma2
Parameters
----------
y : array, 1d
normally distributed random variable
params: array, (nobs, 2)
array of mean, variance (mu, sigma2) with observations in rows
Returns
-------
grad : array (nobs, 2)
derivative of loglikelihood for each observation wrt mean in first
column, and wrt variance in second column
Notes
-----
this is actually the derivative wrt sigma not sigma**2, but evaluated
with parameter sigma2 = sigma**2
'''
mu, sigma2 = params.T
dllsdmu = (y-mu)/sigma2
dllsdsigma2 = ((y-mu)**2/sigma2 - 1)/np.sqrt(sigma2)
return np.column_stack((dllsdmu, dllsdsigma2))
[docs]def mean_grad(x, beta):
'''gradient/Jacobian for d (x*beta)/ d beta
'''
return x
[docs]def normgrad(y, x, params):
'''Jacobian of normal loglikelihood wrt mean mu and variance sigma2
Parameters
----------
y : array, 1d
normally distributed random variable with mean x*beta, and variance sigma2
x : array, 2d
explanatory variables, observation in rows, variables in columns
params: array_like, (nvars + 1)
array of coefficients and variance (beta, sigma2)
Returns
-------
grad : array (nobs, 2)
derivative of loglikelihood for each observation wrt mean in first
column, and wrt scale (sigma) in second column
assume params = (beta, sigma2)
Notes
-----
TODO: for heteroscedasticity need sigma to be a 1d array
'''
beta = params[:-1]
sigma2 = params[-1]*np.ones((len(y),1))
dmudbeta = mean_grad(x, beta)
mu = np.dot(x, beta)
#print(beta, sigma2)
params2 = np.column_stack((mu,sigma2))
dllsdms = norm_lls_grad(y,params2)
grad = np.column_stack((dllsdms[:,:1]*dmudbeta, dllsdms[:,:1]))
return grad
[docs]def tstd_lls(y, params, df):
'''t loglikelihood given observations and mean mu and variance sigma2 = 1
Parameters
----------
y : array, 1d
normally distributed random variable
params: array, (nobs, 2)
array of mean, variance (mu, sigma2) with observations in rows
df : integer
degrees of freedom of the t distribution
Returns
-------
lls : array
contribution to loglikelihood for each observation
Notes
-----
parameterized for garch
'''
mu, sigma2 = params.T
df = df*1.0
#lls = gammaln((df+1)/2.) - gammaln(df/2.) - 0.5*np.log((df-2)*np.pi)
#lls -= (df+1)/2. * np.log(1. + (y-mu)**2/(df-2.)/sigma2) + 0.5 * np.log(sigma2)
lls = gammaln((df+1)/2.) - gammaln(df/2.) - 0.5*np.log((df-2)*np.pi)
lls -= (df+1)/2. * np.log(1. + (y-mu)**2/(df-2)/sigma2) + 0.5 * np.log(sigma2)
return lls
[docs]def norm_dlldy(y):
'''derivative of log pdf of standard normal with respect to y
'''
return -y
def ts_dlldy(y, df):
'''derivative of log pdf of standardized (?) t with respect to y
Notes
-----
parameterized for garch, with mean 0 and variance 1
'''
#(df+1)/2. / (1 + y**2/(df-2.)) * 2.*y/(df-2.)
#return -(df+1)/(df-2.) / (1 + y**2/(df-2.)) * y
return -(df+1)/(df) / (1 + y**2/(df)) * y
[docs]def tstd_pdf(x, df):
'''pdf for standardized (not standard) t distribution, variance is one
'''
r = np.array(df*1.0)
Px = np.exp(special.gammaln((r+1)/2.)-special.gammaln(r/2.))/np.sqrt((r-2)*pi)
Px /= (1+(x**2)/(r-2))**((r+1)/2.)
return Px
[docs]def ts_lls(y, params, df):
'''t loglikelihood given observations and mean mu and variance sigma2 = 1
Parameters
----------
y : array, 1d
normally distributed random variable
params: array, (nobs, 2)
array of mean, variance (mu, sigma2) with observations in rows
df : integer
degrees of freedom of the t distribution
Returns
-------
lls : array
contribution to loglikelihood for each observation
Notes
-----
parameterized for garch
normalized/rescaled so that sigma2 is the variance
>>> df = 10; sigma = 1.
>>> stats.t.stats(df, loc=0., scale=sigma.*np.sqrt((df-2.)/df))
(array(0.0), array(1.0))
>>> sigma = np.sqrt(2.)
>>> stats.t.stats(df, loc=0., scale=sigma*np.sqrt((df-2.)/df))
(array(0.0), array(2.0))
'''
print(y, params, df)
mu, sigma2 = params.T
df = df*1.0
#lls = gammaln((df+1)/2.) - gammaln(df/2.) - 0.5*np.log((df-2)*np.pi)
#lls -= (df+1)/2. * np.log(1. + (y-mu)**2/(df-2.)/sigma2) + 0.5 * np.log(sigma2)
lls = gammaln((df+1)/2.) - gammaln(df/2.) - 0.5*np.log((df)*np.pi)
lls -= (df+1.)/2. * np.log(1. + (y-mu)**2/(df)/sigma2) + 0.5 * np.log(sigma2)
return lls
[docs]def ts_dlldy(y, df):
'''derivative of log pdf of standard t with respect to y
Parameters
----------
y : array_like
data points of random variable at which loglike is evaluated
df : array_like
degrees of freedom,shape parameters of log-likelihood function
of t distribution
Returns
-------
dlldy : array
derivative of loglikelihood wrt random variable y evaluated at the
points given in y
Notes
-----
with mean 0 and scale 1, but variance is df/(df-2)
'''
df = df*1.
#(df+1)/2. / (1 + y**2/(df-2.)) * 2.*y/(df-2.)
#return -(df+1)/(df-2.) / (1 + y**2/(df-2.)) * y
return -(df+1)/(df) / (1 + y**2/(df)) * y
[docs]def tstd_dlldy(y, df):
'''derivative of log pdf of standardized t with respect to y
Parameters
----------
y : array_like
data points of random variable at which loglike is evaluated
df : array_like
degrees of freedom,shape parameters of log-likelihood function
of t distribution
Returns
-------
dlldy : array
derivative of loglikelihood wrt random variable y evaluated at the
points given in y
Notes
-----
parameterized for garch, standardized to variance=1
'''
#(df+1)/2. / (1 + y**2/(df-2.)) * 2.*y/(df-2.)
return -(df+1)/(df-2.) / (1 + y**2/(df-2.)) * y
#return (df+1)/(df) / (1 + y**2/(df)) * y
[docs]def locscale_grad(y, loc, scale, dlldy, *args):
'''derivative of log-likelihood with respect to location and scale
Parameters
----------
y : array_like
data points of random variable at which loglike is evaluated
loc : float
location parameter of distribution
scale : float
scale parameter of distribution
dlldy : function
derivative of loglikelihood fuction wrt. random variable x
args : array_like
shape parameters of log-likelihood function
Returns
-------
dlldloc : array
derivative of loglikelihood wrt location evaluated at the
points given in y
dlldscale : array
derivative of loglikelihood wrt scale evaluated at the
points given in y
'''
yst = (y-loc)/scale #ystandardized
dlldloc = -dlldy(yst, *args) / scale
dlldscale = -1./scale - dlldy(yst, *args) * (y-loc)/scale**2
return dlldloc, dlldscale
if __name__ == '__main__':
verbose = 0
if verbose:
sig = 0.1
beta = np.ones(2)
rvs = np.random.randn(10,3)
x = rvs[:,1:]
y = np.dot(x,beta) + sig*rvs[:,0]
params = [1,1,1]
print(normgrad(y, x, params))
dllfdbeta = (y-np.dot(x, beta))[:,None]*x #for sigma = 1
print(dllfdbeta)
print(locscale_grad(y, np.dot(x, beta), 1, norm_dlldy))
print(y-np.dot(x, beta))
from scipy import stats, misc
def llt(y,loc,scale,df):
return np.log(stats.t.pdf(y, df, loc=loc, scale=scale))
def lltloc(loc,y,scale,df):
return np.log(stats.t.pdf(y, df, loc=loc, scale=scale))
def lltscale(scale,y,loc,df):
return np.log(stats.t.pdf(y, df, loc=loc, scale=scale))
def llnorm(y,loc,scale):
return np.log(stats.norm.pdf(y, loc=loc, scale=scale))
def llnormloc(loc,y,scale):
return np.log(stats.norm.pdf(y, loc=loc, scale=scale))
def llnormscale(scale,y,loc):
return np.log(stats.norm.pdf(y, loc=loc, scale=scale))
if verbose:
print('\ngradient of t')
print(misc.derivative(llt, 1, dx=1e-6, n=1, args=(0,1,10), order=3))
print('t ', locscale_grad(1, 0, 1, tstd_dlldy, 10))
print('ts', locscale_grad(1, 0, 1, ts_dlldy, 10))
print(misc.derivative(llt, 1.5, dx=1e-10, n=1, args=(0,1,20), order=3),)
print('ts', locscale_grad(1.5, 0, 1, ts_dlldy, 20))
print(misc.derivative(llt, 1.5, dx=1e-10, n=1, args=(0,2,20), order=3),)
print('ts', locscale_grad(1.5, 0, 2, ts_dlldy, 20))
print(misc.derivative(llt, 1.5, dx=1e-10, n=1, args=(1,2,20), order=3),)
print('ts', locscale_grad(1.5, 1, 2, ts_dlldy, 20))
print(misc.derivative(lltloc, 1, dx=1e-10, n=1, args=(1.5,2,20), order=3),)
print(misc.derivative(lltscale, 2, dx=1e-10, n=1, args=(1.5,1,20), order=3))
y,loc,scale,df = 1.5, 1, 2, 20
print('ts', locscale_grad(y,loc,scale, ts_dlldy, 20))
print(misc.derivative(lltloc, loc, dx=1e-10, n=1, args=(y,scale,df), order=3),)
print(misc.derivative(lltscale, scale, dx=1e-10, n=1, args=(y,loc,df), order=3))
print('\ngradient of norm')
print(misc.derivative(llnorm, 1, dx=1e-6, n=1, args=(0,1), order=3))
print(locscale_grad(1, 0, 1, norm_dlldy))
y,loc,scale = 1.5, 1, 2
print('ts', locscale_grad(y,loc,scale, norm_dlldy))
print(misc.derivative(llnormloc, loc, dx=1e-10, n=1, args=(y,scale), order=3),)
print(misc.derivative(llnormscale, scale, dx=1e-10, n=1, args=(y,loc), order=3))
y,loc,scale = 1.5, 0, 1
print('ts', locscale_grad(y,loc,scale, norm_dlldy))
print(misc.derivative(llnormloc, loc, dx=1e-10, n=1, args=(y,scale), order=3),)
print(misc.derivative(llnormscale, scale, dx=1e-10, n=1, args=(y,loc), order=3))
#print('still something wrong with handling of scale and variance'
#looks ok now
print('\nloglike of t')
print(tstd_lls(1, np.array([0,1]), 100), llt(1,0,1,100), 'differently standardized')
print(tstd_lls(1, np.array([0,1]), 10), llt(1,0,1,10), 'differently standardized')
print(ts_lls(1, np.array([0,1]), 10), llt(1,0,1,10))
print(tstd_lls(1, np.array([0,1.*10./8.]), 10), llt(1.,0,1.,10))
print(ts_lls(1, np.array([0,1]), 100), llt(1,0,1,100))
print(tstd_lls(1, np.array([0,1]), 10), llt(1,0,1.*np.sqrt(8/10.),10))
from numpy.testing import assert_almost_equal
params =[(0, 1), (1.,1.), (0.,2.), ( 1., 2.)]
yt = np.linspace(-2.,2.,11)
for loc,scale in params:
dlldlo = misc.derivative(llnormloc, loc, dx=1e-10, n=1, args=(yt,scale), order=3)
dlldsc = misc.derivative(llnormscale, scale, dx=1e-10, n=1, args=(yt,loc), order=3)
gr = locscale_grad(yt, loc, scale, norm_dlldy)
assert_almost_equal(dlldlo, gr[0], 5, err_msg='deriv loc')
assert_almost_equal(dlldsc, gr[1], 5, err_msg='deriv scale')
for df in [3, 10, 100]:
for loc,scale in params:
dlldlo = misc.derivative(lltloc, loc, dx=1e-10, n=1, args=(yt,scale,df), order=3)
dlldsc = misc.derivative(lltscale, scale, dx=1e-10, n=1, args=(yt,loc,df), order=3)
gr = locscale_grad(yt, loc, scale, ts_dlldy, df)
assert_almost_equal(dlldlo, gr[0], 4, err_msg='deriv loc')
assert_almost_equal(dlldsc, gr[1], 4, err_msg='deriv scale')
assert_almost_equal(ts_lls(yt, np.array([loc, scale**2]), df),
llt(yt,loc,scale,df), 5,
err_msg='loglike')
assert_almost_equal(tstd_lls(yt, np.array([loc, scale**2]), df),
llt(yt,loc,scale*np.sqrt((df-2.)/df),df), 5,
err_msg='loglike')