<s> import statsmodels . api as sm <EOL> import pandas as pd <EOL> import patsy <EOL> import numpy as np <EOL> import warnings <EOL> """<STR_LIT>""" <EOL> _predict_functional_doc = """<STR_LIT>""" <EOL> def _make_exog_from_formula ( result , focus_var , summaries , values , num_points ) : <EOL> """<STR_LIT>""" <EOL> model = result . model <EOL> exog = model . data . frame <EOL> if summaries is None : <EOL> summaries = { } <EOL> if values is None : <EOL> values = { } <EOL> if exog [ focus_var ] . dtype is np . dtype ( '<STR_LIT:O>' ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> colnames = list ( summaries . keys ( ) ) + list ( values . keys ( ) ) + [ focus_var ] <EOL> dtypes = [ exog [ x ] . dtype for x in colnames ] <EOL> varl = set ( exog . columns . tolist ( ) ) - set ( [ model . endog_names ] ) <EOL> unmatched = varl - set ( colnames ) <EOL> unmatched = list ( unmatched ) <EOL> if len ( unmatched ) > <NUM_LIT:0> : <EOL> warnings . warn ( "<STR_LIT>" <EOL> % "<STR_LIT:U+002CU+0020>" . join ( [ "<STR_LIT>" % x for x in unmatched ] ) ) <EOL> ix = range ( num_points ) <EOL> fexog = pd . DataFrame ( index = ix , columns = colnames ) <EOL> for d , x in zip ( dtypes , colnames ) : <EOL> fexog [ x ] = pd . Series ( index = ix , dtype = d ) <EOL> pctls = np . linspace ( <NUM_LIT:0> , <NUM_LIT:100> , num_points ) . tolist ( ) <EOL> fvals = np . percentile ( exog [ focus_var ] , pctls ) <EOL> fvals = np . asarray ( fvals ) <EOL> fexog . loc [ : , focus_var ] = fvals <EOL> for ky in summaries . keys ( ) : <EOL> fexog . loc [ : , ky ] = summaries [ ky ] ( exog . loc [ : , ky ] ) <EOL> for ky in values . keys ( ) : <EOL> fexog . loc [ : , ky ] = values [ ky ] <EOL> dexog = patsy . dmatrix ( model . data . design_info . builder , fexog , return_type = '<STR_LIT>' ) <EOL> return dexog , fexog , fvals <EOL> def _make_exog_from_arrays ( result , focus_var , summaries , values , num_points ) : <EOL> """<STR_LIT>""" <EOL> model = result . model <EOL> model_exog = model . exog <EOL> exog_names = model . exog_names <EOL> if summaries is None : <EOL> summaries = { } <EOL> if values is None : <EOL> values = { } <EOL> exog = np . zeros ( ( num_points , model_exog . shape [ <NUM_LIT:1> ] ) ) <EOL> colnames = list ( values . keys ( ) ) + list ( summaries . keys ( ) ) + [ focus_var ] <EOL> unmatched = set ( exog_names ) - set ( colnames ) <EOL> unmatched = list ( unmatched ) <EOL> if len ( unmatched ) > <NUM_LIT:0> : <EOL> warnings . warn ( "<STR_LIT>" <EOL> % "<STR_LIT:U+002CU+0020>" . join ( [ "<STR_LIT>" % x for x in unmatched ] ) ) <EOL> pctls = np . linspace ( <NUM_LIT:0> , <NUM_LIT:100> , num_points ) . tolist ( ) <EOL> ix = exog_names . index ( focus_var ) <EOL> fvals = np . percentile ( model_exog [ : , ix ] , pctls ) <EOL> exog [ : , ix ] = fvals <EOL> for ky in summaries . keys ( ) : <EOL> ix = exog_names . index ( ky ) <EOL> exog [ : , ix ] = summaries [ ky ] ( model_exog [ : , ix ] ) <EOL> for ky in values . keys ( ) : <EOL> ix = exog_names . index ( ky ) <EOL> exog [ : , ix ] = values [ ky ] <EOL> return exog , fvals <EOL> def _make_exog ( result , focus_var , summaries , values , num_points ) : <EOL> if hasattr ( result . model . data , "<STR_LIT>" ) : <EOL> dexog , fexog , fvals = _make_exog_from_formula ( result , focus_var , <EOL> summaries , values , num_points ) <EOL> else : <EOL> exog , fvals = _make_exog_from_arrays ( result , focus_var , summaries , <EOL> values , num_points ) <EOL> dexog , fexog = exog , exog <EOL> return dexog , fexog , fvals <EOL> def _check_args ( values , summaries , values2 , summaries2 ) : <EOL> if values is None : <EOL> values = { } <EOL> if values2 is None : <EOL> values2 = { } <EOL> if summaries is None : <EOL> summaries = { } <EOL> if summaries2 is None : <EOL> summaries2 = { } <EOL> for ( s , v ) in ( summaries , values ) , ( summaries2 , values2 ) : <EOL> ky = set ( v . keys ( ) ) & set ( s . keys ( ) ) <EOL> ky = list ( ky ) <EOL> if len ( ky ) > <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" + <EOL> "<STR_LIT:U+002CU+0020>" . join ( ky ) ) <EOL> return values , summaries , values2 , summaries2 <EOL> def predict_functional ( result , focus_var , summaries = None , values = None , <EOL> summaries2 = None , values2 = None , alpha = <NUM_LIT> , <EOL> ci_method = "<STR_LIT>" , linear = True , num_points = <NUM_LIT:10> , <EOL> exog = None , exog2 = None , ** kwargs ) : <EOL> if ci_method not in ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> contrast = ( values2 is not None ) or ( summaries2 is not None ) <EOL> if contrast and not linear : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> model = result . model <EOL> if exog is not None : <EOL> if any ( x is not None for x in [ summaries , summaries2 , values , values2 ] ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> fexog = exog <EOL> dexog = patsy . dmatrix ( model . data . design_info . builder , <EOL> fexog , return_type = '<STR_LIT>' ) <EOL> fvals = exog [ focus_var ] <EOL> if exog2 is not None : <EOL> fexog2 = exog <EOL> dexog2 = patsy . dmatrix ( model . data . design_info . builder , <EOL> fexog2 , return_type = '<STR_LIT>' ) <EOL> fvals2 = fvals <EOL> else : <EOL> values , summaries , values2 , summaries2 = _check_args ( values , <EOL> summaries , values2 , summaries2 ) <EOL> dexog , fexog , fvals = _make_exog ( result , focus_var , summaries , values , num_points ) <EOL> if len ( summaries2 ) + len ( values2 ) > <NUM_LIT:0> : <EOL> dexog2 , fexog2 , fvals2 = _make_exog ( result , focus_var , summaries2 , values2 , num_points ) <EOL> from statsmodels . genmod . generalized_linear_model import GLM <EOL> from statsmodels . genmod . generalized_estimating_equations import GEE <EOL> if isinstance ( result . model , ( GLM , GEE ) ) : <EOL> kwargs_pred = kwargs . copy ( ) <EOL> kwargs_pred . update ( { "<STR_LIT>" : True } ) <EOL> else : <EOL> kwargs_pred = kwargs <EOL> pred = result . predict ( exog = fexog , ** kwargs_pred ) <EOL> if contrast : <EOL> pred2 = result . predict ( exog = fexog2 , ** kwargs_pred ) <EOL> pred = pred - pred2 <EOL> dexog = dexog - dexog2 <EOL> if ci_method == '<STR_LIT>' : <EOL> t_test = result . t_test ( dexog ) <EOL> cb = t_test . conf_int ( alpha = alpha ) <EOL> elif ci_method == '<STR_LIT>' : <EOL> t_test = result . t_test ( dexog ) <EOL> sd = t_test . sd <EOL> cb = np . zeros ( ( num_points , <NUM_LIT:2> ) ) <EOL> from scipy . stats . distributions import f as fdist <EOL> df1 = result . model . exog . shape [ <NUM_LIT:1> ] <EOL> df2 = result . model . exog . shape [ <NUM_LIT:0> ] - df1 <EOL> qf = fdist . cdf ( <NUM_LIT:1> - alpha , df1 , df2 ) <EOL> fx = sd * np . sqrt ( df1 * qf ) <EOL> cb [ : , <NUM_LIT:0> ] = pred - fx <EOL> cb [ : , <NUM_LIT:1> ] = pred + fx <EOL> elif ci_method == '<STR_LIT>' : <EOL> sigma , c = _glm_basic_scr ( result , dexog , alpha ) <EOL> cb = np . zeros ( ( dexog . shape [ <NUM_LIT:0> ] , <NUM_LIT:2> ) ) <EOL> cb [ : , <NUM_LIT:0> ] = pred - c * sigma <EOL> cb [ : , <NUM_LIT:1> ] = pred + c * sigma <EOL> if not linear : <EOL> link = result . family . link <EOL> pred = link . inverse ( pred ) <EOL> cb = link . inverse ( cb ) <EOL> return pred , cb , fvals <EOL> predict_functional . __doc__ = _predict_functional_doc <EOL> def _glm_basic_scr ( result , exog , alpha ) : <EOL> """<STR_LIT>""" <EOL> model = result . model <EOL> n = model . exog . shape [ <NUM_LIT:0> ] <EOL> cov = result . cov_params ( ) <EOL> hess = np . linalg . inv ( cov ) <EOL> A = hess / n <EOL> B = np . linalg . cholesky ( A ) . T <EOL> sigma2 = ( np . dot ( exog , cov ) * exog ) . sum ( <NUM_LIT:1> ) <EOL> sigma = np . sqrt ( sigma2 ) <EOL> bz = np . linalg . solve ( B . T , exog . T ) . T <EOL> bz /= np . sqrt ( n ) <EOL> bz /= sigma [ : , None ] <EOL> bzd = np . diff ( bz , <NUM_LIT:1> , axis = <NUM_LIT:0> ) <EOL> bzdn = ( bzd ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) <EOL> kappa_0 = np . sqrt ( bzdn ) . sum ( ) <EOL> from scipy . stats . distributions import norm <EOL> def func ( c ) : <EOL> return kappa_0 * np . exp ( - c ** <NUM_LIT:2> / <NUM_LIT:2> ) / np . pi + <NUM_LIT:2> * ( <NUM_LIT:1> - norm . cdf ( c ) ) - alpha <EOL> from scipy . optimize import brentq <EOL> c , rslt = brentq ( func , <NUM_LIT:1> , <NUM_LIT:10> , full_output = True ) <EOL> if rslt . converged == False : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> return sigma , c </s>
