← Back to branch summary

~scipystats/+junk/skipper-revjp

« back to all changes in this revision

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

  • Committer: joep
  • Date: 2009-08-18 02:04:29 UTC
  • mfrom: (1799.2.14 skipper-working)
  • Revision ID: josef.pktd@gmail.com-20090818020429-td9xt6ay2acg5dr9
merge from skipper, have criss-cross with lot's of conflicts in ols_fstat.py and model.py

Show diffs side-by-side

added added

removed removed

Lines of Context:
nipy/fixes/scipy/stats/models/model.py
1
1
import numpy as np
2
2
from scipy.stats import t, norm
3
 
from scipy import optimize
 
3
from scipy import optimize, derivative
4
4
from models.contrast import ContrastResults
5
5
from models.tools import recipr
6
6
 
39
39
    """
40
40
 
41
41
    def __init__(self, endog, exog=None):
42
 
        self._endog = endog
43
 
        self._exog = exog
 
42
        self._endog = np.asarray(endog)
 
43
        self._exog = np.asarray(exog)
44
44
        self.initialize()
45
45
        
46
46
    def initialize(self):
59
59
 
60
60
    def score(self, params):
61
61
        """
62
 
        Score function of model = gradient of logL with respect to
63
 
        params.
 
62
        Score function of model.
 
63
        
 
64
        The gradient of logL with respect to params.
64
65
        """
65
 
        raise NotImplementedError
66
 
 
 
66
        return derivative(self.loglike, params, dx=1e-04, n=1, order=3)
 
67
                
67
68
    def information(self, params):
68
69
        """
69
70
        Fisher information function of model = - Hessian of logL with respect
84
85
#       so supplied or default guess?
85
86
    def newton(self, params):
86
87
        #JP this is not newton, it's fmin
87
 
# optimize.newton is only for singlevariate?
88
 
# fmin can take multivariate
89
 
# probably called newton bc it's the well known 
90
 
# root finding for MLE
 
88
# it also doesn't work
 
89
# tried fmin, fmin_ncg, fmin_powell
 
90
# converges to wrong estimates
 
91
# will probably need to write own root-finding routine
 
92
#FIXME: should we call it MLE for Maximum Likelihood Estimator?
91
93
# SS no this isn't used anywhere right now, but needs attention for
92
94
# MLE
93
95
        # is this used anywhere
94
96
        # results should be attached to self
95
97
        f = lambda params: -self.loglike(params)
96
 
        xopt, fopt, iter, funcalls, warnflag =\
97
 
          optimize.fmin(f, params, full_output=True)
 
98
        score = lambda params: -self.score(params)
 
99
#        xopt, fopt, iter, funcalls, warnflag =\
 
100
#          optimize.fmin(f, params, full_output=True)
 
101
#        xopt, fopt, fcalls, gcalls, hcalls, warnflag = \
 
102
#                optimize.fmin_ncg(f, params, score) 
98
103
        converge = not warnflag
99
104
        extras = dict(iter=iter, evaluations=funcalls, converge=converge)
100
105
#        return LikelihoodModelResults(self, params, llf=fopt, **extras)
221
226
        -----------
222
227
        r_matrix : array-like
223
228
            Can be 1d, or 2d.  Can be used alone or with other.
224
 
        column :  array-like
 
229
        column :  array-like, optional
225
230
            Must be used on its own.  Can be 0d or 1d see below.
226
 
        scale : float
 
231
        scale : float, optional
227
232
            Can be specified or not.  Default is None, which means that
228
233
            the scale argument is taken from the model.
229
 
        other : array-like
 
234
        other : array-like, optional
230
235
            Can be used when r_matrix is specified.
231
236
 
232
237
        Returns
239
244
        (scale)*(X.T X)^(-1)
240
245
 
241
246
        If contrast is specified it pre and post-multiplies as follows
242
 
        (scale) * contrast (X.T X)^(-1) contrast.T
 
247
        (scale) * r_matrix (X.T X)^(-1) r_matrix.T
243
248
 
244
249
        If contrast and other are specified returns
245
 
        (scale) * contrast (X.T X)^(-1) other.T
 
250
        (scale) * r_matrix (X.T X)^(-1) other.T
246
251
 
247
252
        If column is specified returns
248
253
        (scale) * (X.T X)^(-1)[column,column] if column is 0d 
296
301
        ----------
297
302
        r_matrix : array-like
298
303
            A length p row vector specifying the linear restrictions.
 
304
        scale : float, optional
 
305
            An optional `scale` to use.  Default is the scale specified
 
306
            by the model fit.
299
307
            
300
308
        scale : scalar
301
309
        
338
346
        if self.normalized_cov_params is None:
339
347
            raise ValueError, 'Need covariance of parameters for computing \
340
348
T statistics'
341
 
        if r_matrix.shape[-1] != self.params.shape[0]:
342
 
            raise ValueError, 'r_matrix and params are not aligned'
 
349
        if r_matrix.ndim == 1:
 
350
            if r_matrix.shape[0] != self.params.shape[0]:
 
351
                raise ValueError, 'r_matrix and params are not aligned'
 
352
        elif r_matrix.ndim >1:
 
353
            if r_matrix.shape[1] != self.params.shape[0]:
 
354
                raise ValueError, 'r_matrix and params are not aligned'
343
355
 
344
356
        _t = _sd = None
345
357
 
366
378
        r_matrix : array-like
367
379
            q x p array where q is the number of restrictions to test and
368
380
            p is the number of regressors in the full model fit.
369
 
            If q is 1 then f_test is equivalent to the square of t_test.
370
 
        scale : float
 
381
            If q is 1 then f_test(r_matrix).fvalue is equivalent to the square of 
 
382
            t_test(r_matrix).t
 
383
        scale : float, optional
371
384
            Default is 1.0 for no scaling.
372
 
        invcov : array-like
373
 
            Optional, a qxq matrix to specify an inverse covariance 
 
385
        invcov : array-like, optional
 
386
            A qxq matrix to specify an inverse covariance 
374
387
            matrix based on a restrictions matrix.
375
388
 
376
389
        Examples
430
443
 
431
444
    def conf_int(self, alpha=.05, cols=None):
432
445
        """
433
 
        Returns the confidence interval of the specified params estimates.
 
446
        Returns the confidence interval of the fitted parameters.
434
447
 
435
448
        Parameters
436
449
        ----------
437
450
        alpha : float, optional
438
451
            The `alpha` level for the confidence interval.
439
 
            ie., `alpha` = .05 returns a 95% confidence interval.
440
 
        cols : tuple, optional
 
452
            ie., The default `alpha` = .05 returns a 95% confidence interval.
 
453
        cols : array-like, optional
441
454
            `cols` specifies which confidence intervals to return
442
455
                
443
 
        Returns : array
444
 
            Each item contains [lower, upper]
 
456
        Returns
 
457
        --------
 
458
        conf_int : array
 
459
            Each row contains [lower, upper] confidence interval
445
460
        
446
461
        Example
447
462
        -------
448
 
        >>>import numpy as np
449
 
        >>>from numpy.random import standard_normal as stan
450
 
        >>>import nipy.fixes.scipy.stats.models as SSM
451
 
        >>>x = np.hstack((stan((30,1)),stan((30,1)),stan((30,1))))
452
 
        >>>params=np.array([3.25, 1.5, 7.0])
453
 
        >>>y = np.dot(x,params) + stan((30))
454
 
        >>>model = SSM.regression.OLSModel(x, hascons=False).fit(y)
455
 
        >>>model.conf_int(cols=(1,2))
 
463
        >>>import models
 
464
        >>>from models.datasets.longley.data import load
 
465
        >>>data = load()
 
466
        >>>data.exog = models.tools.add_constant(data.exog)
 
467
        >>>results = models.OLS(data.endog, data.exog).fit()
 
468
        >>>results.conf_int()
 
469
        array([[ -1.77029035e+02,   2.07152780e+02],
 
470
       [ -1.11581102e-01,   3.99427438e-02],
 
471
       [ -3.12506664e+00,  -9.15392966e-01],
 
472
       [ -1.51794870e+00,  -5.48505034e-01],
 
473
       [ -5.62517214e-01,   4.60309003e-01],
 
474
       [  7.98787515e+02,   2.85951541e+03],
 
475
       [ -5.49652948e+06,  -1.46798779e+06]])
 
476
 
 
477
        >>>results.conf_int(cols=(1,2))
 
478
        array([[-0.1115811 ,  0.03994274],
 
479
       [-3.12506664, -0.91539297]])
456
480
 
457
481
        Notes
458
482
        -----
459
 
        TODO:
460
 
        tails : string, optional
461
 
            `tails` can be "two", "upper", or "lower"
 
483
        The confidence interval is based on Student's t distribution for all 
 
484
        models except RLM and GLM, which uses the standard normal distribution.
462
485
        """
463
486
        if self.__class__.__name__ in ['RLMResults','GLMResults']:
464
487
            dist = norm
473
496
            lower = self.params - dist.ppf(1-alpha/2)*self.bse
474
497
            upper = self.params + dist.ppf(1-alpha/2)*self.bse
475
498
        elif cols is not None and dist == t:
476
 
            lower=[]
477
 
            upper=[]
478
 
            for i in cols:
479
 
                lower.append(self.params[i] - dist.ppf(1-\
480
 
                        alpha/2,self.df_resid) *self.bse) 
481
 
                upper.append(self.params[i] + dist.ppf(1-\
482
 
                        alpha/2,self.df_resid) *self.bse)
 
499
            cols = np.asarray(cols)
 
500
            lower = self.params[cols] - dist.ppf(1-\
 
501
                        alpha/2,self.df_resid) *self.bse[cols]
 
502
            upper = self.params[cols] + dist.ppf(1-\
 
503
                        alpha/2,self.df_resid) *self.bse[cols]
483
504
        elif cols is not None and dist == norm:
484
 
            lower = self.params - dist.ppf(1-alpha/2)*self.bse
485
 
            upper = self.params + dist.ppf(1-alpha/2)*self.bse
 
505
            cols = np.asarray(cols)
 
506
            lower = self.params[cols] - dist.ppf(1-alpha/2)*self.bse[cols]
 
507
            upper = self.params[cols] + dist.ppf(1-alpha/2)*self.bse[cols]
486
508
        return np.asarray(zip(lower,upper))
487
509
 
488
510