Python api.GLS屬性代碼示例

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def _get_sigma(sigma, nobs):
    """
    Returns sigma (matrix, nobs by nobs) for GLS and the inverse of its
    Cholesky decomposition.  Handles dimensions and checks integrity.
    If sigma is None, returns None, None. Otherwise returns sigma,
    cholsigmainv.
    """
    if sigma is None:
        return None, None
    sigma = np.asarray(sigma).squeeze()
    if sigma.ndim == 0:
        sigma = np.repeat(sigma, nobs)
    if sigma.ndim == 1:
        if sigma.shape != (nobs,):
            raise ValueError("Sigma must be a scalar, 1d of length %s or a 2d "
                             "array of shape %s x %s" % (nobs, nobs, nobs))
        cholsigmainv = 1/np.sqrt(sigma)
    else:
        if sigma.shape != (nobs, nobs):
            raise ValueError("Sigma must be a scalar, 1d of length %s or a 2d "
                             "array of shape %s x %s" % (nobs, nobs, nobs))
        cholsigmainv = np.linalg.cholesky(np.linalg.inv(sigma)).T
    return sigma, cholsigmainv 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def reg_m(y, x, estimator, weights=None):
    ones = np.ones(len(x[0]))
    X = sm.add_constant(np.column_stack((x[0], ones)))
    for ele in x[1:]:
        X = sm.add_constant(np.column_stack((ele, X)))

    if estimator=='ols':
        return sm.OLS(y, X).fit()

    elif estimator=='wls':
        return sm.WLS(y, X, weights).fit()

    elif estimator=='gls':
        return sm.GLS(y, X).fit()

    return None

############################
#Run general linear regression
####func array contains the array of functions consdered in the regression
####params coefficients are reversed; the first param coefficient corresponds to the last function in func array
####notice the independent vectors are given in dictionary format, egs:{'bob':[1,2,3,4,5],'mary':[1,2,3,4,5]} 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def _get_sigma(sigma, nobs):
    """
    Returns sigma (matrix, nobs by nobs) for GLS and the inverse of its
    Cholesky decomposition.  Handles dimensions and checks integrity.
    If sigma is None, returns None, None. Otherwise returns sigma,
    cholsigmainv.
    """
    if sigma is None:
        return None, None
    sigma = np.asarray(sigma).squeeze()
    if sigma.ndim == 0:
        sigma = np.repeat(sigma, nobs)
    if sigma.ndim == 1:
        if sigma.shape != (nobs,):
            raise ValueError("Sigma must be a scalar, 1d of length %s or a 2d "
                             "array of shape %s x %s" % (nobs, nobs, nobs))
        cholsigmainv = 1/np.sqrt(sigma)
    else:
        if sigma.shape != (nobs, nobs):
            raise ValueError("Sigma must be a scalar, 1d of length %s or a 2d "
                             "array of shape %s x %s" % (nobs, nobs, nobs))
        cholsigmainv = np.linalg.cholesky(np.linalg.pinv(sigma)).T

    return sigma, cholsigmainv 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def test_pipeline(self):
        from sklearn.feature_selection import SelectKBest
        from sklearn.feature_selection import f_regression
        from sklearn.pipeline import Pipeline

        diabetes = datasets.load_diabetes()
        models = ['OLS', 'GLS', 'WLS', 'GLSAR', 'QuantReg', 'GLM', 'RLM']

        for model in models:
            klass = getattr(sm, model)

            selector = SelectKBest(f_regression, k=5)
            estimator = Pipeline([('selector', selector),
                                  ('reg', base.StatsModelsRegressor(klass))])

            estimator.fit(diabetes.data, diabetes.target)
            result = estimator.predict(diabetes.data)

            data = SelectKBest(f_regression, k=5).fit_transform(diabetes.data, diabetes.target)
            expected = klass(diabetes.target, data).fit().predict(data)
            self.assert_numpy_array_almost_equal(result, expected) 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def __init__(self, endog, exog, sigma=None, missing='none', hasconst=None,
                 **kwargs):
        # TODO: add options igls, for iterative fgls if sigma is None
        # TODO: default if sigma is none should be two-step GLS
        sigma, cholsigmainv = _get_sigma(sigma, len(endog))

        super(GLS, self).__init__(endog, exog, missing=missing,
                                  hasconst=hasconst, sigma=sigma,
                                  cholsigmainv=cholsigmainv, **kwargs)

        # store attribute names for data arrays
        self._data_attr.extend(['sigma', 'cholsigmainv']) 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def whiten(self, X):
        """
        GLS whiten method.

        Parameters
        -----------
        X : array-like
            Data to be whitened.

        Returns
        -------
        np.dot(cholsigmainv,X)

        See Also
        --------
        regression.GLS
        """
        X = np.asarray(X)
        if self.sigma is None or self.sigma.shape == ():
            return X
        elif self.sigma.ndim == 1:
            if X.ndim == 1:
                return X * self.cholsigmainv
            else:
                return X * self.cholsigmainv[:, None]
        else:
            return np.dot(self.cholsigmainv, X) 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def loglike(self, params):
        """
        Returns the value of the Gaussian log-likelihood function at params.

        Given the whitened design matrix, the log-likelihood is evaluated
        at the parameter vector `params` for the dependent variable `endog`.

        Parameters
        ----------
        params : array-like
            The parameter estimates

        Returns
        -------
        loglike : float
            The value of the log-likelihood function for a GLS Model.


        Notes
        -----
        The log-likelihood function for the normal distribution is

        .. math:: -\\frac{n}{2}\\log\\left(\\left(Y-\\hat{Y}\\right)^{\\prime}\\left(Y-\\hat{Y}\\right)\\right)-\\frac{n}{2}\\left(1+\\log\\left(\\frac{2\\pi}{n}\\right)\\right)-\\frac{1}{2}\\log\\left(\\left|\\Sigma\\right|\\right)

        Y and Y-hat are whitened.

        """
        # TODO: combine this with OLS/WLS loglike and add _det_sigma argument
        nobs2 = self.nobs / 2.0
        SSR = np.sum((self.wendog - np.dot(self.wexog, params))**2, axis=0)
        llf = -np.log(SSR) * nobs2      # concentrated likelihood
        llf -= (1+np.log(np.pi/nobs2))*nobs2  # with likelihood constant
        if np.any(self.sigma):
            # FIXME: robust-enough check? unneeded if _det_sigma gets defined
            if self.sigma.ndim == 2:
                det = np.linalg.slogdet(self.sigma)
                llf -= .5*det[1]
            else:
                llf -= 0.5*np.sum(np.log(self.sigma))
            # with error covariance matrix
        return llf 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def test_statsmodels_unknown_constant_position():
    estimator = utils.StatsmodelsSklearnLikeWrapper(
        sm.GLS,
        dict(init=dict(hasconst=True)))
    _, __, estimator = utils.get_regression_model_trainer()(estimator)

    assembler = assemblers.StatsmodelsLinearModelAssembler(estimator)
    assembler.assemble() 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def gls(data, xseq, **params):
    """
    Fit GLS
    """
    if params['formula']:
        return gls_formula(data, xseq, **params)

    X = sm.add_constant(data['x'])
    Xseq = sm.add_constant(xseq)

    init_kwargs, fit_kwargs = separate_method_kwargs(
        params['method_args'], sm.OLS, sm.OLS.fit)
    model = sm.GLS(data['y'], X, **init_kwargs)
    results = model.fit(**fit_kwargs)

    data = pd.DataFrame({'x': xseq})
    data['y'] = results.predict(Xseq)

    if params['se']:
        alpha = 1 - params['level']
        prstd, iv_l, iv_u = wls_prediction_std(
            results, Xseq, alpha=alpha)
        data['se'] = prstd
        data['ymin'] = iv_l
        data['ymax'] = iv_u

    return data 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def gls_formula(data, xseq, **params):
    """
    Fit GLL using a formula
    """
    eval_env = params['enviroment']
    formula = params['formula']
    init_kwargs, fit_kwargs = separate_method_kwargs(
        params['method_args'], sm.GLS, sm.GLS.fit)
    model = smf.gls(
        formula,
        data,
        eval_env=eval_env,
        **init_kwargs
    )
    results = model.fit(**fit_kwargs)
    data = pd.DataFrame({'x': xseq})
    data['y'] = results.predict(data)

    if params['se']:
        _, predictors = dmatrices(formula, data, eval_env=eval_env)
        alpha = 1 - params['level']
        prstd, iv_l, iv_u = wls_prediction_std(
            results, predictors, alpha=alpha)
        data['se'] = prstd
        data['ymin'] = iv_l
        data['ymax'] = iv_u
    return data 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def __init__(self, endog, exog, sigma=None, missing='none', hasconst=None,
                 **kwargs):
    #TODO: add options igls, for iterative fgls if sigma is None
    #TODO: default if sigma is none should be two-step GLS
        sigma, cholsigmainv = _get_sigma(sigma, len(endog))

        super(GLS, self).__init__(endog, exog, missing=missing,
                                  hasconst=hasconst, sigma=sigma,
                                  cholsigmainv=cholsigmainv, **kwargs)

        #store attribute names for data arrays
        self._data_attr.extend(['sigma', 'cholsigmainv']) 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def loglike(self, params):
        """
        Returns the value of the Gaussian log-likelihood function at params.

        Given the whitened design matrix, the log-likelihood is evaluated
        at the parameter vector `params` for the dependent variable `endog`.

        Parameters
        ----------
        params : array-like
            The parameter estimates

        Returns
        -------
        loglike : float
            The value of the log-likelihood function for a GLS Model.


        Notes
        -----
        The log-likelihood function for the normal distribution is

        .. math:: -\\frac{n}{2}\\log\\left(\\left(Y-\\hat{Y}\\right)^{\\prime}\\left(Y-\\hat{Y}\\right)\\right)-\\frac{n}{2}\\left(1+\\log\\left(\\frac{2\\pi}{n}\\right)\\right)-\\frac{1}{2}\\log\\left(\\left|\\Sigma\\right|\\right)

        Y and Y-hat are whitened.

        """
        #TODO: combine this with OLS/WLS loglike and add _det_sigma argument
        nobs2 = self.nobs / 2.0
        SSR = np.sum((self.wendog - np.dot(self.wexog, params))**2, axis=0)
        llf = -np.log(SSR) * nobs2      # concentrated likelihood
        llf -= (1+np.log(np.pi/nobs2))*nobs2  # with likelihood constant
        if np.any(self.sigma):
        #FIXME: robust-enough check?  unneeded if _det_sigma gets defined
            if self.sigma.ndim==2:
                det = np.linalg.slogdet(self.sigma)
                llf -= .5*det[1]
            else:
                llf -= 0.5*np.sum(np.log(self.sigma))
            # with error covariance matrix
        return llf 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def test_Regressions(self):
        diabetes = datasets.load_diabetes()
        models = ['OLS', 'GLS', 'WLS', 'GLSAR', 'QuantReg', 'GLM', 'RLM']

        for model in models:
            klass = getattr(sm, model)

            estimator = base.StatsModelsRegressor(klass)
            estimator.fit(diabetes.data, diabetes.target)
            result = estimator.predict(diabetes.data)

            expected = klass(diabetes.target, diabetes.data).fit().predict(diabetes.data)
            self.assert_numpy_array_almost_equal(result, expected) 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def test_gridsearch(self):
        import sklearn.model_selection as ms
        tuned_parameters = {'statsmodel': [sm.OLS, sm.GLS]}
        diabetes = datasets.load_diabetes()

        cv = ms.GridSearchCV(base.StatsModelsRegressor(sm.OLS), tuned_parameters, cv=5, scoring=None)
        fitted = cv.fit(diabetes.data, diabetes.target)
        self.assertTrue(fitted.best_estimator_.statsmodel is sm.OLS) 

# 需要導入模塊: from statsmodels import api [as 別名]
# 或者: from statsmodels.api import GLS [as 別名]
def iterative_fit(self, maxiter=3, rtol=1e-4, **kwds):
        """
        Perform an iterative two-stage procedure to estimate a GLS model.

        The model is assumed to have AR(p) errors, AR(p) parameters and
        regression coefficients are estimated iteratively.

        Parameters
        ----------
        maxiter : integer, optional
            the number of iterations
        rtol : float, optional
            Relative tolerance between estimated coefficients to stop the
            estimation.  Stops if

            max(abs(last - current) / abs(last)) < rtol

        """
        # TODO: update this after going through example.
        converged = False
        i = -1  # need to initialize for maxiter < 1 (skip loop)
        history = {'params': [], 'rho': [self.rho]}
        for i in range(maxiter - 1):
            if hasattr(self, 'pinv_wexog'):
                del self.pinv_wexog
            self.initialize()
            results = self.fit()
            history['params'].append(results.params)
            if i == 0:
                last = results.params
            else:
                diff = np.max(np.abs(last - results.params) / np.abs(last))
                if diff < rtol:
                    converged = True
                    break
                last = results.params
            self.rho, _ = yule_walker(results.resid,
                                      order=self.order, df=None)
            history['rho'].append(self.rho)

        # why not another call to self.initialize
        # Use kwarg to insert history
        if not converged and maxiter > 0:
            # maxiter <= 0 just does OLS
            if hasattr(self, 'pinv_wexog'):
                del self.pinv_wexog
            self.initialize()

        # if converged then this is a duplicate fit, because we didn't
        # update rho
        results = self.fit(history=history, **kwds)
        results.iter = i + 1
        # add last fit to history, not if duplicate fit
        if not converged:
            results.history['params'].append(results.params)
            results.iter += 1

        results.converged = converged

        return results