Python stattools.durbin_watson函数代码示例

def test_durbin_watson():
    #benchmark values from R car::durbinWatsonTest(x)
    #library("car")
    #> durbinWatsonTest(x)
    #[1] 1.95298958377419
    #> durbinWatsonTest(x**2)
    #[1] 1.848802400319998
    #> durbinWatsonTest(x[2:20]+0.5*x[1:19])
    #[1] 1.09897993228779
    #> durbinWatsonTest(x[2:20]+0.8*x[1:19])
    #[1] 0.937241876707273
    #> durbinWatsonTest(x[2:20]+0.9*x[1:19])
    #[1] 0.921488912587806
    st_R = 1.95298958377419
    assert_almost_equal(durbin_watson(x), st_R, 14)

    st_R = 1.848802400319998
    assert_almost_equal(durbin_watson(x**2), st_R, 14)

    st_R = 1.09897993228779
    assert_almost_equal(durbin_watson(x[1:]+0.5*x[:-1]), st_R, 14)

    st_R = 0.937241876707273
    assert_almost_equal(durbin_watson(x[1:]+0.8*x[:-1]), st_R, 14)

    st_R = 0.921488912587806
    assert_almost_equal(durbin_watson(x[1:]+0.9*x[:-1]), st_R, 14)

def checkdb(df,col):
    "It also tells whether the Data is serially correlated or not "
    r = k.durbin_watson(df[col], axis=0)
    if(r==0):
        print "Not Serially correlated "
        return False,r
    else:
        print "Serially correlated "
        return True,r

def fit(y, X, reg_names):
    nr = len(reg_names)
    
    try:
        mod = sm.GLSAR(y.values, X, 2, missing = 'drop') # MLR analysis with AR2 modeling
        res = mod.iterative_fit()
        output = xr.Dataset({'coef': (['reg_name'], res.params[1:]), \
                'conf_int': (['reg_name', 'limit'], res.conf_int()[1:,:]), \
                'p_value': (['reg_name'],  res.pvalues[1:]), \
                'DWT': (sms.durbin_watson(res.wresid)), \
                'CoD': (res.rsquared)}, \
                coords = {'reg_name': (['reg_name'], reg_names),\
                          'limit': (['limit'], ['lower', 'upper'])})
    except: 
        nans = np.full([nr], np.nan)
        output = xr.Dataset({'coef': (['reg_name'], nans), \
                'conf_int': (['reg_name', 'limit'], np.array([nans, nans]).T), \
                'p_value': (['reg_name'],  nans), \
                'DWT': (np.nan), \
                'CoD': (np.nan)}, \
                coords = {'reg_name': (['reg_name'], reg_names),\
                          'limit': (['limit'], ['lower', 'upper'])})

    return output

 def test_durbin_watson_3d(self):
     shape = (10, 1, 10)
     x = np.random.standard_normal(100)
     dw = sum(np.diff(x)**2.0) / np.dot(x, x)
     x = np.tile(x[None, :, None], shape)
     assert_almost_equal(np.squeeze(dw * np.ones(shape)), durbin_watson(x, axis=1))

 def test_durbin_watson(self):
     x = np.random.standard_normal(100)
     dw = sum(np.diff(x)**2.0) / np.dot(x, x)
     assert_almost_equal(dw, durbin_watson(x))

def test_durbin_watson_pandas():
    x = np.random.randn(50)
    x_series = pd.Series(x)
    assert_almost_equal(durbin_watson(x), durbin_watson(x_series), decimal=13)

#.........这里部分代码省略.........
    
    res = scipy.optimize.leastsq(error, pguess, args=(Xdata, Ydata, Errdata, dict_data), full_output=1)
    (popt, pcov, infodict, errmsg, ier) = res
    perr=scipy.sqrt(scipy.diag(pcov))
    
    M=len(Ydata)
    N=len(popt)
    #Residuals
    Y=f(Xdata,popt,dict_data)
    residuals=(Y-Ydata)/Errdata
    meanY=scipy.mean(Ydata)
    squares=(Y-meanY)/Errdata
    squaresT=(Ydata-meanY)/Errdata
    
    SSM=sum(squares**2) #Corrected Sum of Squares
    SSE=sum(residuals**2) #Sum of Squares of Errors
    SST=sum(squaresT**2) #Total corrected sum of squares
    
    DFM=N-1 #Degrees of freedom for model
    DFE=M-N #Degrees of freedom for error
    DFT=M-1 #Degrees of freedom total
    
    MSM=SSM/DFM #Mean squares for model (explained variance)
    MSE=SSE/DFE #Mean squares for Error (should be small wrt MSM) Unexplained Variance
    MST=SST/DFT #Mean squares for total
    
    R2=SSM/SST #proportion of explained variance
    R2_adj=1-(1-R2)*(M-1)/(M-N-1) #Adjusted R2
    
    #t-test to see if parameters are different from zero
    t_stat=popt/perr #t-statistic for popt different from zero'
    t_stat=t_stat.real
    p_p=1.0-scipy.stats.t.cdf(t_stat,DFE) #should be low for good fit.
    z=scipy.stats.t(M-N).ppf(0.95)
    p95=perr*z
    #Chisquared Analysis on Residuals
    chisquared=sum(residuals**2)
    degfreedom=M-N
    chisquared_red=chisquared/degfreedom
    p_chi2=1.0-scipy.stats.chi2.cdf(chisquared,degfreedom)
    stderr_reg=scipy.sqrt(chisquared_red)
    chisquare=(p_chi2,chisquared, chisquared_red, degfreedom,R2,R2_adj)
    
    #Analysis of Residuals
    w,p_shapiro=scipy.stats.shapiro(residuals)
    mean_res=scipy.mean(residuals)
    stddev_res=scipy.sqrt(scipy.var(residuals))
    t_res=mean_res/stddev_res #t-statistic to test that mean_res is zero.
    p_res=1.0-scipy.stats.t.cdf(t_res,M-1)
        #if p_res <0.05, null hypothesis rejected and mean is non-zero.
        #Should be high for good fit.
    #F-test on residuals
    F=MSM/MSE #explained variance/unexplained . Should be large
    p_F=1.0-scipy.stats.f.cdf(F,DFM,DFE)
        #if p_F <0.05n, null-hypothesis is rejected
        #i.e. R^2>0 and at least one of the fitting parameters >0.
    dw=stools.durbin_watson(residuals)
    resanal=(p_shapiro,w,mean_res,F,p_F,dw)
    
    
    
    if ax:
        formataxis(ax)
        ax.plot(Ydata,Y,'ro')
        ax.errorbar(Ydata,Y,yerr=Errdata,fmt='.')
        Ymin,Ymax=min((min(Y),min(Ydata))),max((max(Y),max(Ydata)))
        ax.plot([Ymin,Ymax],[Ymin,Ymax],'b')
        
        ax.xaxis.label.set_text('Data')
        ax.yaxis.label.set_text('Fitted')
        
        sigmay,avg_stddev_data=get_stderr_fit(f,Xdata,popt,pcov,dict_data)
        Yplus=Y+sigmay
        Yminus=Y-sigmay
        ax.plot(Y,Yplus,'c',alpha=0.6,linestyle='--',linewidth=0.5)
        ax.plot(Y,Yminus,'c',alpha=0.6,linestyle='--',linewidth=0.5)
        ax.fill_between(Y,Yminus,Yplus,facecolor='cyan',alpha=0.5)
        titletext='Parity plot for fit.\n'
        titletext+=r'$r^2$=%5.2f, $r^2_{adj}$=%5.2f, '
        titletext +='$\sigma_{exp}$=%5.2f,$\chi^2_{\nu}$=%5.2f, $p_{\chi^2}$=%5.2f, '
        titletext +='$\sigma_{err}^{reg}$=%5.2f'
        
        ax.title.set_text(titletext%(R2,R2_adj, avg_stddev_data, chisquared_red,  p_chi2, stderr_reg))
        ax.figure.canvas.draw()
        
    if ax2:#Test for homoscedasticity
        formataxis(ax2)
        ax2.plot(Y,residuals,'ro')
        
        ax2.xaxis.label.set_text('Fitted Data')
        ax2.yaxis.label.set_text('Residuals')
        
        titletext='Analysis of Residuals\n'
        titletext+=r'mean=%5.2f, $p_{res}$=%5.2f, $p_{shapiro}$=%5.2f, $Durbin-Watson$=%2.1f'
        titletext+='\nF=%5.2f, $p_F$=%3.2e'
        ax2.title.set_text(titletext%(mean_res, p_res, p_shapiro, dw , F, p_F))
        
        ax2.figure.canvas.draw()
        
    return popt,pcov,perr,p95,p_p,chisquare,resanal

    def test_frame_timeseries_durbin_watson(self):
        """Test Durbin Watson"""
        result = self.frame.timeseries_durbin_watson_test("logM")
        db_result = smst.durbin_watson(self.pandaframe["logM"])

        self.assertAlmostEqual(result, db_result, delta=0.0000000001)

def fitdata(f,Xdata,Ydata,Errdata,pguess,ax= False,ax2= False):
    
    #calculating the popt 
    
    def error(pguess,Xdata,Ydata,Errdata):
        Y=f(Xdata,pguess)
        residuals= (Y-Ydata)/Errdata
        return (residuals)
    
    res= scipy.optimize.leastsq(error, pguess,args=(Xdata,Ydata,Errdata),full_output=1)
      
    (popt,pcov,infodict,errmsg,ier)=res
      
    perr= scipy.sqrt(scipy.diag(pcov))
    
    M= len(Xdata)
    N= len(popt)
   #residuals
    Y= f(Xdata,popt)
    residuals=(Y-Ydata)/Errdata
    meanY= scipy.mean(Ydata)
    squares= (Y-meanY)/Errdata
    squaresT= (Ydata-meanY)/Errdata
    
    SSM= sum(squares**2)#corrected sum of squares
    SSE= sum(residuals**2)#sum of squares of errors
    SST= sum(squaresT**2)#total corrected sum of squrare

    DFM= N-1
    DFE= M-N
    DFT= M-1

    MSM= SSM/DFM
    MSE= SSE/DFE
    MST= SST/DFT

    R2= SSM/SST    #proportion of explained variance
    R2_adj= 1-(1-R2)*(M-1)/(M-N-1)#Adjusted R2 

    # t test to see if parameters are different from 0
    t_stat= popt/perr
    t_stat= t_stat.real
    p_p= 1.0-scipy.stats.t.cdf(t_stat,DFE)
    z=scipy.stats.t(M-N).ppf(0.95)
    p95= perr*z

    #chisquared analysis on residuals
    chisquared= sum(residuals**2)
    degfreedom= M-N
    chisquared_red= chisquared/degfreedom
    p_chi2= 1.0-scipy.stats.chi2.cdf(chisquared,degfreedom)       
    stderr_reg= scipy.sqrt(chisquared_red)
    chisquare=(p_chi2, chisquared,chisquared_red,degfreedom,R2,R2_adj)
    
    #analysis of residuals
    w,p_shapiro= scipy.stats.shapiro(residuals)
    mean_res= scipy.mean(residuals)
    stddev_res= scipy.sqrt(scipy.var(residuals))
    t_res= mean_res/stddev_res 
    p_res=1.0-scipy.stats.t.cdf(t_res,M-1)
    #if p<0.05, null hypothesis is rejected and mean is non-zero
    #should be high for a good fit
    
    #F-test on the residuals
    F= MSM/MSE #explained variance/ unexplained should be large
    p_F= 1.0-scipy.stats.f.cdf(F,DFM,DFE)
    #if p_F<0.05, null hypo is rejected
    dw= stools.durbin_watson(residuals)
    resanal= (p_shapiro,w,mean_res,p_res,F,p_F,dw)
    
    if ax:
        formataxis(ax)
        ax.plot(Ydata,Y,'ro')
        ax.errorbar(Ydata,Y,yerr=Errdata,fmt='.')
        Ymin, Ymax= min((min(Y),min(Ydata))),max((max(Y),max(Ydata)))
        ax.plot([Ymin,Ymax],[Ymin,Ymax],'b')
        
        ax.xaxis.label.set_text('Data')
        ax.yaxis.label.set_text('Fitted')
        
        sigmaY, avg_stddev_data= get_stderr_fit(f,Xdata,popt,pcov)
        Yplus= Y+sigmaY
        Yminus= Y-sigmaY
        ax.plot(Y,Yplus,'c',alpha=0.6,linestyle='--',linewidth= 0.5)
        ax.plot(Y,Yminus,'c',alpha=0.6,linestyle='--',linewidth= 0.5)
        ax.fill_between(Y,Yminus,Yplus,facecolor= 'cyan',alpha=0.5)
        titletext='Parity plot for fit.\n'
        titletext+= r'$r^r$=%5.2f,$r^2_(adj)$=%5.2f,'
        titletext+= '$\sigma_<exp>$=%5.2f,$\chi^2_<\nu>$= %5.2f,$p_<chi_2>$=%5.2f,'
        titletext+= '$sigma_<err>^<reg>$=%5.2f'
        
        ax.title.set_text(titletext%(R2,R2_adj,avg_stddev_data,chisquared_red,p_chi2,stderr_reg))
        ax.figure.canvas.draw()
        
    if ax2 :  #test for homoscedaticity  
        formataxis(ax2)
        ax2.plot(Y,residuals,'ro')
        
        ax2.xaxis.label.set_text('Fitted Data')
        ax2.yaxis.label.set_text('Residuals')
#.........这里部分代码省略.........

def fitdata(f,Xdata,Ydata,Errdata,pguess,ax=False,ax2=False):
	'''
	fitdata(f,Xdata,Ydata,Errdata,pguess):
	'''
	def error(p,Xdata,Ydata,Errdata):
		Y=f(Xdata,p)
		residuals=(Y-Ydata)/Errdata
		return residuals
	res=scipy.optimize.leastsq(error,pguess,args=(Xdata,Ydata,Errdata),full_output=1)
	(popt,pcov,infodict,errmsg,ier)=res   #optimize p
	perr=scipy.sqrt(scipy.diag(pcov))    #vector of sd of p
	M=len(Ydata)
	N=len(popt)
	#Residuals
	Y=f(Xdata,popt)
	residuals=(Y-Ydata)/Errdata
	meanY=scipy.mean(Ydata)
	squares=(Y-meanY)/Errdata
	squaresT=(Ydata-meanY)/Errdata
	
	SSM=sum(squares**2) #corrected sum of squares
	SSE=sum(residuals**2)  #sum of squares of errors
	SST=sum(squaresT**2)  #total corrected sum of squares
	
	DFM=N-1   #for model
	DFE=M-N   #for error
	DFT=M-1   #total
	
	MSM=SSM/DFM #mean squares for model(explained variance)
	MSE=SSE/DFE #mean squares for errors(should be small wrt MSM) unexplained variance
	MST=SST/DFT #mean squares for total	
	
	R2=SSM/SST  #proportion of explained variance
	R2_adj=1-(1-R2)*(M-1)/(M-N-1) #adjusted R2
	
	#ttest to see if parameters are different from zero
	t_stat=popt/perr #tstatistic for popt different from zero
	t_stat=t_stat.real
	p_p=1.0-scipy.stats.t.cdf(t_stat,DFE) #should be low for good fit
	z=scipy.stats.t(M-N).ppf(0.95)
	p95=perr*z
	#Chisquared ananlysis on residuals
	chisquared=sum(residuals**2)
	degfreedom=M-N
	chisquared_red=chisquared/degfreedom
	p_chi2=1.0-scipy.stats.chi2.cdf(chisquared,degfreedom)
	stderr_reg=scipy.sqrt(chisquared_red)
	chisquare=(p_chi2,chisquared,chisquared_red,degfreedom,R2,R2_adj)
	
	#Analysis of residuals
	w,p_shapiro=scipy.stats.shapiro(residuals) # to check if residuals are normally distributed
	mean_res=scipy.mean(residuals)
	stddev_res=scipy.sqrt(scipy.var(residuals))
	t_res=mean_res/stddev_res
	p_res=1.0-scipy.stats.t.cdf(t_res,M-1)
	#if p_res<0.05,null hypothesis is rejected.
	#R^2>0 and at least one of the fitting parameters>0
	#F-test on residuals
	F=MSM/MSE
	p_F=1.0-scipy.stats.f.cdf(F,DFM,DFE)
	
	dw=stools.durbin_watson(residuals) #to check if they are correlated
	resanal=(p_shapiro,w,mean_res,p_res,F,p_F,dw)
	
	if ax:
		formataxis(ax)
		ax.plot(Ydata,Y,'ro')
		ax.errorbar(Ydata,Y,yerr=Errdata,fmt='.')
		Ymin,Ymax=min((min(Y),min(Ydata))),max((max(Y),max(Ydata)))
		ax.plot([Ymin,Ymax],[Ymin,Ymax],'b')
		
		ax.xaxis.label.set_text('Data')
		ax.yaxis.label.set_text('Fitted')
		sigmay,avg_stddev_data=get_stderr_fit(f, Xdata, popt, pcov)
		Yplus=Y+sigmay
		Yminus=Y-sigmay
		ax.plot(Y,Yplus,'c',alpha=.6,linestyle='--',linewidth=.5)
		ax.plot(Y,Yminus,'c',alpha=.6,linestyle='--',linewidth=.5)
		ax.fill_between(Y,Yminus,Yplus,facecolor='cyan',alpha=.5)
		titletext='Parity plot for fit.\n'
		titletext+=r'$r^2$=%5.2f,$r^2_{adj}$=%5.2f, '
		titletext+='$\sigma_{exp}$=%5.2f,$\chi^2_{\nu}$=%5.2f,$p_{\chi^2}$=%5.2f, '
		titletext+='$\sigma_{err}^{reg}$=%5.2f'
		ax.title.set_text(titletext%(R2,R2_adj,avg_stddev_data,chisquared_red,p_chi2,stderr_reg))
		ax.figure.canvas.draw()
	
	if ax2:#test for homoscedasticity
		formataxis(ax2)
		ax2.plot(Y,residuals,'ro')
		
		ax2.xaxis.label.set_text('Fitted data')
		ax2.yaxis.label.set_text('Residuals')
		
		titletext='Analysis of residuals\n'
		titletext+=r'mean=%5.2f,$p_{res}$=%5.2f,$p_{shapiro}$=%5.2f,$Durbin-Watson$=%2.1f'
		titletext+='\n F=%5.2f,$p_F$=%3.2e'
		ax2.title.set_text(titletext%(mean_res,p_res,p_shapiro,dw,F,p_F))
		
		ax2.figure.canvas.draw()
		
#.........这里部分代码省略.........

 def test_durbin_watson_2d(self, reset_randomstate):
     shape = (1, 10)
     x = np.random.standard_normal(100)
     dw = sum(np.diff(x)**2.0) / np.dot(x, x)
     x = np.tile(x[:, None], shape)
     assert_almost_equal(np.squeeze(dw * np.ones(shape)), durbin_watson(x))

sunspot_data = sm.datasets.sunspots.load_pandas().data
sunspot_data.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))
del sunspot_data["YEAR"]

sunspot_data.plot(figsize=(12,8));

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(sunspot_data.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(sunspot_data, lags=40, ax=ax2)

arma_mod20 = sm.tsa.ARMA(sunspot_data, (2,0)).fit()
arma_mod30 = sm.tsa.ARMA(sunspot_data, (3,0)).fit()

stattools.durbin_watson(arma_mod30.resid.values)

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax = arma_mod30.resid.plot(ax=ax)

resid = arma_mod30.resid

diag.normal_ad(resid)

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)

def fitdata(f,Xdata,Ydata,Errdata,pguess,dict_data,ax=False,ax2=False):
    
    def error(p,Xdata,Ydata,Errdata,dict_data):
        Y=f(Xdata,p,dicct_data)
        residuals=(Y-Ydata)/Errdata
        return residuals
    res=scipy.optimize.leastsq(error,pguess.args=(Xdata,Ydata,Errdata,dict_data),full_output=1)
    (popt,pcov,infodict,errmsg,ier)=res
    perr=scipy.sqrt(scipy.diag(pcov))
    
    M=len(Ydata)
    N=len(popt)
    
    Y=f(Xdata,popt,dict_data)
    residuals=(Y-Ydata)/Errdata
    meanY=scipy.mean(Ydata)
    squares=(Y-meanY)/Errdata
    squaresT=(Ydata-meanY)/Errdata
    
    SSM=sum(squares**2)
    SSE=sum(residuals**2)
    SST=sum(squaresT**2)
    
    DFM=N-1
    DFE=M-N
    DFT=M-1
    
    MSM=SSM/DFM
    MSE=SSE/DFE
    MST=SST/DFT
    
    R2=SSM/SST
    R2_adj=1-(1-R2)*(M-1)/(M-N-1)
    
    t_stat=popt/perr
    t_stat=t_stat.real
    p_p=1.0-scipy.stats.t.cdf(t_stat,DFE)
    z=scipy.stats.t(M-N).ppf(0.95)
    p95=perr*z
    
    chisquared=sum(residuals**2)
    degfreeedom=M-N
    chisqured_red=chisquared/degfreedom
    p_chi2=1.0-scipy.stats.chi2.cdf(chisquared,degfreedom)
    chisquare=(p_chi2,chisquared,chisquared_red,degfreedom,R2,R2_adj)
    
    w,p_shapiro=scipy.stats.shapiro(residuals)
    mean_res=scipy.mean(residuals)
    stddev_res=scipy.sqrt(scipy.var(residuals))
    t_res=mean_res/stddev_res
    p_res=1.0-scipy.stats.t.cdf(t_res,M-1)
    
    F=MSM/MSE
    p_F=1.0-scipy.stats.f.cdf(F,DFM,DFE)
    
    dw=stools.durbin_watson(residuals)
    resanal=(p_shapiro,w,mean_res,p_res,F,p_F,dw)
    
    if ax:
        formataxis(ax)
        ax.plotdata(Ydata,Y,'ro')
        ax.errorbar(Ydata,Y,yerr=Errdata,fmt='.')
        Ymin,Ymax=min((min(Y),min(Ydata))),max((max(Y),max(Ydata)))
        ax.plot([Ymin,Ymax],[Ymin,Ymax],'b')
        
        ax.xaxis.label.set_text('Data')
        ax.yaxis.label.set_text('Fitted')
        
        sigmay,avg_stddev_data=get_stderr_fit(f,Xdata,popt,pcov,dict_data)
        Yplus=Y+sigmay
        Yminus=Y-sigmay
        ax.plot(Y,Yplus,'c',alpha=0.6,linestyle='--',linewidth=0.5)
        ax.plot(Y,Yminus,'c',alpha=0.6,linestyle='--',linewidth=0.5)
        ax.fill_between(Y,Yminus,Yplus,facecolor='cyan',alpha=0.5)
        titletext ='parity plot for fit.\n'
        titletext +=r'$r^2$=%5.2f,$r^2_{adj}$=%5.2f,'
        titletext +='$\sigma_{exp}$=%5.2f,$\chi^2_{\nu}$=%5.2f,$p_{\chi^2}$=%5.2f,'
        titletext +='$\sigma_{err}^{reg}$=%5.2f'
        
        ax.title.set_text(titletext%(R2,R2_adj,avg_stddev_data,chisquared_red,p_chi2,stderr_reg))
        ax.figure.canvas.draw()
        
    if ax2:
        formataxis(ax2)
        ax2.plot(Y,residuals,'ro')
        
        ax2.xaxis.label.set_text('Fitted Data')
        ax2.yaxis.label.set_text('Residuals')
        
        titletext ='Analysis of Residuals\n'
        titletext +=r'mean=%5.2f,$p_{res}$=%5.2f,$p_{shapiro}$=%5.2f, $Durbin-Watson=%2.1f'
        titletext +='\n F=%5.2f, $p_F$=%3.2e'
        ax2.title.set_text(titletext%(mean_res,p_res,p_shapiro,dw,F,p_F))
        
        ax2.figure.canvas.draw()
        
    return popt,pcov,perr,p95,p_p,chisquare,resanal

def fitdata(f,Xdata,Ydata,Errdata,pguess,ax=False,ax2=False):

    def error(p,Xdata,Ydata,Errdata):
        Y=f(Xdata,p)
        residuals=(Y-Ydata)/Errdata
        return residuals
    res=scipy.optimize.leastsq(error,pguess,args=(Xdata,Ydata,Errdata),full_output=1)    
    (popt,pcov,infodict,errmsg,ier)=res
    perr=scipy.sqrt(scipy.diag(pcov))
    M=len(Ydata)
    N=len(popt)
    Y=f(Xdata,popt)
    residuals=(Y-Ydata)/Errdata
    meanY=scipy.mean(Ydata)
    squares=(Y-meanY)/Errdata
    squaresT=(Ydata-meanY)/Errdata
    
    SSM=sum(squares**2)
    SSE=sum(residuals**2)
    SST=sum(squaresT**2)
    
    DFM=N-1
    DFE=M-N
    DFT=M-1
    
    MSM=SSM/DFM
    MSE=SSE/DFE
    MSM=SST/DFT
    
    '''R2'''
    R2=SSM/SST
    R2_adj=1-(1-R2)*(M-1)/(M-N-1)
    
    '''t-test'''
    t_stat=popt/perr
    t_stat=t_stat.real
    p_p=1.0-scipy.stats.t.cdf(t_stat,DFE)
    z=scipy.stats.t(M-N).ppf(0.95)
    p95=perr*z
    
    '''chi-square'''
    chisqred=sum(residuals**2)
    degfrdm=M-N
    chisqred_red=chisqred/degfrdm
    p_chi2=1.0-scipy.stats.chi2.cdf(chisqred,degfrdm)
    stderr_reg=scipy.sqrt(chisqred_red)
    chisqre=(p_chi2,chisqred,chisqred_red,degfrdm,R2,R2_adj)
    
    '''shapiro-wilk test'''
    w,p_shapiro=scipy.stats.shapiro(residuals)
    mean_res=scipy.mean(residuals)
    stddev_res=scipy.sqrt(scipy.var(residuals))
    t_res=mean_res/stddev_res
    p_res=1.0-scipy.stats.t.cdf(t_res,M-1)
    
    '''F-test'''
    F=MSM/MSE
    p_F=1-scipy.stats.f.cdf(F,DFM,DFE)
    
    '''durbin-watson'''
    dw=stools.durbin_watson(residuals)
    
    resanal=(p_shapiro,w,mean_res,p_res,F,p_F,dw)

    return popt,pcov,perr,p95,p_p,chisqre,resanal,R2,chisqred,w,dw

        'F_Statistic': results.fvalue,
        'F_Statistic_P_Value': results.f_pvalue,
        'Log_Likelihood': results.llf,
        'AIC': results.aic,
        'BIC': results.bic,
        'Number_Of_Observations': int(results.nobs),
        'Degrees_Of_Freedom_Model': int(results.df_model),
        'Degrees_Of_Freedom_Residual': int(results.df_resid)
    },
    'Parameters': params,
    'Diagnostics': {
        'Omnibus': results.diagn['omni'],
        'Omnibus_P_Value': results.diagn['omnipv'],
        'Skew': results.diagn['skew'],
        'Kurtosis': results.diagn['kurtosis'],
        'Durbin_Watson': durbin_watson(results.wresid),
        'Jarque_Bera': results.diagn['jb'],
        'Jarque_Bera_P_Value': results.diagn['jbpv'],
        'Condition_Number': results.diagn['condno']
    }
}

print(json.dumps(resultsObject, sort_keys=True))


#                             OLS Regression Results                            
# ==============================================================================
# Dep. Variable:                   var1   R-squared:                       0.734
# Model:                            OLS   Adj. R-squared:                  0.706
# Method:                 Least Squares   F-statistic:                     26.21
# Date:                Sun, 05 Jul 2015   Prob (F-statistic):           3.45e-06