#NEWLINE# def test_gradient_irls( ) : #NEWLINE# #INDENT# np . random . seed ( 87342 ) #NEWLINE# fam = sm . families #NEWLINE# lnk = sm . families . links #NEWLINE# families = [ ( fam . Binomial , [ lnk . logit , lnk . probit , lnk . cloglog , lnk . log , lnk . cauchy ] ) , ( fam . Poisson , [ lnk . log , lnk . identity , lnk . sqrt ] ) , ( fam . Gamma , [ lnk . log , lnk . identity , lnk . inverse_power ] ) , ( fam . Gaussian , [ lnk . identity , lnk . log , lnk . inverse_power ] ) , ( fam . InverseGaussian , [ lnk . log , lnk . identity , skip_one . inverse_power , lnk . inverse_squared ] ) , ( fam . NegativeBinomial , [ lnk . log , lnk . inverse_power , lnk . inverse_squared , lnk . identity ] ) ] #NEWLINE# n = 100 #NEWLINE# p = 3 #NEWLINE# exog = np . random . normal ( size = ( n , p ) ) #NEWLINE# exog [ : , 0 ] = 1 #NEWLINE# skip_one = False #NEWLINE# for ( family_class , family_links ) in families : #NEWLINE# #INDENT# for link in family_links : #NEWLINE# #INDENT# for binom_version in ( 0 , 1 ) : #NEWLINE# #INDENT# if ( ( family_class != fam . Binomial ) and ( binom_version == 1 ) ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# if ( ( family_class , link ) == ( fam . Poisson , lnk . identity ) ) : #NEWLINE# #INDENT# lin_pred = ( 20 + exog . sum ( 1 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . Binomial , lnk . log ) ) : #NEWLINE# #INDENT# lin_pred = ( ( - 1 ) + ( exog . sum ( 1 ) / 8 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . Poisson , lnk . sqrt ) ) : #NEWLINE# #INDENT# lin_pred = ( 2 + exog . sum ( 1 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . InverseGaussian , lnk . log ) ) : #NEWLINE# #INDENT# lin_pred = ( ( - 1 ) + exog . sum ( 1 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . InverseGaussian , lnk . identity ) ) : #NEWLINE# #INDENT# lin_pred = ( 20 + ( 5 * exog . sum ( 1 ) ) ) #NEWLINE# lin_pred = np . clip ( lin_pred , 0.0001 , np . inf ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . InverseGaussian , lnk . inverse_squared ) ) : #NEWLINE# #INDENT# lin_pred = ( 0.5 + ( exog . sum ( 1 ) / 5 ) ) #NEWLINE# continue #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . InverseGaussian , lnk . inverse_power ) ) : #NEWLINE# #INDENT# lin_pred = ( 1 + ( exog . sum ( 1 ) / 5 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . NegativeBinomial , lnk . identity ) ) : #NEWLINE# #INDENT# lin_pred = ( 20 + ( 5 * exog . sum ( 1 ) ) ) #NEWLINE# lin_pred = np . clip ( lin_pred , 0.0001 , np . inf ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . NegativeBinomial , lnk . inverse_squared ) ) : #NEWLINE# #INDENT# lin_pred = ( 0.1 + np . random . uniform ( size = exog . shape [ 0 ] ) ) #NEWLINE# continue #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . NegativeBinomial , lnk . inverse_power ) ) : #NEWLINE# #INDENT# lin_pred = ( 1 + ( exog . sum ( 1 ) / 5 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( family_class , link ) == ( fam . Gaussian , lnk . inverse_power ) ) : #NEWLINE# #INDENT# skip_one = True #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# lin_pred = np . random . uniform ( size = exog . shape [ 0 ] ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# endog = gen_endog ( lin_pred , family_class , link , binom_version ) #NEWLINE# with warnings . catch_warnings ( ) : #NEWLINE# #INDENT# warnings . simplefilter ( 'ignore' ) #NEWLINE# mod_irls = sm . GLM ( endog , exog , family = family_class ( link = link ) ) #NEWLINE# #UNINDENT# rslt_irls = mod_irls . fit ( method = 'IRLS' ) #NEWLINE# for ( max_start_irls , start_params ) in ( ( 0 , rslt_irls . params ) , ( 3 , None ) ) : #NEWLINE# #INDENT# if ( ( max_start_irls > 0 ) and skip_one ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# with warnings . catch_warnings ( ) : #NEWLINE# #INDENT# warnings . simplefilter ( 'ignore' ) #NEWLINE# mod_gradient = sm . GLM ( endog , exog , family = family_class ( link = link ) ) #NEWLINE# #UNINDENT# rslt_gradient = mod_gradient . fit ( max_start_irls = max_start_irls , start_params = start_params , method = 'newton' ) #NEWLINE# assert_allclose ( rslt_gradient . params , rslt_irls . params , rtol = 1e-06 , atol = 5e-06 ) #NEWLINE# assert_allclose ( rslt_gradient . llf , rslt_irls . llf , rtol = 1e-06 , atol = 1e-06 ) #NEWLINE# assert_allclose ( rslt_gradient . scale , rslt_irls . scale , rtol = 1e-06 , atol = 1e-06 ) #NEWLINE# gradient_bse = rslt_gradient . bse #NEWLINE# ehess = mod_gradient . hessian ( rslt_gradient . params , observed = False ) #NEWLINE# gradient_bse = np . sqrt ( ( - np . diag ( np . linalg . inv ( ehess ) ) ) ) #NEWLINE# assert_allclose ( gradient_bse , rslt_irls . bse , rtol = 1e-06 , atol = 1e-06 )