# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [default]
#     language: python
#     name: python3
# ---

# +
# %pylab
# %matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as st 

import statsmodels.api as sm

pd.set_option('display.mpl_style', 'default') 
pd.set_option('display.max_columns',1000)
plt.rcParams['figure.figsize'] = (10,6)

# -

Filename = "cities_r2.csv"
df = pd.read_csv(Filename)

df.head()

df.dtypes

df.shape


# ## Normality Tests

def Normality(data, name, ddof =1 , assume_discrete = False): 
    IQR = np.percentile(data, 75) - np.percentile(data, 25) 
    bins_width = 2 * IQR * (len(data)**(-1/3)) 
    suggest_bins = np.ceil((data.max() - data.min())/bins_width)
    
    if assume_discrete: 
        bins = arange(data.min(), data.max() + 2) - 0.5
        freq, bins, _ = plt.hist(data, bins = bins, normed = True , alpha=0.6) 
    else: 
        freq, bins, _ = plt.hist(data, bins= suggest_bins, normed = True, alpha = 0.6)
    X = st.norm(mean(data) , std(data, ddof = ddof))
    x1 = (bins[1:] + bins[:-1])/2
    y1 = X.pdf(x1) 
    plt.plot(x1,y1,color='red', marker = 'o')
    plt.title("Histogram")
    plt.xlabel(name) 
    plt.ylabel("Freq")
    plt.show() 
    
    x2 = sort(data) 
    probs = linspace(1, len(x2) , num = len(x2))/len(x2) 
    y2 = X.ppf(probs) 
    plt.plot(x2[:-1],y2[:-1],'o', color = 'blue') 
    plt.plot([x2.min() , x2.max()],[x2.min() , x2.max()],color= 'yellow',linewidth = 2) 
    plt.title("Line plot") 
    plt.xlabel(name) 
    plt.show() 
    
    print("Empirical Rule:") 
    x_bar, x_s = mean(data) , std(data, ddof = ddof) 
    for i in range(1,4): 
        within_std = np.sum(np.abs(data - x_bar) < i * x_s)
        print("Data within {} std: {:.2f}%".format(i, within_std/data.size * 100))


x  = df['effective_literacy_rate_total'].values

Normality(x,'effective_literacy_rate_total ')

y = df['sex_ratio'].values

Normality(y, 'sex Ratio')

# ## Linear Regression 

# +
x  = df['effective_literacy_rate_total'].values
X = sm.add_constant(x)
#Use X for determining coefficients and predicted values



y = df['sex_ratio'].values


# +
def fitted_model(X,y): 
    fitted_model = sm.OLS(y,X) 
    return fitted_model.fit() 

def coeff(X,y): 
    coeff = fitted_model(X,y).params
    return np.array(coeff)
    
def prediction(X,y): 
    n = len(X)
    predictions = fitted_model(X,y).predict(X) 
    return np.array(predictions.reshape(n,1))


# -

# Linear regression line 

print("y= {:.3f} + {:.3f}x".format(coeff(X,y)[0] , coeff(X,y)[1]))

# Prediction

prediction(X,y)[:10]


# ## Correlation coefficient 

# +
def cov(x,y,bias=0): 
    n = float(len(x))
    return np.sum((x - np.mean(x)) * (y - np.mean(y)))/(n - min(bias,1))

def z_score(x):
    return (x-np.mean(x))/np.std(x)

def correlation(x,y,bias = 0): 
    return cov(z_score(x),z_score(y),bias)


# -

print("Correlation coefficient = {:.3f}".format(correlation(x,y)))


# ## Determination Coefficient

def R_squared(x,y):
    n = len(x)
    X = np.mat([np.ones(n), x]).T
    Y = np.mat(y).T
    M = X.T * X

    coef1 = M.I * (X.T * Y)

    YR = X * coef1
    Ygem = Y.mean()

    SSE = np.dot((Y - YR).T, Y - YR)
    SSY = np.dot((Y - Ygem).T, Y - Ygem)
    SST = np.dot((YR - Ygem).T, YR - Ygem)

    Rkwadraat = 1 - SSE / SSY
    return(Rkwadraat[0,0])


print("Determination coefficient = {:.4f}".format(R_squared(x,y)))


# ## P-value

# +

def p_value(X,y):
    n = len(X) 
    R  = st.rdist(n - 2) 
    pvalue = 2*(1 - R.cdf(abs(correlation(x,y))))
    return pvalue


# -

print("pvalue = {:.10f} ".format(p_value(X,y)))


# ## Regression plot

def regression_plot(dataset,xname,yname,fitted_values): 
    x_range = [dataset[xname].min(), dataset[xname].max()]
    y_range = [dataset[yname].min(), dataset[yname].max()]
    
    scatter_plot = dataset.plot(kind='scatter', x = xname, y = yname, xlim=x_range,ylim=y_range, linewidth=10)
    
    y_hat = scatter_plot.plot(x_range,[dataset[yname].mean(), dataset[yname].mean()],'--', color='red', linewidth=1.5)
    x_hat = scatter_plot.plot([dataset[xname].mean(),dataset[xname].mean()],y_range, '--', color= 'red', linewidth=1.5) 
    
    regression_line = scatter_plot.plot(dataset[xname],fitted_values,color='yellow', linewidth=1.5)
    plt.title('Regression of {} and {}'.format(xname,yname))
    plt.show()


regression_plot(df,'effective_literacy_rate_total','sex_ratio',prediction(X,y) )

print("Correlation coefficient = {:.3f}".format(correlation(x,y)))
print("Coefficient of determination = {:.3f}".format(R_squared(x,y)))
print("pvalue = {:.10f} ".format(p_value(X,y)))

# ## Conclusion
#
#
# With n = 493 observations, the correlationcoefficient r = 0.256 with p-value = 0.0000000084. Effective literacy rate total and sex ratio correlate with each-other if assumed that 6.54% of the sex ratio are linearly dependent with Effective literacy rate total.

StarcoderdataJupyterScriptsDedupFiltered