<commit_before>"""

Test functions for models.GLM
"""

import numpy as np
from numpy.testing import *
import scikits.statsmodels as models
from scikits.statsmodels.glm import GLM
from scikits.statsmodels.tools import add_constant
from nose import SkipTest
from check_for_rpy import skip_rpy

DECIMAL = 4
DECIMAL_less = 3
DECIMAL_lesser = 2
DECIMAL_least = 1
DECIMAL_none = 0
skipR = skip_rpy()
if not skipR:
    from rpy import r
    from rmodelwrap import RModel

class CheckModelResults(object):
    '''
    res2 should be either the results from RModelWrap
    or the results as defined in model_results_data
    '''
    def test_params(self):
        self.check_params(self.res1.params, self.res2.params)

    def test_standard_errors(self):
        assert_almost_equal(self.res1.bse, self.res2.bse, DECIMAL)

    def test_residuals(self):
        if 'rmodelwrap' in self.res2.__module__ and not hasattr(self.res2, 'resids'):
           assert_almost_equal(self.res1.resid_dev, self.res2.resid_dev,
                DECIMAL)
        else:
            resids = np.column_stack((self.res1.resid_pearson,
            self.res1.resid_dev, self.res1.resid_working,
            self.res1.resid_anscombe, self.res1.resid_response))
            self.check_resids(resids, self.res2.resids)

    def test_aic_R(self):
        # R includes the estimation of the scale as a lost dof
        # Doesn't with Gamma though
        if self.res1.scale != 1:
            dof = 2
        else: dof = 0
        self.check_aic_R(self.res1.aic+dof,
                self.res2.aic_R)

    def test_aic_Stata(self):
        if 'rmodelwrap' in self.res2.__module__:
            raise SkipTest("Results are from RModel wrapper")
        aic = self.res1.aic/self.res1.nobs
        self.check_aic_Stata(aic, self.res2.aic_Stata)

    def test_deviance(self):
        assert_almost_equal(self.res1.deviance, self.res2.deviance, DECIMAL)

    def test_scale(self):
        assert_almost_equal(self.res1.scale, self.res2.scale, DECIMAL)

    def test_loglike(self):
        self.check_loglike(self.res1.llf, self.res2.llf)

    def test_null_deviance(self):
        assert_almost_equal(self.res1.null_deviance, self.res2.null_deviance,
                    DECIMAL)

    def test_bic(self):
        if 'rmodelwrap' in self.res2.__module__ and not hasattr(self.res2, 'bic'):
            raise SkipTest("Results are from RModel wrapper")
        self.check_bic(self.res1.bic,
            self.res2.bic)

    def test_degrees(self):
        if not 'rmodelwrap' in self.res2.__module__:
            assert_almost_equal(self.res1.model.df_model,self.res2.df_model, DECIMAL)
        assert_almost_equal(self.res1.model.df_resid,self.res2.df_resid, DECIMAL)

    def test_pearsonX2(self):
        if 'rmodelwrap' in self.res2.__module__:
            raise SkipTest("Results are from RModel wrapper")
        self.check_pearsonX2(self.res1.pearsonX2, self.res2.pearsonX2)

class TestGlmGaussian(CheckModelResults):
    def __init__(self):
        '''
        Test Gaussian family with canonical identity link
        '''

        from scikits.statsmodels.datasets.longley import Load
        self.data = Load()
        self.data.exog = add_constant(self.data.exog)
        self.res1 = GLM(self.data.endog, self.data.exog,
                        family=models.family.Gaussian()).fit()
                                            # I think this is a bug in Rpy

    def setup(self):
        if skipR:
            raise SkipTest, "Rpy not installed."
        Gauss = r.gaussian
        self.res2 = RModel(self.data.endog, self.data.exog, r.glm, family=Gauss)
        self.res2.resids = np.array(self.res2.resid)[:,None]*np.ones((1,5))
        self.res2.null_deviance = 185008826 # taken from R.


    def check_params(self, params1, params2):
        assert_almost_equal(params1, params2, DECIMAL)

    def check_resids(self, resids1, resids2):
        assert_almost_equal(resids1, resids2, DECIMAL)

    def check_aic_R(self, aic1, aic2):
        assert_almost_equal(aic1, aic2, DECIMAL)

    def check_aic_Stata(self, aic1, aic2):
        assert_almost_equal(aic1, aic2, DECIMAL)

    def check_loglike(self, llf1, llf2):
        assert_almost_equal(llf1, llf2, DECIMAL)

    def check_bic(self, bic1, bic2):
        assert_almost_equal(bic1, bic2, DECIMAL)

    def check_pearsonX2(self, pearsonX21, pearsonX22):
        assert_almost_equal(pearsonX21, pearsonX22, DECIMAL)

class TestGaussianLog(CheckModelResults):
    def __init__(self):
        nobs = 100
        x = np.arange(nobs)
        np.random.seed(54321)
#        y = 1.0 - .02*x - .001*x**2 + 0.001 * np.random.randn(nobs)
        X = np.c_[np.ones((nobs,1)),x,x**2]
        lny = np.exp(-(-1.0 + 0.02*x + 0.0001*x**2)) +\
                        0.001 * np.random.randn(nobs)

        GaussLog_Model = GLM(lny, X, \
                family=models.family.Gaussian(models.family.links.log))
        GaussLog_Res = GaussLog_Model.fit()
        self.res1 = GaussLog_Res

    def setup(self):
        if skipR:
            raise SkipTest, "Rpy not installed"
        GaussLogLink = r.gaussian(link = "log")
        GaussLog_Res_R = RModel(lny, X, r.glm, family=GaussLogLink)


    def test_null_deviance(self):
        assert_almost_equal(self.res1.null_deviance, self.res2.null_deviance,
                    DECIMAL_least)<commit_msg>Tests now run with or without rpy installed<commit_after>"""

Test functions for models.GLM
"""

import numpy as np
from numpy.testing import *
import scikits.statsmodels as models
from scikits.statsmodels.glm import GLM
from scikits.statsmodels.tools import add_constant
from nose import SkipTest
from check_for_rpy import skip_rpy

DECIMAL = 4
DECIMAL_less = 3
DECIMAL_lesser = 2
DECIMAL_least = 1
DECIMAL_none = 0
skipR = skip_rpy()
if not skipR:
    from rpy import r
    from rmodelwrap import RModel

class CheckModelResults(object):
    '''
    res2 should be either the results from RModelWrap
    or the results as defined in model_results_data
    '''
    def test_params(self):
        self.check_params(self.res1.params, self.res2.params)

    def test_standard_errors(self):
        assert_almost_equal(self.res1.bse, self.res2.bse, DECIMAL)

    def test_residuals(self):
        if 'rmodelwrap' in self.res2.__module__ and not hasattr(self.res2, 'resids'):
           assert_almost_equal(self.res1.resid_dev, self.res2.resid_dev,
                DECIMAL)
        else:
            resids = np.column_stack((self.res1.resid_pearson,
            self.res1.resid_dev, self.res1.resid_working,
            self.res1.resid_anscombe, self.res1.resid_response))
            self.check_resids(resids, self.res2.resids)

    def test_aic_R(self):
        # R includes the estimation of the scale as a lost dof
        # Doesn't with Gamma though
        if self.res1.scale != 1:
            dof = 2
        else: dof = 0
        self.check_aic_R(self.res1.aic+dof,
                self.res2.aic_R)

    def test_aic_Stata(self):
        if 'rmodelwrap' in self.res2.__module__:
            raise SkipTest("Results are from RModel wrapper")
        aic = self.res1.aic/self.res1.nobs
        self.check_aic_Stata(aic, self.res2.aic_Stata)

    def test_deviance(self):
        assert_almost_equal(self.res1.deviance, self.res2.deviance, DECIMAL)

    def test_scale(self):
        assert_almost_equal(self.res1.scale, self.res2.scale, DECIMAL)

    def test_loglike(self):
        self.check_loglike(self.res1.llf, self.res2.llf)

    def test_null_deviance(self):
        assert_almost_equal(self.res1.null_deviance, self.res2.null_deviance,
                    DECIMAL)

    def test_bic(self):
        if 'rmodelwrap' in self.res2.__module__ and not hasattr(self.res2, 'bic'):
            raise SkipTest("Results are from RModel wrapper")
        self.check_bic(self.res1.bic,
            self.res2.bic)

    def test_degrees(self):
        if not 'rmodelwrap' in self.res2.__module__:
            assert_almost_equal(self.res1.model.df_model,self.res2.df_model, DECIMAL)
        assert_almost_equal(self.res1.model.df_resid,self.res2.df_resid, DECIMAL)

    def test_pearsonX2(self):
        if 'rmodelwrap' in self.res2.__module__:
            raise SkipTest("Results are from RModel wrapper")
        self.check_pearsonX2(self.res1.pearsonX2, self.res2.pearsonX2)

class TestGlmGaussian(CheckModelResults):
    def __init__(self):
        '''
        Test Gaussian family with canonical identity link
        '''

        from scikits.statsmodels.datasets.longley import Load
        self.data = Load()
        self.data.exog = add_constant(self.data.exog)
        self.res1 = GLM(self.data.endog, self.data.exog,
                        family=models.family.Gaussian()).fit()
                                            # I think this is a bug in Rpy

    def setup(self):
        if skipR:
            raise SkipTest, "Rpy not installed."
        Gauss = r.gaussian
        self.res2 = RModel(self.data.endog, self.data.exog, r.glm, family=Gauss)
        self.res2.resids = np.array(self.res2.resid)[:,None]*np.ones((1,5))
        self.res2.null_deviance = 185008826 # taken from R.


    def check_params(self, params1, params2):
        assert_almost_equal(params1, params2, DECIMAL)

    def check_resids(self, resids1, resids2):
        assert_almost_equal(resids1, resids2, DECIMAL)

    def check_aic_R(self, aic1, aic2):
        assert_almost_equal(aic1, aic2, DECIMAL)

    def check_aic_Stata(self, aic1, aic2):
        assert_almost_equal(aic1, aic2, DECIMAL)

    def check_loglike(self, llf1, llf2):
        assert_almost_equal(llf1, llf2, DECIMAL)

    def check_bic(self, bic1, bic2):
        assert_almost_equal(bic1, bic2, DECIMAL)

    def check_pearsonX2(self, pearsonX21, pearsonX22):
        assert_almost_equal(pearsonX21, pearsonX22, DECIMAL)

class TestGaussianLog(CheckModelResults):
    def __init__(self):
        nobs = 100
        x = np.arange(nobs)
        np.random.seed(54321)
#        y = 1.0 - .02*x - .001*x**2 + 0.001 * np.random.randn(nobs)
        self.X = np.c_[np.ones((nobs,1)),x,x**2]
        self.lny = np.exp(-(-1.0 + 0.02*x + 0.0001*x**2)) +\
                        0.001 * np.random.randn(nobs)

        GaussLog_Model = GLM(self.lny, self.X, \
                family=models.family.Gaussian(models.family.links.log))
        GaussLog_Res = GaussLog_Model.fit()
        self.res1 = GaussLog_Res

    def setup(self):
        if skipR:
            raise SkipTest, "Rpy not installed"
        GaussLogLink = r.gaussian(link = "log")
        GaussLog_Res_R = RModel(self.lny, self.X, r.glm, family=GaussLogLink)
        self.res2 = GaussLog_Res_R


    def test_null_deviance(self):
        assert_almost_equal(self.res1.null_deviance, self.res2.null_deviance,
                    DECIMAL_least)
{'fraction_non_alphanumeric': 0.10234107522901822, 'fraction_numerical': 0.026556861293605995, 'mean_word_length': 2.5788079470198677, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5613092', 'n_tokens_mistral': 4027, 'n_tokens_neox': 3718, 'n_words': 736}