Source code for statsmodels.miscmodels.tmodel
"""
Linear Model with Student-t distributed errors
Because the t distribution has fatter tails than the normal distribution, it
can be used to model observations with heavier tails and observations that have
some outliers. For the latter case, the t-distribution provides more robust
estimators for mean or mean parameters (what about var?).
Created on 2010-09-24
Author: josef-pktd
License: BSD
References
----------
Kenneth L. Lange, Roderick J. A. Little, Jeremy M. G. Taylor (1989)
Robust Statistical Modeling Using the t Distribution
Journal of the American Statistical Association
Vol. 84, No. 408 (Dec., 1989), pp. 881-896
Published by: American Statistical Association
Stable URL: http://www.jstor.org/stable/2290063
Notes
-----
Reference above not read yet.
TODO:
* add starting values based on OLS
* bugs: store_params does not seem to be defined, I think this was a module
global for debugging - commented out
* parameter restriction: check whether version with some fixed parameters works
"""
# mostly copied from the examples directory written for trying out generic mle.
import numpy as np
from scipy import special, stats
from statsmodels.base.model import GenericLikelihoodModel
# redefine some shortcuts
np_log = np.log
np_pi = np.pi
sps_gamln = special.gammaln
[docs]
class TLinearModel(GenericLikelihoodModel):
"""
Maximum Likelihood Estimation of Linear Model with t-distributed errors
This is an example for generic MLE.
Except for defining the negative log-likelihood method, all
methods and results are generic. Gradients and Hessian
and all resulting statistics are based on numerical
differentiation.
"""
[docs]
def initialize(self):
"""Initialize the model, setting up parameter names and start values"""
print("running Tmodel initialize")
# TODO: here or in __init__
self.k_vars = self.exog.shape[1]
if not hasattr(self, "fix_df"):
self.fix_df = False
if self.fix_df is False:
# df will be estimated, no parameter restrictions
self.fixed_params = None
self.fixed_paramsmask = None
self.k_params = self.exog.shape[1] + 2
extra_params_names = ["df", "scale"]
else:
# df fixed
self.k_params = self.exog.shape[1] + 1
fixdf = np.nan * np.zeros(self.exog.shape[1] + 2)
fixdf[-2] = self.fix_df
self.fixed_params = fixdf
self.fixed_paramsmask = np.isnan(fixdf)
extra_params_names = ["scale"]
super().initialize()
# Note: this needs to be after super initialize
# super initialize sets default df_resid,
# _set_extra_params_names adjusts it
self._set_extra_params_names(extra_params_names)
self._set_start_params()
def _set_start_params(self, start_params=None, use_kurtosis=False):
"""
Set starting values for the parameters
Parameters
----------
start_params : array_like, optional
Starting values to use directly. If None, starting values are
constructed from an OLS fit of `endog` on `exog`.
use_kurtosis : bool
If True and `start_params` is None and `df` is not fixed, use
the kurtosis of the OLS residuals to construct a starting
value for the degrees of freedom parameter. Otherwise a
starting value of 5 is used.
"""
if start_params is not None:
self.start_params = start_params
else:
from statsmodels.regression.linear_model import OLS
res_ols = OLS(self.endog, self.exog).fit()
start_params = 0.1*np.ones(self.k_params)
start_params[:self.k_vars] = res_ols.params
if self.fix_df is False:
if use_kurtosis:
kurt = stats.kurtosis(res_ols.resid)
df = 6./kurt + 4
else:
df = 5
start_params[-2] = df
# TODO adjust scale for df
start_params[-1] = np.sqrt(res_ols.scale)
self.start_params = start_params
[docs]
def loglike(self, params):
"""
Loglikelihood of the model evaluated at params
Parameters
----------
params : ndarray
The parameters of the model. The last 2 parameters are degrees
of freedom and scale.
Returns
-------
float
The log likelihood of the model evaluated at `params`, summed
over all observations.
"""
return -self.nloglikeobs(params).sum(0)
[docs]
def nloglikeobs(self, params):
"""
Negative loglikelihood of linear model with t distributed errors
Parameters
----------
params : ndarray
The parameters of the model. The last 2 parameters are degrees of
freedom and scale.
Returns
-------
ndarray
The negative log likelihood of the model evaluated at `params`
for each observation defined by self.endog and self.exog.
Notes
-----
.. math:: \\ln f(x)=\\ln\\Gamma\\left(\\frac{df+1}{2}\\right)-\\ln\\Gamma\\left(\\frac{df}{2}\\right)-\\frac{1}{2}\\ln(df\\pi)-\\frac{df+1}{2}\\ln\\left(1+\\frac{x^{2}}{df}\\right)-\\ln(scale)
where :math:`x=(y-X\\beta)/scale`.
The t distribution is the standard t distribution and not a standardized
t distribution, which means that the scale parameter is not equal to the
standard deviation.
self.fixed_params and self.expandparams can be used to fix some
parameters. (I doubt this has been tested in this model.)
"""
# print len(params),
# store_params.append(params)
if self.fixed_params is not None:
# print 'using fixed'
params = self.expandparams(params)
beta = params[:-2]
df = params[-2]
scale = np.abs(params[-1]) # TODO check behavior around zero
loc = np.dot(self.exog, beta)
endog = self.endog
x = (endog - loc)/scale
# next part is stats.t._logpdf
lPx = sps_gamln((df+1)/2) - sps_gamln(df/2.)
lPx -= 0.5*np_log(df*np_pi) + (df+1)/2.*np_log(1+(x**2)/df)
lPx -= np_log(scale) # correction for scale
return -lPx
[docs]
def predict(self, params, exog=None):
"""
Return predicted mean values
Parameters
----------
params : ndarray
The parameters of the model. Only the first `exog.shape[1]`
elements (the regression coefficients) are used.
exog : array_like, optional
Explanatory variables to use for prediction. If None,
`self.exog` is used.
Returns
-------
ndarray
The predicted mean, ``exog @ beta``.
"""
if exog is None:
exog = self.exog
return np.dot(exog, params[:self.exog.shape[1]])
