← Back to branch summary

~scipystats/+junk/skipper-revjp

« back to all changes in this revision

Viewing changes to nipy/fixes/scipy/stats/models/regression.py

  • Committer: joep
  • Date: 2009-08-17 03:43:53 UTC
  • mfrom: (1799.1.26 skipper-working)
  • Revision ID: josef.pktd@gmail.com-20090817034353-ktg588ju3qhmruvq
merge from Skipper, trying to clear up move, keeping history

Show diffs side-by-side

added added

removed removed

Lines of Context:
nipy/fixes/scipy/stats/models/regression.py
24
24
import numpy as np
25
25
from scipy.linalg import norm, toeplitz
26
26
from models.model import LikelihoodModel, LikelihoodModelResults
27
 
from models import utils                # rank utils in np?
 
27
from models import tools
28
28
from models.tools import add_constant
29
29
from scipy import stats, derivative
30
30
from scipy.stats.stats import ss        # could be avoided to eliminate overhead
129
129
        self.calc_params = np.linalg.pinv(self.wdesign)
130
130
        self.normalized_cov_params = np.dot(self.calc_params,
131
131
                                         np.transpose(self.calc_params))
132
 
        self.df_resid = self.wdesign.shape[0] - utils.rank(self._exog)
 
132
        self.df_resid = self.wdesign.shape[0] - tools.rank(self._exog)
133
133
#       Below assumes that we will always have a constant
134
 
        self.df_model = utils.rank(self._exog)-1
 
134
        self.df_model = tools.rank(self._exog)-1
135
135
 
136
136
    def whiten(self, Y):
137
137
        if np.any(self.sigma) and not self.sigma==() :
207
207
        lfit.Z = Z   # not a good name wendog analogy to wdesign
208
208
        lfit.df_resid = self.df_resid
209
209
        lfit.df_model = self.df_model
210
 
        lfit.calc_params = self.calc_params # needed for cov_params()
 
210
        lfit.calc_params = self.calc_params
211
211
        self._summary(lfit)
212
212
        return lfit
213
213
 
214
 
#TODO: make results a property
215
 
# this throws up a set attribute error when running old glm
216
214
    @property
217
215
    def results(self):
218
216
        if self._results is None:
231
229
        lfit.cTSS = ss(lfit.Z-np.mean(lfit.Z))
232
230
#TODO: Z or Y here?  Need to have tests in GLS.
233
231
#JP what does c and u in front of TSS stand for?
 
232
#c is centered and u is uncentered
234
233
#JP I think, it should be Y instead of Z, are the following results correct, with Z?
 
234
#TODO: more robust tests for WLS or GLS, to see if Y or Z is used.
 
235
# I think Y as well, but Z = Y for OLS
235
236
 
236
237
        lfit.uTSS = ss(lfit.Z)
237
238
# Centered R2 for models with intercepts 
250
251
        lfit.F = lfit.MSE_model/lfit.MSE_resid                  
251
252
        lfit.F_p = 1 - stats.f.cdf(lfit.F, lfit.df_model, lfit.df_resid)
252
253
        lfit.bse = np.sqrt(np.diag(lfit.cov_params()))
 
254
        lfit.llf = self.logLike(lfit.params)
 
255
        lfit.aic = -2 * lfit.llf + 2*(self.df_model+1)
 
256
        lfit.bic = -2 * lfit.llf + np.log(lfit.nobs)*(self.df_model+1)
253
257
 
254
 
    def llf(self, params):
 
258
    def logLike(self, params):
255
259
        """
256
260
        Returns the value of the gaussian loglikelihood function at b.
257
261
        
636
640
 
637
641
    It handles the output of contrasts, estimates of covariance, etc.
638
642
    """
639
 
    _llf = None  #this makes it a class attribute - bad, move to init
640
 
    
641
643
    def __init__(self, model, params, normalized_cov_params=None, scale=1.):
642
644
        super(RegressionResults, self).__init__(model, params,
643
645
                                                 normalized_cov_params,
644
646
                                                 scale)
645
 
# should this go up to likelihood model results?
646
 
 
647
 
    @property
648
 
    def llf(self):
649
 
        if self._llf is None:
650
 
            self._llf = self.model.llf(self.params)
651
 
        return self._llf
652
 
 
653
 
    def information_criteria(self):
654
 
        llf = self.llf
655
 
        aic = -2 * llf + 2*(self.df_model + 1)
656
 
        bic = -2 * llf + np.log(self.nobs) * (self.df_model + 1)
657
 
        return dict(aic=aic, bic=bic)
658
 
# could be added as properties to results class.
659
 
 
 
647
#TODO: this needs a test 
660
648
    def norm_resid(self):
661
649
        """
662
650
        Residuals, normalized to have unit length.
678
666
        """
679
667
        if not hasattr(self, 'resid'):
680
668
            raise ValueError, 'need normalized residuals to estimate standard deviation'
681
 
        return self.resid * utils.recipr(np.sqrt(self.scale))
 
669
        return self.resid * tools.recipr(np.sqrt(self.scale))
682
670
 
683
671
    def predictors(self, design):
684
672
        """