def plotFit(fit):
  """Create's the 2x2 panel of plots that plot(fit) would create in R"""
  resid = fit.resid
  mu = resid.mean()
  std = resid.std(axis=0)

  #had to write my own normalize function
  def _normalize(resid):
      return (resid-mu)/std
  norm_resid = resid.apply(_normalize)


  f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex='col', sharey='row')

  ax1.scatter(fit.fittedvalues, fit.resid)
  ax1.set_xlabel('Fitted Values')
  ax1.set_ylabel('Residuals')
  ax1.set_title('Residuals vs Fitted')


  sm.qqplot(fit.resid, ax=ax2)
  ax2.set_title('QQ plot')

  ax3.scatter(fit.fittedvalues, norm_resid)
  ax3.set_xlabel('Fitted Values')
  ax3.set_ylabel('Standardized Residuals')
  ax3.set_title('Scale-Location')

  sm.graphics.influence_plot(fit, ax=ax4, criterion="cooks")

  plt.show()

def return_qqplot(data):
    ''' Generates a Q-Q plot of the returns.'''
    plt.figure(figsize=(9, 5))
    sm.qqplot(data['returns'], line='s')
    plt.grid(True)
    plt.xlabel('theoretical quantiles')
    plt.ylabel('sample quantiles')

def plot_single_peak(peak, ff = False, num_bins = 50, qq = scipy.stats.norm):
    '''Plotte fuer einen Peak das Histogramm sowie qq-Plot zur Verteilung qq
    Besser plot_simlist verwenden, wenn nicht nur gezielt ein Peak angeschaut werden soll, oder Histogrammdarstellung erwuenscht'''
    data = peak
    # Falls from_file gewaehlt, oeffne file
    if ff:
        with open (peak, 'rb') as daten:
            data = pickle.load(daten)
    #Normales Hist plotten
    n, bins, patches = plt.hist(data.times, num_bins, normed=1, alpha=0.5 )
    plt.suptitle("params:" + str(data.params))
    # Jetzt noch ein qq-Plot
    x = np.arange(1, 250, 0.5)
    if qq == scipy.stats.invgauss:
        mu, loc, scale =  scipy.stats.invgauss.fit(data.times)
        logging.log(20, "ig-paramss, %s, %s, %s", str(mu), str(loc), str(scale))
        plt.plot(x,scipy.stats.invgauss.pdf(x,mu, loc, scale))
        logging.log(20,'skew, %s', str(scipy.stats.skew(data.times)))
        sm.qqplot(np.array(data.times), qq, distargs=(mu,),  line = 'r')
        plt.suptitle("params:" + str(data.params) + " qq-Plot mit Normalverteilung" )
    elif qq == scipy.stats.norm:
        sm.qqplot(np.array(data.times), qq, line='r')
        plt.suptitle("params:" + str(data.params) + " qq-Plot mit Inverser Gauss Verteilung: ")
    else: 
        print("not yet implemented, distribution:", qq)
    plt.show()

def hist(request, sym):
    """create a histogram plot"""
    data = Data(syms=[sym], start=start)
        
    r = data.panel.r.copy()
    r = r.dropna()
    
    fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(12, 7))
    
    ax = axes[0,0]
    ax.hist(r[sym].values, bins=30)
    r.plot(kind="kde", ax=ax,grid=True)
    r.boxplot(ax=axes[0,1],grid=True)
    r.plot(kind="kde", ax=axes[1,0],grid=True)
    sm.qqplot(r[sym], line='r', fit=True, ax=axes[1,1])
    
    r['mean'] = pandas.rolling_mean(r[sym], 12)
    r['std'] = pandas.rolling_std(r[sym], 12)
    r['cum_ret'] = r[sym].cumsum()
    r[['mean', 'std']].plot(ax=axes[0,2], grid=True, rot=45)
    
    r[['cum_ret']].plot(ax=axes[1,2], grid=True, rot=45)
    
    fig.tight_layout()
    
    fig.set_facecolor((1,.8,.6,0))
    canvas = FigureCanvas(fig)
    response = HttpResponse(content_type='image/png')
    canvas.print_png(response)
    return response

def plot (sim_liste, histogram_separate, histogram_spec, qq_Plot, fit_qq_Plot, num_bins = 50, vergleich= scipy.stats.invgauss):
    startzeit = time.clock()   
    if histogram_spec:
        print "Erstelle Spektrum"
        fig, ax = plt.subplots()
        fig.suptitle("Laenge: "+str(sim_liste[0].length)+" Anz Teilchen: " +str(sim_liste[1].number)) #TODO, gehe hier davon aus, dass gleiche sim-bedingungen vorliegen
        for sim in sim_liste:
            ax.hist(sim.times, num_bins, alpha=0.5, normed = 1, label = str(sim.params) )
       # plt.show()  
        legend = ax.legend(loc='upper right', shadow=True)

    
    # Je Simulation ein Ausgabefenster mit separatem Histogramm/qq-Plot mit gewählten Params/qq mit automatischem Fit 
    number_stats = sum([histogram_separate, qq_Plot, fit_qq_Plot])
    print number_stats
    if histogram_separate or qq_Plot or fit_qq_Plot:
	print "Erstelle separate Dinge"
	for sim in sim_liste:
	    fig = plt.figure(figsize=(4*number_stats, 4))
            gs1 = gridspec.GridSpec(1, number_stats)
            ax_list = [fig.add_subplot(ss) for ss in gs1]
           
	    akt = 0
	    fig.suptitle("ps, pm"+str(sim.params)+str(round(sim.params[0]-sim.params[1],5)), size = 15)
	    if histogram_separate:
		ax_list[akt].hist(sim.times, num_bins)
		ax_list[akt].set_title("Histogramm")
                akt+=1
                
            #print "hist sep", time.clock()-startzeit
	    if qq_Plot:
                sm.qqplot (np.array(sim.times), scipy.stats.norm,  line = 'r', ax=ax_list[akt])
		ax_list[akt].set_title("qq-Plot; norm!! Params: 0.05")
                akt+=1
            #print 'qq 0.05', time.clock()-startzeit
	    if fit_qq_Plot:
		                
                #mu, loc, scale = scipy.stats.invgauss.fit(sim.times)
                #mean, var = scipy.stats.invgauss.stats(mu, loc, scale, moments='mv')
                #print  "params", sim.params, '(mu, loc, scale), mean, var', round(mu, 5), round(loc, 2), round(scale, 2), '\n',  mean, '\n', var
                
                #sm.qqplot (np.array(sim.times), vergleich, fit = True,  line = 'r', ax=ax_list[akt])
		#ax_list[akt].set_title("qq-Plot mit auto Fit")
                #akt+=1 
                sm.qqplot (np.array(sim.times), vergleich, distargs= (sim.mu, ),  line = 'r', ax=ax_list[akt])
		ax_list[akt].set_title("qq-Plot mit mu:" + str(sim.mu))
                akt+=1
            #print "qq plus rechnen", time.clock()-startzeit                

                #fig.subplots_adjust(top=5.85)
            gs1.tight_layout(fig, rect=[0, 0.03, 1, 0.95]) 
            print time.clock()-startzeit
            #plt.tight_layout()
    plt.show()    
       

   
   
    '''x = np.linspace(0, 2*np.pi, 400)

def plot_single_histqq_ff(datei, num_bins=50):
    with open(datei, 'rb') as daten:
        sim = pickle.load(daten)
        n, bins, patches = plt.hist(sim.times, num_bins, normed=1, alpha=0.5 )
        x = np.arange(50000, 250000, 100)
        print "ig-params", scipy.stats.invgauss.fit(sim.times)
        mu, loc, scale =  scipy.stats.invgauss.fit(sim.times)
        plt.plot(x,scipy.stats.invgauss.pdf(x,mu, loc, scale))
        print 'skew', scipy.stats.skew(sim.times)
        
        sm.qqplot(np.array(sim.times), scipy.stats.invgauss, distargs=(mu,),  line = 'r')

def plot_qq_checkout():

	path = './qq_checkout'
	if os.path.exists(path) == False:
		os.mkdir(path)

	global number_attribute_remove_lost_arr
	for k, v in number_attribute_remove_lost_arr.iteritems():
		sm.qqplot(np.array(v), line='r')
		#plt.xlabel(k)
		plt.title(k)
		plt.grid(True)
		#plt.show()
		plt.savefig(path + '/' + k + '.png')
		plt.close()

def plot_ic_qq(ic, theoretical_dist=stats.norm, ax=None):
    """
    Plots Spearman Rank Information Coefficient "Q-Q" plot relative to
    a theoretical distribution.

    Parameters
    ----------
    ic : pd.DataFrame
        DataFrame indexed by date, with IC for each forward return.
    theoretical_dist : scipy.stats._continuous_distns
        Continuous distribution generator. scipy.stats.norm and
        scipy.stats.t are popular options.
    ax : matplotlib.Axes, optional
        Axes upon which to plot.

    Returns
    -------
    ax : matplotlib.Axes
        The axes that were plotted on.
    """

    ic = ic.copy()

    num_plots = len(ic.columns)

    v_spaces = ((num_plots - 1) // 3) + 1

    if ax is None:
        f, ax = plt.subplots(v_spaces, 3, figsize=(18, v_spaces * 6))
        ax = ax.flatten()

    if isinstance(theoretical_dist, stats.norm.__class__):
        dist_name = 'Normal'
    elif isinstance(theoretical_dist, stats.t.__class__):
        dist_name = 'T'
    else:
        dist_name = 'Theoretical'

    for a, (period_num, ic) in zip(ax, ic.iteritems()):
        sm.qqplot(ic.replace(np.nan, 0.).values, theoretical_dist, fit=True,
                  line='45', ax=a)
        a.set(title="{} Period IC {} Dist. Q-Q".format(
              period_num, dist_name),
              ylabel='Observed Quantile',
              xlabel='{} Distribution Quantile'.format(dist_name))

    return ax

def plot_model(prediction, y, x):
    fig, axs = sns.plt.subplots(2, 2, figsize=(16, 10))
    axs = axs.flatten()

    resid = pd.Series(y - prediction, index=y.index, name='Residuals')

    resid.hist(bins=40, ax=axs[0])
    axs[0].set_xlabel('Residuals')
    sm.qqplot(resid, line='q', ax=axs[1])
    axs[1].set_xlabel('Residuals')
    tbpd.hist2d(resid, prediction, ax=axs[2],
                vlabel='Residuals', hlabel='Predicted value',
                integer_aligned_bins=True)
    tbpd.hist2d(y, prediction, ax=axs[3],
                vlabel='True value', hlabel='Predicted value',
                integer_aligned_bins=True, sqrt=True)
    fig.tight_layout()

def q_q_plot(filepath, parameter):
    df = pandas.read_csv(filepath)
    array = df[parameter]
    try:
        fig = sm.qqplot(array, scipy.stats.t, fit=True, line='45')
        plt.show()
    except:
        print "There was an error."

def test_qqplot():
    #just test that it runs
    data = sm.datasets.longley.load()
    data.exog = sm.add_constant(data.exog, prepend=False)
    mod_fit = sm.OLS(data.endog, data.exog).fit()
    res = mod_fit.resid
    fig = sm.qqplot(res, line='r')

    plt.close('all')

def test_qqplot():
    #just test that it runs
    data = sm.datasets.longley.load()
    data.exog = sm.add_constant(data.exog)
    mod_fit = sm.OLS(data.endog, data.exog).fit()
    res = mod_fit.resid
    fig = sm.qqplot(res)

    plt.close(fig)

def plot(file_name,negative_control_gRNAs=None,wald_only=False):
    data=open(file_name,'rb')
    short_file_name=file_name[:file_name.index(".gene_summary.txt")]
    data.readline()
    permute_p_value_list=[]
    wald_p_value_list=[]
    beta_value_list=[]

    if negative_control_gRNAs!=None:
        negative_control_permute_p_value_list=[]
        negative_control_wald_p_value_list=[]
        negative_control_beta_value_list=[]


    for line in data:
        elements=line.decode().strip().split("\t")
        if negative_control_gRNAs!=None and elements[0] in negative_control_gRNAs:
            negative_control_beta_value_list.append(float(elements[2]))
            if wald_only==True:
                negative_control_wald_p_value_list.append(float(elements[4]))
            else:
                negative_control_permute_p_value_list.append(float(elements[4]))
                negative_control_wald_p_value_list.append(float(elements[6]))
        else:
            beta_value_list.append(float(elements[2]))
            if wald_only==True:
                wald_p_value_list.append(float(elements[4]))
            else:
                permute_p_value_list.append(float(elements[4]))
                wald_p_value_list.append(float(elements[6]))
    beta_value_list=[x for x in beta_value_list if str(x) != 'nan' and abs(x)<3]
    wald_p_value_list=[x for x in wald_p_value_list if str(x) != 'nan']
    if negative_control_gRNAs!=None:
        negative_control_beta_value_list=[x for x in beta_value_list if str(x) != 'nan' and abs(x)<3]
        negative_control_wald_p_value_list=[x for x in wald_p_value_list if str(x) != 'nan']

    if wald_only!=True:
        permute_p_value_list=[x for x in permute_p_value_list if str(x) != 'nan']
        stats.probplot(permute_p_value_list, dist="uniform",plot=pylab)
        pylab.savefig("QQplot of permute_p value %s.png" %short_file_name)
        pylab.close()

    pylab.hist(beta_value_list,bins=1000)
    pylab.savefig("Hist of beta value %s.png" %short_file_name)
    pylab.close()

    #stats.probplot(wald_p_value_list, dist="uniform",plot=pylab)
    fig=sm.qqplot(np.array(wald_p_value_list),stats.uniform,fit=True, line='45')
    pylab.xlim(0,1)
    pylab.ylim(0,1)
    #fig.set_xlim(0,1)
    pylab.savefig("QQplot of wald_p value %s.png" %short_file_name)
    pylab.close()
    '''

 def qqPlot(self):
     """ Plots sample signals against theorethical distribution"""
     import statsmodels.api as sm #pandas, patsy
     import matplotlib.pyplot as plt
     data = self.array.probes[:, 2 + self.number]  # add log2
     plt.figure(self.number)
     fig = sm.qqplot(data)
     plt.xlabel('Theoretical quantiles')
     plt.ylabel('Sample quantiles')
     plt.title('Probe intensities for %s' % (self.name))
     plt.savefig("%s_qqprob.png" % (self.name))

    def tsplot(y, lags=None, figsize=(10, 8), style='bmh'):
        if not isinstance(y, pd.Series):
            y = pd.Series(y)
        with plt.style.context(style):    
            fig = plt.figure(figsize=figsize)
            #mpl.rcParams['font.family'] = 'Ubuntu Mono'
            layout = (3, 2)
            ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
            acf_ax = plt.subplot2grid(layout, (1, 0))
            pacf_ax = plt.subplot2grid(layout, (1, 1))
            qq_ax = plt.subplot2grid(layout, (2, 0))
            pp_ax = plt.subplot2grid(layout, (2, 1))

            y.plot(ax=ts_ax)
            ts_ax.set_title('Time Series Analysis Plots')
            smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)
            smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)
            sm.qqplot(y, line='s', ax=qq_ax)
            qq_ax.set_title('QQ Plot')        
            scs.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)

            plt.tight_layout()
        return

def print_qqplot_and_residuals_plot(model):
    # qq-plot
    ax1 = plt.subplot(1, 3, 1)
    qq_plot = sm.qqplot(model.resid, line = 'r', ax = ax1)
    
    # Residuals plot
    ax2 = plt.subplot(1, 3, 2)
    stdres = pandas.DataFrame(model.resid_pearson)
    residuals_plot = plt.plot(stdres, 'o', ls = 'None')
    plt.axhline(y = 0, color = 'r')
    plt.ylabel('Standarized Residual')
    plt.xlabel('Observation Number')
    
    plt.show()

def mult_regression(wine_set):
    # center quantitative IVs for regression analysis
    w = wine_set['quality']
    wine_set = wine_set - wine_set.mean()
    wine_set['quality'] = w

    print ("OLS multivariate regression model")
    # first i have run with all columns; than chose the most significant for each wine set and rerun:

    if len(wine_set) < 2000:
        # for red
        model1 = smf.ols(
            formula="quality ~ volatile_acidity + chlorides + pH + sulphates + alcohol",
            data=wine_set)
    else:
        # for white
        model1 = smf.ols(
            formula="quality ~ volatile_acidity + density + pH + sulphates + alcohol",
            data=wine_set)

    results1 = model1.fit()
    print(results1.summary())

    # q-q plot for normality
    qq = sm.qqplot(results1.resid, line = 'r')
    plt.show()

    # plot of residuals
    stdres = pd.DataFrame(results1.resid_pearson)
    plt.plot(stdres, 'o', ls = 'None')
    l = plt.axhline(y=0, color = 'r')
    plt.ylabel('Standardized redisual')
    plt.xlabel('Observation number')
    plt.show()

    # # diagnostic plots
    # figure1 = plt.figure(figsize=(12, 8))
    # figure1 = sm.graphics.plot_regress_exog(results1, "alcohol", fig = figure1)
    # plt.show()
    #
    # figure1 = plt.figure(figsize=(12, 8))
    # figure1 = sm.graphics.plot_regress_exog(results1, "sulphates", fig = figure1)
    # plt.show()

    # leverage plot
    figure1 = sm.graphics.influence_plot(results1, size=8)
    plt.show()

def plot_box_resids(fit_model, y_pred, subset=None):
    '''More than you ever wanted to know about your residuals'''
    s_resid = (fit_model.resid - np.mean(fit_model.resid)) /\
               np.var(fit_model.resid)
    if subset:
        s_resid = np.random.choice(s_resid,
                                  replace=False,
                                  size=math.floor(len(s_resid) * subset))
    df = pd.DataFrame(s_resid, columns=['resids'])
    temp_df = pd.DataFrame(y_pred, columns=['target'])
    df = df.join(temp_df)

    if min(y_pred) < -1:
        df['turnout_bucket'] = df['target']\
        .apply(lambda x: int(math.floor(10 * np.exp(x))))
        y = df['target'].apply(lambda x: np.exp(x))
    else:
        df['turnout_bucket'] = df['target']\
        .apply(lambda x: int(math.floor(10 * x)))
        y = df['target']

    posit = sorted(df['turnout_bucket'].unique())

    plt.scatter(y, s_resid, alpha=.2)
    slope, intercept = np.polyfit(y, s_resid, 1)
    plt.plot(y, np.poly1d(np.polyfit(y, s_resid, 1))(y))
    plt.title('Studentized Residuals vs Prediction')
    plt.xlabel('Predicted Value')
    plt.ylabel('Studentized Residual')
    print 'Slope of best fit line: %s' % slope
    plt.show()

    ax1 = df[['resids', 'turnout_bucket']]\
        .boxplot(by='turnout_bucket', positions=posit, widths=.5)
    plt.title('Residuals versus Turnout')
    plt.xlabel('Turnout Bucket')
    plt.ylabel('Studentized Residuals')
    plt.suptitle('')
    plt.show()

    fig = sm.qqplot(s_resid, line='s')
    plt.title('Q-Q Plot')
    plt.show()

    w, p_val = shapiro(s_resid)
    print 'Shapiro-Wilk P_val is %s, larger the better' % p_val

    k, p_val = normaltest(s_resid)
    print 'D’Agostino and Pearson’s P_val is %s, larger the better' % p_val

    k, p_val = kstest(s_resid, 'norm')
    print 'Kolmogorov–Smirnov P_val is %s, larger the better' % p_val

    A, critical, sig = anderson(s_resid)
    print 'Anderson-Darling A2 is %s, smaller the better' % A
    print critical
    print sig

    n, bins, patches = plt.hist(s_resid, 75, normed=1)
    mu = np.mean(s_resid)
    sigma = np.std(s_resid)
    plt.plot(bins, mlab.normpdf(bins, mu, sigma))
    plt.title('Residuals versus a Normal Dist')
    plt.show()

    df['turnout_bucket'].hist(bins=posit, align='left', color='b')
    plt.title('Histogram of Turnout Bucket')
    plt.ylabel('Count')
    plt.xlim(-.5, - .5 + len(posit))

    temp = df[['resids', 'turnout_bucket']].groupby('turnout_bucket').count()
    temp.columns = ['Count']
    plt.show()
    print temp

reg2 = smf.ols('lifeexpectancy ~ breastcancerper100th_c + I(breastcancerper100th_c**2)', data=sub1).fit()
print (reg2.summary())


####################################################################################
# EVALUATING MODEL FIT
####################################################################################

# adding alcohol consumption
reg3 = smf.ols('lifeexpectancy ~ breastcancerper100th_c + I(breastcancerper100th_c**2) + breastcancerper100th_c', 
               data=sub1).fit()
print (reg3.summary())


#Q-Q plot for normality
fig4=sm.qqplot(reg3.resid, line='r')

# simple plot of residuals
stdres=pandas.DataFrame(reg3.resid_pearson)
plt.plot(stdres, 'o', ls='None')
l = plt.axhline(y=0, color='r')
plt.ylabel('Standardized Residual')
plt.xlabel('Observation Number')


# additional regression diagnostic plots
fig2 = plt.figure(figsize=(12,8))
fig2 = sm.graphics.plot_regress_exog(reg3,  "breastcancerper100th_c", fig=fig2)

# leverage plot
fig3=sm.graphics.influence_plot(reg3, size=8)

comb.boxplot(column=[0])


## Q-Q Plot 

##### In statistics, a Q–Q plot ("Q" stands for quantile) is a probability plot, which is a graphical method for comparing two probability distributions by plotting their quantiles against each other. If the two distributions being compared are similar, the points in the Q–Q plot will approximately lie on the line y = x. If the distributions are linearly related, the points in the Q–Q plot will approximately lie on a line, but not necessarily on the line y = x.

# In[266]:

import statsmodels.api as sm


# In[269]:

sm.qqplot(comb[1],line='45')


# In[275]:

os.getcwd()


# In[287]:

for i in np.arange(0,40,1):
    pieces1='histograms/histogram',format(i),'.jpg'
    hist=comb[i].hist()
    fig = hist.get_figure()
    fig.savefig(''.join(pieces1))
    fig.clear()