<s> """<STR_LIT>""" <EOL> from statsmodels . compat . python import next , range <EOL> import numpy as np <EOL> from statsmodels . genmod import families <EOL> from statsmodels . sandbox . nonparametric . smoothers import PolySmoother <EOL> from statsmodels . genmod . generalized_linear_model import GLM <EOL> from statsmodels . tools . sm_exceptions import IterationLimitWarning , iteration_limit_doc <EOL> import warnings <EOL> DEBUG = False <EOL> def default_smoother ( x , s_arg = None ) : <EOL> '''<STR_LIT:U+0020>''' <EOL> _x = np . sort ( x ) <EOL> n = x . shape [ <NUM_LIT:0> ] <EOL> if n < <NUM_LIT> : <EOL> nknots = n <EOL> else : <EOL> a1 = np . log ( <NUM_LIT:50> ) / np . log ( <NUM_LIT:2> ) <EOL> a2 = np . log ( <NUM_LIT:100> ) / np . log ( <NUM_LIT:2> ) <EOL> a3 = np . log ( <NUM_LIT> ) / np . log ( <NUM_LIT:2> ) <EOL> a4 = np . log ( <NUM_LIT:200> ) / np . log ( <NUM_LIT:2> ) <EOL> if n < <NUM_LIT:200> : <EOL> nknots = <NUM_LIT:2> ** ( a1 + ( a2 - a1 ) * ( n - <NUM_LIT:50> ) / <NUM_LIT> ) <EOL> elif n < <NUM_LIT> : <EOL> nknots = <NUM_LIT:2> ** ( a2 + ( a3 - a2 ) * ( n - <NUM_LIT:200> ) / <NUM_LIT> ) <EOL> elif n < <NUM_LIT> : <EOL> nknots = <NUM_LIT:2> ** ( a3 + ( a4 - a3 ) * ( n - <NUM_LIT> ) / <NUM_LIT> ) <EOL> else : <EOL> nknots = <NUM_LIT:200> + ( n - <NUM_LIT> ) ** <NUM_LIT> <EOL> knots = _x [ np . linspace ( <NUM_LIT:0> , n - <NUM_LIT:1> , nknots ) . astype ( np . int32 ) ] <EOL> if s_arg is None : <EOL> order = <NUM_LIT:3> <EOL> else : <EOL> order = s_arg <EOL> s = PolySmoother ( order , x = x . copy ( ) ) <EOL> return s <EOL> class Offset ( object ) : <EOL> def __init__ ( self , fn , offset ) : <EOL> self . fn = fn <EOL> self . offset = offset <EOL> def __call__ ( self , * args , ** kw ) : <EOL> return self . fn ( * args , ** kw ) + self . offset <EOL> class Results ( object ) : <EOL> def __init__ ( self , Y , alpha , exog , smoothers , family , offset ) : <EOL> self . nobs , self . k_vars = exog . shape <EOL> self . Y = Y <EOL> self . alpha = alpha <EOL> self . smoothers = smoothers <EOL> self . offset = offset <EOL> self . family = family <EOL> self . exog = exog <EOL> self . offset = offset <EOL> self . mu = self . linkinversepredict ( exog ) <EOL> def __call__ ( self , exog ) : <EOL> '''<STR_LIT>''' <EOL> return self . linkinversepredict ( exog ) <EOL> def linkinversepredict ( self , exog ) : <EOL> '''<STR_LIT>''' <EOL> return self . family . link . inverse ( self . predict ( exog ) ) <EOL> def predict ( self , exog ) : <EOL> '''<STR_LIT>''' <EOL> exog_smoothed = self . smoothed ( exog ) <EOL> if exog_smoothed . shape [ <NUM_LIT:0> ] == self . k_vars : <EOL> import warnings <EOL> warnings . warn ( "<STR_LIT>" , <EOL> FutureWarning ) <EOL> return np . sum ( self . smoothed ( exog ) , axis = <NUM_LIT:0> ) + self . alpha <EOL> if exog_smoothed . shape [ <NUM_LIT:1> ] == self . k_vars : <EOL> return np . sum ( exog_smoothed , axis = <NUM_LIT:1> ) + self . alpha <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> def smoothed ( self , exog ) : <EOL> '''<STR_LIT>''' <EOL> return np . array ( [ self . smoothers [ i ] . predict ( exog [ : , i ] ) + self . offset [ i ] <EOL> for i in range ( exog . shape [ <NUM_LIT:1> ] ) ] ) . T <EOL> def smoothed_demeaned ( self , exog ) : <EOL> components = self . smoothed ( exog ) <EOL> means = components . mean ( <NUM_LIT:0> ) <EOL> constant = means . sum ( ) + self . alpha <EOL> components_demeaned = components - means <EOL> return components_demeaned , constant <EOL> class AdditiveModel ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __init__ ( self , exog , smoothers = None , weights = None , family = None ) : <EOL> self . exog = exog <EOL> if not weights is None : <EOL> self . weights = weights <EOL> else : <EOL> self . weights = np . ones ( self . exog . shape [ <NUM_LIT:0> ] ) <EOL> self . smoothers = smoothers or [ default_smoother ( exog [ : , i ] ) for i in range ( exog . shape [ <NUM_LIT:1> ] ) ] <EOL> for i in range ( exog . shape [ <NUM_LIT:1> ] ) : <EOL> self . smoothers [ i ] . df = <NUM_LIT:10> <EOL> if family is None : <EOL> self . family = families . Gaussian ( ) <EOL> else : <EOL> self . family = family <EOL> def _iter__ ( self ) : <EOL> '''<STR_LIT>''' <EOL> self . iter = <NUM_LIT:0> <EOL> self . dev = np . inf <EOL> return self <EOL> def next ( self ) : <EOL> '''<STR_LIT>''' <EOL> _results = self . results <EOL> Y = self . results . Y <EOL> mu = _results . predict ( self . exog ) <EOL> offset = np . zeros ( self . exog . shape [ <NUM_LIT:1> ] , np . float64 ) <EOL> alpha = ( Y * self . weights ) . sum ( ) / self . weights . sum ( ) <EOL> for i in range ( self . exog . shape [ <NUM_LIT:1> ] ) : <EOL> tmp = self . smoothers [ i ] . predict ( ) <EOL> bad = np . isnan ( Y - alpha - mu + tmp ) . any ( ) <EOL> if bad : <EOL> print ( Y , alpha , mu , tmp ) <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> self . smoothers [ i ] . smooth ( Y - mu + tmp , <EOL> weights = self . weights ) <EOL> tmp2 = self . smoothers [ i ] . predict ( ) <EOL> self . results . offset [ i ] = - ( tmp2 * self . weights ) . sum ( ) / self . weights . sum ( ) <EOL> if DEBUG : <EOL> print ( self . smoothers [ i ] . params ) <EOL> mu += tmp2 - tmp <EOL> offset = self . results . offset <EOL> return Results ( Y , alpha , self . exog , self . smoothers , self . family , offset ) <EOL> def cont ( self ) : <EOL> '''<STR_LIT>''' <EOL> self . iter += <NUM_LIT:1> <EOL> if DEBUG : <EOL> print ( self . iter , self . results . Y . shape ) <EOL> print ( self . results . predict ( self . exog ) . shape , self . weights . shape ) <EOL> curdev = ( ( ( self . results . Y - self . results . predict ( self . exog ) ) ** <NUM_LIT:2> ) * self . weights ) . sum ( ) <EOL> if self . iter > self . maxiter : <EOL> return False <EOL> if np . fabs ( ( self . dev - curdev ) / curdev ) < self . rtol : <EOL> self . dev = curdev <EOL> return False <EOL> self . dev = curdev <EOL> return True <EOL> def df_resid ( self ) : <EOL> '''<STR_LIT>''' <EOL> return self . results . Y . shape [ <NUM_LIT:0> ] - np . array ( [ self . smoothers [ i ] . df_fit ( ) for i in range ( self . exog . shape [ <NUM_LIT:1> ] ) ] ) . sum ( ) <EOL> def estimate_scale ( self ) : <EOL> '''<STR_LIT>''' <EOL> return ( ( self . results . Y - self . results ( self . exog ) ) ** <NUM_LIT:2> ) . sum ( ) / self . df_resid ( ) <EOL> def fit ( self , Y , rtol = <NUM_LIT> , maxiter = <NUM_LIT:30> ) : <EOL> '''<STR_LIT>''' <EOL> self . rtol = rtol <EOL> self . maxiter = maxiter <EOL> self . _iter__ ( ) <EOL> mu = <NUM_LIT:0> <EOL> alpha = ( Y * self . weights ) . sum ( ) / self . weights . sum ( ) <EOL> offset = np . zeros ( self . exog . shape [ <NUM_LIT:1> ] , np . float64 ) <EOL> for i in range ( self . exog . shape [ <NUM_LIT:1> ] ) : <EOL> self . smoothers [ i ] . smooth ( Y - alpha - mu , <EOL> weights = self . weights ) <EOL> tmp = self . smoothers [ i ] . predict ( ) <EOL> offset [ i ] = ( tmp * self . weights ) . sum ( ) / self . weights . sum ( ) <EOL> tmp -= tmp . sum ( ) <EOL> mu += tmp <EOL> self . results = Results ( Y , alpha , self . exog , self . smoothers , self . family , offset ) <EOL> while self . cont ( ) : <EOL> self . results = self . next ( ) <EOL> if self . iter >= self . maxiter : <EOL> warnings . warn ( iteration_limit_doc , IterationLimitWarning ) <EOL> return self . results <EOL> class Model ( GLM , AdditiveModel ) : <EOL> def __init__ ( self , endog , exog , smoothers = None , family = families . Gaussian ( ) ) : <EOL> AdditiveModel . __init__ ( self , exog , smoothers = smoothers , family = family ) <EOL> GLM . __init__ ( self , endog , exog , family = family ) <EOL> assert self . family is family <EOL> def next ( self ) : <EOL> _results = self . results <EOL> Y = _results . Y <EOL> if np . isnan ( self . weights ) . all ( ) : <EOL> print ( "<STR_LIT>" ) <EOL> _results . mu = self . family . link . inverse ( _results . predict ( self . exog ) ) <EOL> weights = self . family . weights ( _results . mu ) <EOL> if np . isnan ( weights ) . all ( ) : <EOL> self . weights = weights <EOL> print ( "<STR_LIT>" ) <EOL> self . weights = weights <EOL> if DEBUG : <EOL> print ( '<STR_LIT>' , np . isnan ( self . family . link . deriv ( _results . mu ) ) . any ( ) ) <EOL> Z = _results . predict ( self . exog ) + self . family . link . deriv ( _results . mu ) * ( Y - _results . mu ) <EOL> m = AdditiveModel ( self . exog , smoothers = self . smoothers , <EOL> weights = self . weights , family = self . family ) <EOL> _results = m . fit ( Z ) <EOL> self . history . append ( [ Z , _results . predict ( self . exog ) ] ) <EOL> _results . Y = Y <EOL> _results . mu = self . family . link . inverse ( _results . predict ( self . exog ) ) <EOL> self . iter += <NUM_LIT:1> <EOL> self . results = _results <EOL> return _results <EOL> def estimate_scale ( self , Y = None ) : <EOL> """<STR_LIT>""" <EOL> if Y is None : <EOL> Y = self . Y <EOL> resid = Y - self . results . mu <EOL> return ( np . power ( resid , <NUM_LIT:2> ) / self . family . variance ( self . results . mu ) ) . sum ( ) / self . df_resid <EOL> def fit ( self , Y , rtol = <NUM_LIT> , maxiter = <NUM_LIT:30> ) : <EOL> self . rtol = rtol <EOL> self . maxiter = maxiter <EOL> self . Y = np . asarray ( Y , np . float64 ) <EOL> self . history = [ ] <EOL> self . _iter__ ( ) <EOL> alpha = self . Y . mean ( ) <EOL> mu0 = self . family . starting_mu ( Y ) <EOL> Z = self . family . link ( alpha ) + self . family . link . deriv ( alpha ) * ( Y - mu0 ) <EOL> m = AdditiveModel ( self . exog , smoothers = self . smoothers , family = self . family ) <EOL> self . results = m . fit ( Z ) <EOL> self . results . mu = self . family . link . inverse ( self . results . predict ( self . exog ) ) <EOL> self . results . Y = Y <EOL> while self . cont ( ) : <EOL> self . results = self . next ( ) <EOL> self . scale = self . results . scale = self . estimate_scale ( ) <EOL> if self . iter >= self . maxiter : <EOL> import warnings <EOL> warnings . warn ( iteration_limit_doc , IterationLimitWarning ) <EOL> return self . results </s>
