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

# # Anscombe's Quartet Tutorial
#
# Author: <NAME>
#
# ## Tutorial goals:
#
# * explore some simple statistics
# * illustrate the importance of visualization
# * robust vs non-robust statistics
# * create and use pandas data frames
# * derive and fit simple linear regression
# * the role of outliers
# * explore bayesian approach

# ## Background
#
# The Anscombe data set shows multiple groups of data which have similiar summary statistics, **but** as we will see, there are some differences in the simple statistics, which will made clear by the distinction of robust and non-robust statistics. 
#
# This data set was created by the statistician <NAME> in 1973 to demonstrate both the importance of graphing data before analyzing it and the effect of outliers on statistical properties.
#
# There is also an interesting paper by Chattergee and Firat on how to derive arbitrary amounts of data that fit these criteria:
#
# http://www.tandfonline.com/doi/abs/10.1198/000313007X220057#.VjDpsZeUirM
#

# ### Load libraries

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

# %matplotlib inline

# ### Load the Seaborn Anscombe data set

anscombe = sns.load_dataset("anscombe")
#anscombe
#anscombe.describe() # does not work need to break up

sns.lmplot(x="x", y="y", col="dataset", data=anscombe, aspect=.5);

# ### Visualization
#
# Plotting the data allows one to see that while each quartet is fit by the same linear function, the grouping and clustering of the data is different in each case. This can be quickly seen by plotting.
#
# Recall:
#
# * a picture is worth a thousand words

# ## Build up our own Pandas data frame
#
# start with numpy arrays

# +
x1 = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5])
x2 = np.array([10,8,13,9,11,14,6,4,12,7,5])
x3 = np.array([10,8,13,9,11,14,6,4,12,7,5])
x4 = np.array([8, 8, 8,8,8,8,8,19,8,8,8])

y1 = np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])
y2 = np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])
y3 = np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])
x4 = np.array([8,8,8,8,8,8,8,19,8,8,8])
y4 = np.array([6.58,5.76,7.71,8.84,8.47,7.04,5.25,12.50,5.56,7.91,6.89])

a_xs = [x1,x2,x3,x4]
a_ys = [y1,y2,y3,y4]
# -

sns.regplot(x=x1, y=y1)

# ### Load local pdapt libraries

import os
os.chdir('/home/matej/develop/pdapt')
os.getcwd()

# import my statistical library (pdapt stats)
import pdapt_lib.machine_learning.stats as pstats

pstats.standard_deviation(x1)

# ## Summary statistics
#
# Notice that the robust statistics like the median, interquartile range (IQR), etc, **do** vary while the non-robust statisics: mean and standard deviation do not and very similar for each x and y data set in the quartet.

# summary stats for each set should be same eg compare means
for x in a_xs:
    pstats.summary(x)
    print()

for y in a_ys:
    pstats.summary(y)
    print()

# This value shows a close grouping in the y ranges of the data
pstats.interquartile_range(y3)

# This value shows a close grouping in the x ranges of the data
pstats.interquartile_range(x4)

# print out y4 and standarized equivalent..
# note standardized value of 12.5 is above abs(2) 
#   which is another definition of and an outlier
s = list(zip(y4, pstats.standardize(y4)))
sorted(s,key=lambda x: abs(x[1]))


def winsorization(v, edge):
    """
    input: v (vector), edge (percentile of vector to replace with quantile value)
    output winsorized vector
    """
    vs = sorted(v)
    lower_value_limit = pstats.quantile(v, edge)
    upper_value_limit = pstats.quantile(v, 1.0 - edge)
    lower_value = vs[vs.index(lower_value_limit)+1]
    upper_value = vs[vs.index(upper_value_limit)]
    w = []
    for i in v:
        if i < lower_value:
            w.append(lower_value)
        elif i > upper_value:
            w.append(upper_value)
        else:
            w.append(i)
    return w


a = [92, 19, 101, 58, 1053, 91, 26, 78, 10, 13, -40, 101, 86, 85, 15, 89, 89, 28, -5, 41]

# should replace -40 with -5 and 1053 with 101
w = winsorization(a,0.05)
w

print(pstats.mean(a))
print(pstats.mean(w))

# ### Outliers
#
# Some notes:
#
# * The 3rd and 4th parts of the dataset (quartet) have outliers which skew the mean and standard deviation to be equal to the other cases.
#
# * Especially in cases of small sample sizes, it is not necessarily true that the outliers should be removed, they may be part of the naturally occuring underlying distribution of the population of interest. On the other hand, the outliers could be due to instrument failure, or the result of someone answering the wrong question.
#
# * *one must investigate the outliers!*

# ## Linear Regression
#
# Here we will solve for the coefficients a few different ways, using:
#
# 1. general solution via linear algebra
# 2. a statistical relation
# 3. analytical solution to simple regression 

# ### 1. General Solution
# To fit a linear regression
#
# https://en.wikipedia.org/wiki/Linear_regression
#
# we need to solve a general linear equation, 
#
#  $\mathbf{Ax} = \mathbf{b}$
#  
# https://en.wikipedia.org/wiki/System_of_linear_equations
#  
# but here we are solving for the coefficients (weights) so looks more like
#
# $\mathbf{A w}=\mathbf{y}$
#
#
#  

# +
from numpy import arange,array,ones,linalg
from pylab import plot,show

# simple example
xi = arange(0,9)
A = array([ xi, ones(9)])
print(A)
y = [19, 20, 20.5, 21.5, 22, 23, 23, 25.5, 24]
# solve it
w = linalg.lstsq(A.T,y)[0]
print(w)
# -

# ###  2.  Statistical Relation for Simple Regression
#
# the theory of regression to the mean has some nice relationships. For example,
# to solve the following:
#
# $\mathbf{y} =  \mathrm{w_0} \mathbf{x} + \mathrm{w_1}$,
#
# use the statistical relations:
#
# $\mathrm{w_0} = \frac{\hat{\rho}(\mathbf{x},\mathbf{y}) \hat{\sigma}_y} {\hat{\sigma}_x}$
#
# and
#
# $\mathrm{w_1} = \bar{\mathbf{y}} - \mathrm{w_0} \bar{\mathbf{x}}$
#
#
# where $\hat{\rho}$ is the sample correlation between **x** and **y** and $\hat{\sigma}$ is the sample standard deviation.
#

# test the above 
w0_stat = pstats.correlation(xi,y)*pstats.standard_deviation(y)/pstats.standard_deviation(xi)
w1_stat = pstats.mean(y) - w0_stat*pstats.mean(xi)
print(w0_stat,w1_stat)

# +
# could have used scipy
from scipy import stats
slope, intercept, r_value, p_value, std_err = stats.linregress(xi,y)
print(slope, intercept, r_value, p_value, std_err)

line = w[0]*xi+w[1]
plot(xi,line,'r-',xi,y,'o')
show()


# -

# ### 3. Analytical Solution for Simple Regression
#
# The analytical solution for the optimal fit is done by first creating a cost function, here we use a residual sum of squares cost function,
#
# $J_{RSS} = \sum_i^n \left( y_i - \hat{y}_i \right) ^2$
#
# which pits the outcomes ($y_i$) against our predictions ($\hat{y}_i$).
#
# $J_{RSS} = \sum_i^n \left( y_i - \left( \mathrm{w_0} \mathbf{x} + \mathrm{w_1} \right) \right) ^2$
#
# Now to obtain the optimal fit parameters we take the gradient of the cost function
#
#
# $\frac{\partial{J_{RSS}}}{\partial{\mathrm{w_0}}} = -2 \sum_i^n \left( y_i - \left( \mathrm{w_0} \mathbf{x} + \mathrm{w_1} \right) \right)\mathbf{x}  $
#
# and 
#
# $\frac{\partial{J_{RSS}}}{\partial{\mathrm{w_1}}} = -2 \sum_i^n \left( y_i - \left( \mathrm{w_0} \mathbf{x} + \mathrm{w_1} \right) \right) $
#
# and by setting the gradient to zero and solving, we obtain the slope $\mathrm{w_0}$,
#
# $\mathrm{w_0} =\frac{ y_i x_i -  y_i  x_i / n }{  x_i^2 -  x_i  x_i / n}$
#
#
# and intercept ($\mathrm{w_1}$):
#
# $\mathrm{w_1} = \sum_i y_i/n - \mathrm{w_0}\sum_i x_i/n$
#
#

# we will need a dot product
def dot(x,y):
    return sum(x_i * y_i for x_i, y_i in zip(x,y))


def solve_simple_regression(x,y):
    """ returns w_0 (slope) and w_1 (intercept)"""
    n = float(len(x))
    slope = (dot(y,x) - ((sum(y)*sum(x))/n))/(dot(x,x)-((sum(x)*sum(x))/n))
    intercept = sum(y)/n - slope*sum(x)/n
    return slope, intercept


slope, intercept = solve_simple_regression(xi,y)
print(slope, intercept)

# ## Bayesian analysis

# ## Summary
#
# * Some simple statistics **do** differ for the anscombe dataset (median etc)
#
# * outliers can complicate interpretation and should be carefully considered

#  note see file mkm_notebooks/license.txt for license of this notebook.

StarcoderdataJupyterScriptsDedupFiltered