import os
os.chdir('/Users/toshan/dev/pysal/pysal/contrib/spint')
from count_model import CountModel
import numpy as np
import pandas as pd
from gravity import Gravity, Production, Attraction, Doubly, BaseGravity
import statsmodels.formula.api as smf
from statsmodels.api import families
os.chdir('/Users/toshan/dev/pysal/pysal/contrib/glm')
from glm import GLM
from family import Poisson, QuasiPoisson

import pysal
import os
os.chdir('/Users/toshan/dev/pysal/pysal/contrib/spint')
from dispersion import alpha_disp, phi_disp



rec = pd.read_csv('/Users/toshan/Documents/RecreationDemand.csv')

rec.head()

Unnamed: 0 trips quality ski income userfee costC costS costH
0 1 0 0 yes 4 no 67.59 68.620 76.800
1 2 0 0 no 9 no 68.86 70.936 84.780
2 3 0 0 yes 5 no 58.12 59.465 72.110
3 4 0 0 no 2 no 15.79 13.750 23.680
4 5 0 0 yes 3 no 24.02 34.033 34.547
y = rec['trips'].values.reshape((-1,1))

X = rec[['quality', 'income', 'costC', 'costS', 'costH']].values

test = CountModel(y, X, constant=False)

glm_results = test.fit(framework='glm')

phi_disp(glm_results)

array([  7.30811593e+00,   2.71035909e+00,   3.36051997e-03])
alpha_disp(glm_results)

array([  6.30811593e+00,   2.71035909e+00,   3.36051997e-03])
alpha_disp(glm_results, lambda x: x**2)

array([  1.55402055e+00,   3.38253708e+00,   3.59097912e-04])
#Prepare some test data - columbus example
db = pysal.open(pysal.examples.get_path('columbus.dbf'),'r')
y = np.array(db.by_col("HOVAL"))
y = np.reshape(y, (49,1))
X = []
#X.append(np.ones(len(y)))
X.append(db.by_col("INC"))
X.append(db.by_col("CRIME"))
X = np.array(X).T

poisson_y = np.round(y).astype(int)

test = CountModel(poisson_y, X)

glm_results = test.fit(framework='glm')

phi_disp(glm_results)

array([ 5.39968689,  2.3230411 ,  0.01008847])
alpha_disp(glm_results)

array([ 4.39968689,  2.3230411 ,  0.01008847])
alpha_disp(glm_results, lambda x:x**2)

array([ 0.10690133,  2.24709978,  0.01231683])
model1 = GLM(y, X, constant=False, family=Poisson()).fit()

<class 'family.Poisson'>
<class 'family.Poisson'>
<class 'family.Poisson'>
model2 = GLM(y, X, constant=False, family=QuasiPoisson()).fit()

<class 'family.QuasiPoisson'>
<class 'family.QuasiPoisson'>
<class 'family.QuasiPoisson'>
print model1.scale
print model2.scale

1.0
<class 'family.QuasiPoisson'>
7.02573401193