%load_ext watermark

%load_ext lab_black

%matplotlib inline

import pandas as pd
import numpy as np
import seaborn as sn

import math

import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.formula.api import rlm
import statsmodels.api as sm

from statsmodels.sandbox.regression.predstd import (
    wls_prediction_std,
)

from statsmodels.stats.stattools import jarque_bera

DATA_PREFIX = 'd:/StanfordStats191/alr3data/'

data = pd.read_csv(DATA_PREFIX + 'heights.txt', sep=' ')

data.head()

data.describe()

data['Mheight'].hist()

data['Dheight'].hist()

data.boxplot(column=['Mheight', 'Dheight'])

__ = pd.plotting.scatter_matrix(
    data[['Mheight', 'Dheight']]
)

fig, ax = plt.subplots(figsize=(10, 10))

ax.plot(data['Mheight'], data['Dheight'], 'o', label='data')
ax.grid()
ax.set_xlabel("Mother's Height")
ax.set_ylabel("Daughter's Height")

ax.legend(loc='best')

fig, ax = plt.subplots(figsize=(10, 10))

sn.regplot(
    'Mheight',
    'Dheight',
    data=data,
    ax=ax,
    label='Seaborn Regression\n95% Conf Intvl',
)
ax.grid()
ax.set_xlabel("Mother's Height")
ax.set_ylabel("Daughter's Height")

ax.legend(loc='best')

# fig, ax = plt.subplots(figsize=(10, 10))

g = sn.jointplot(
    'Mheight',
    'Dheight',
    kind='reg',
    data=data,
    # ax=ax,
    label='Seaborn Regression\n95% Conf Intvl',
)

g = sn.jointplot(
    'Mheight',
    'Dheight',
    kind='hex',
    data=data,
    # ax=ax,
    label='Seaborn Regression\n95% Conf Intvl',
)

data = pd.read_csv(DATA_PREFIX + 'heights.txt', sep=' ')

daughter = data['Dheight']
mother = data['Mheight']

mother = sm.add_constant(mother)

res = sm.OLS(daughter, mother).fit()

res.summary()

fig, ax = plt.subplots(figsize=(10, 10))

ax.plot(
    data['Mheight'],
    data['Dheight'],
    'o',
    label='data',
    alpha=0.5,
)
ax.plot(
    data['Mheight'],
    res.fittedvalues,
    'r-',
    label='Fitted Values',
)
ax.grid()
ax.set_xlabel("Mother's Height")
ax.set_ylabel("Daughter's Height")

ax.legend(loc='best')

x_ci = np.linspace(
    data['Mheight'].min(), data['Mheight'].max(), 20
)
x_ci_enh = sm.add_constant(x_ci)

gp = res.get_prediction(x_ci_enh)
ci = gp.conf_int(alpha=0.05)
upper = [z[1] for z in ci]
lower = [z[0] for z in ci]

prstd, iv_l, iv_u = wls_prediction_std(res)

fig, ax = plt.subplots(figsize=(10, 10))

ax.plot(
    data['Mheight'],
    data['Dheight'],
    'o',
    label='data',
    alpha=0.5,
)
ax.plot(
    data['Mheight'],
    res.fittedvalues,
    'g-',
    label='Fitted Values',
)

ax.plot(
    x_ci, upper, 'r-', label='95% Conf Int (u)', alpha=0.5
)
ax.plot(
    x_ci, lower, 'r-', label='95% Conf Int (l)', alpha=0.5
)

ax.plot(
    data['Mheight'], iv_l, 'y-', label='95% Pred Int (l)'
)
ax.plot(
    data['Mheight'], iv_u, 'y-', label='95% Pred Int (u)'
)

ax.grid()
ax.set_xlabel("Mother's Height")
ax.set_ylabel("Daughter's Height")

ax.legend(loc='best')

res2 = ols('Dheight ~ Mheight', data=data).fit()

res2.summary()

res.summary()

mean_y = data['Dheight'].mean()

sstot = sum([(y - mean_y) ** 2 for y in data['Dheight']])

print(f'Total sum of squares around the mean {sstot}')

ssreg = sum([(y1 - mean_y) ** 2 for y1 in res.fittedvalues])

ssres = sum(
    [
        (y1 - y2) ** 2
        for y1, y2 in zip(data['Dheight'], res.fittedvalues)
    ]
)

print(
    f'Total sum of square spread around mean {sstot},\n'
    + f'sum of square errors, predictions around mean {ssreg}, \n'
    + f'sum of square errors, predictions against actual {ssres}'
)

print(ssreg + ssres)

r2 = 1 - ssres / sstot
r2

print(
    f'res.ssr = {res.ssr}; res.mse_model = {res.mse_model}'
)

print(
    f'res.mse_total = {res.mse_total}, se_total = {res.mse_total*(1375-1)}'
)

print(res.centered_tss)

print(
    f'ssres {ssres} =  res.mse_resid {res.mse_resid} * 1373 = {res.mse_resid*1373}'
)

DATA_PREFIX = 'd:/StanfordStats191/'

rtw = pd.read_csv(DATA_PREFIX + 'P005.txt', sep='\t')

rtw.head()

fig, ax = plt.subplots(figsize=(10, 10))

props = dict(boxes="Pink")
rtw.boxplot(
    column=['COL'],
    by='RTWL',
    ax=ax,
    color=props,
    patch_artist=True,
)
ax.set_ylabel("Cost Of Living")

__ = pd.plotting.scatter_matrix(
    rtw[
        [
            'City',
            'COL',
            'PD',
            'URate',
            'Pop',
            'Taxes',
            'Income',
            'RTWL',
        ]
    ]
)

rtw.describe()

rtw.median()

fig, ax = plt.subplots(figsize=(10, 10))

__ = pd.plotting.scatter_matrix(
    rtw[['COL', 'PD', 'URate', 'Pop', 'Taxes', 'Income']],
    c=rtw['RTWL'],
    cmap='rainbow',
    ax=ax,
    alpha=0.95,
)

_ = sn.pairplot(
    rtw,
    vars=['COL', 'PD', 'URate', 'Pop', 'Taxes', 'Income'],
    kind='scatter',
    hue='RTWL',
    diag_kind='hist',
    diag_kws={'alpha': 0.5},
)

rtw.iloc[26]

DATA_PREFIX = 'd:/StanfordStats191/'

cal = pd.read_csv(DATA_PREFIX + 'calcium.txt', sep='\t')

cal.head()

cal.describe()

cal[cal['Treatment'] == 'Calcium']['Decrease'].mean()

cal[cal['Treatment'] == 'Placebo']['Decrease'].mean()

fig, ax = plt.subplots(figsize=(10, 10))
cal[cal['Treatment'] == 'Calcium']['Decrease'].hist(
    ax=ax, label='Calcium'
)
cal[cal['Treatment'] == 'Placebo']['Decrease'].hist(
    ax=ax, label='Placebo', alpha=0.5
)
ax.legend(loc='best')
ax.set_ylabel('Count')
ax.set_xlabel('Decrease in blood pressure')

fig, ax = plt.subplots(figsize=(6, 6))
cal.boxplot(column=['Decrease'], by='Treatment', ax=ax)

mean_decrease = cal['Decrease'].mean()
mean_treated = cal[cal['Treatment'] == 'Calcium'][
    'Decrease'
].mean()
mean_placebo = cal[cal['Treatment'] == 'Placebo'][
    'Decrease'
].mean()

print(mean_decrease, mean_placebo, mean_treated)

median_decrease = cal['Decrease'].mean()
median_decrease

median_treated = cal[cal['Treatment'] == 'Calcium'][
    'Decrease'
].median()
median_treated

sorted(cal[cal['Treatment'] == 'Calcium']['Decrease'])

cal[cal['Treatment'] == 'Calcium']['Decrease'].quantile(
    q=[0.25, 0.75]
)

treated = cal[cal['Treatment'] == 'Calcium']['Decrease']
samp_std_dev = sum(
    [(x - mean_treated) ** 2 for x in treated]
) / (len(treated) - 1)

math.sqrt(samp_std_dev)

cal[cal['Treatment'] == 'Calcium']['Decrease'].std()

sampling = [
    sum(
        np.random.choice(
            cal['Decrease'], size=5, replace=True
        )
    )
    / 5
    for i in range(200)
]

sampling_treated = [
    sum(
        np.random.choice(
            cal[cal['Treatment'] == 'Calcium']['Decrease'],
            size=5,
            replace=True,
        )
    )
    / 5
    for i in range(200)
]
sampling_placebo = [
    sum(
        np.random.choice(
            cal[cal['Treatment'] == 'Placebo']['Decrease'],
            size=5,
            replace=True,
        )
    )
    / 5
    for i in range(200)
]

fig, ax = plt.subplots(figsize=(6, 6))

ax.hist(sampling, label='Combined', color='blue', alpha=0.5)


ax.hist(
    sampling_placebo,
    alpha=0.5,
    label='Placebo Only',
    color='green',
)
ax.hist(
    sampling_treated,
    alpha=0.5,
    label='Calcium Only',
    color='red',
)

ax.axvline(mean_decrease, c='blue', label='Total')
ax.axvline(mean_treated, c='red', label='Treated')
ax.axvline(mean_placebo, c='green', label='Placebo')

ax.legend(loc='best')

sample = [
    5,
    -1,
    0,
    2,
    6,
    -1,
    5,
    6,
    3,
    3,
    -3,
    -1,
    -2,
    3,
    0,
    4,
    1,
    4,
    6,
    0,
]
print(np.mean(sample), np.std(sample, ddof=1))

t_value = (np.mean(sample) - 0) / (
    np.std(sample, ddof=1) / math.sqrt(20)
)

t_value

len(sample)

stats.ttest_1samp(sample, 0)

fig, ax = plt.subplots(figsize=(6, 6))

ax.hist(sample, label='Data')
ax.axvline(0, c='blue', label='H0: Mean=0')
ax.axvline(np.mean(sample), c='red', label='Actual Mean')
ax.legend(loc='best')

df = 19

crit1, crit2 = stats.t.ppf([1 - 0.01, 0.01], df)
print(crit1, crit2)

x_rej1 = np.linspace(stats.t.ppf(0.001, df), crit2, 20)
x_rej2 = np.linspace(crit1, stats.t.ppf(0.999, df), 20)

fig, ax = plt.subplots(figsize=(6, 6))

x = np.linspace(
    stats.t.ppf(0.001, df), stats.t.ppf(0.999, df), 50
)

ax.plot(x, stats.t.pdf(x, df), 'r-', label=f'T (df={df})')


ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

ax.axvline(t_value, c='black', label='T value')
ax.axvline(crit1, c='yellow', label='99% range')
ax.axvline(crit2, c='yellow', label='99% range')

ax.fill_between(
    x_rej1,
    0,
    stats.t.pdf(x_rej1, df),
    facecolor='red',
    label='Rejection Region (99%)',
    alpha=0.5,
)
ax.fill_between(
    x_rej2,
    0,
    stats.t.pdf(x_rej2, df),
    facecolor='red',
    alpha=0.5,
)
ax.grid()

ax.legend(loc='right')

stats.t.ppf([0.975, 0.025], df)

treatment = cal[cal['Treatment'] == 'Calcium']['Decrease']
placebo = cal[cal['Treatment'] == 'Placebo']['Decrease']
l_treatment = len(treatment)
l_placebo = len(placebo)

print(
    f'Number of treatments: {l_treatment}, number of placebos: {l_placebo}'
)

t_mean = treatment.mean()
p_mean = placebo.mean()

t_std = np.std(treatment, ddof=1)
p_std = np.std(placebo, ddof=1)

ss2 = (l_treatment - 1) * t_std ** 2 + (
    l_placebo - 1
) * p_std ** 2

mean_ss2 = ss2 / (l_treatment + l_placebo - 2)

t_value = (t_mean - p_mean) / (
    math.sqrt(mean_ss2)
    * math.sqrt(1 / l_treatment + 1 / l_placebo)
)

print(t_mean, p_mean, t_value)

stats.ttest_ind(treatment, placebo, equal_var=True)

df = 19
crit1, crit2 = stats.t.ppf([1 - 0.025, 0.025], df)
print(crit1, crit2)

fig, ax = plt.subplots(figsize=(6, 6))

df = 19
x = np.linspace(
    stats.t.ppf(0.001, df), stats.t.ppf(0.999, df), 50
)

ax.plot(x, stats.t.pdf(x, df), 'r-', label=f'T (df={df})')


ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

ax.axvline(t_value, c='black', label='T value')
ax.axvline(crit1, c='green', label='95% upper cutoff value')
ax.legend(loc='upper left')

stats.ttest_ind(treatment, placebo, equal_var=False)

res = ols(formula='Decrease ~ Treatment', data=cal).fit()
print(res.summary())

x_index = np.linspace(
    1, len(res.fittedvalues) + 1, len(res.fittedvalues)
)
fig, ax = plt.subplots(figsize=(6, 6))

ax.plot(x_index, cal['Decrease'], 'ro', label='Data')
ax.plot(
    x_index, res.fittedvalues, 'g+', label='Fitted values'
)

# we many residual lines, but only label 1 for the legend
# do all residual error lines
for x, y1, y2 in zip(
    x_index, cal['Decrease'], res.fittedvalues
):
    ax.plot([x, x], [y1, y2], 'y-')
# end for

# label the first error line, then bail
for x, y1, y2 in zip(
    x_index, cal['Decrease'], res.fittedvalues
):
    ax.plot(
        [x, x], [y1, y2], 'y-', label='Prediction Error'
    )
    break
# end for

ax.grid()
ax.set_xlabel('Sample Number')
ax.set_ylabel('BP Decrease')
ax.legend(loc='best')

%watermark -h -iv
%watermark