"""
var = sd^2 = summ(Xi^2 - mu)/N = summ(Xi^2)/N - mu^2

Observation is:
Biased_Variance = (n-1)/n Population/Total_Variance
OR,
summ(Xi^2 - mu)/n = (n-1)/n Population/Total_Variance
=> Population_Variance = summ(Xi^2 - mu)/(n-1)
is a better Variance estimate than Sample Variance
"""

from math import sqrt, pi, e
from itertools import izip, imap
from numpy import var
import scipy.stats
from pandas import DataFrame, read_pickle
#zscore, t score, f score
from scipy.stats import norm, t, f, manwhitneyu, shapiro

def mean(x): return reduce(lambda a,b: (a+b), x)/float(len(x))

"""
access like this df[0.01][38][4] - df[critical percent][df][num of samples]
#Make sure the q_table csv is in same parent directory else calc_tukey_hsd won't work
"""
q_table = read_pickle('q_table')

def mean_arr(x):
    sumr = 0.
    prod_list = []
    for i in x:
        sumr += i
        prod_list.append(i*i)
    return sumr/len(x), prod_list


def variance_biased(x):
    """x is a list of numbers"""
    mean_x, prod_list = mean_arr(x)
    return sum(prod_list)/float(len(x)) - (mean_x)**2

def variance_unbiased(x):
    """x is a list of numbers"""
    xbar = mean(x)
    return sum([(i-xbar)**2 for i in x])/float(len(x)-1)


def std_deviation_biased(x):
    return sqrt(variance_biased(x))

def std_deviation_unbiased(x):
    return sqrt(variance_unbiased(x))


def normal_distro_probability(x, xbar, var):
    z2 = float(x-xbar)**2/(2*var)  #z squared
    return 1./sqrt((2*pi*var) * (e**z2))


def normal_sample_variance(original_var, sample_size):
    """
    given original variance but not all data points.
    sqrt(var/n) = Std. Error
    """
    return original_var/sample_size


############# ALL ABOUT CONFIDENCE INTERVAL & INFERENTIAL STATS ################

#standard error
se = lambda sd, sample_size: (sd/sqrt(sample_size))

def se_independent_sample(sd1, sd2, n1, n2):
    """
    1. WELCH's T TEST
    2. USE WHEN WE CMP 2 DIFF MEANS
    ---------------------------------
    Useful for finding out
    CI and comparing Z stats Given a population metric
    """
    return sqrt((sd1**2 + sd2**2)/float(n1)) if n1==n2 else sqrt((sd1**2/float(n1)) + (sd2**2/float(n2)))


def df_two_sampled_welch(var1, var2, n1, n2):
    #IMP dof FUNC!!
    #2 sample Deg of Freedom for 2 Sample Welch's T Test
    df1, df2 = n1-1, n2-1
    return float(var1/float(n1) + var2/float(n2))**2/\
    float((var1**2/float(n1**2 * df1)) + (var2**2/float(n2**2 * df2)))


pooled_variance = lambda x,y,n1,n2: ((var(x)*n1) + (var(y)*n2))/float(n1+n2-2)

def se_pooled_t(x, y):
    """
    Useful when calculating T stats reports
    since sample metrics are given.
    x, y: arrays
    n1, n2: sizes of x, y respectively.
    """
    n1,n2 = len(x), len(y)
    #sd^2 = \summ(x-xbar)^2 + \summ(y-ybar)^2/df1+df2
    pooled_var = pooled_variance(x,y,n1,n2)
    return round(sqrt(pooled_var * (1./n1 + 1./n2)), 2)


df_independent_sample = lambda n1,n2: (n1+n2-2)


#marginal error given z score range
marginal_z = lambda z_score, se: z_score * se

def calc_z(mu, xbar, sd, sample_size):
    """
    Calculate Z score given normal sample distribution mean
    and population mean, its std dev and the sample size for the
    sample distribution.
    """
    #std_error = se(sd, sample_size)
    return float(xbar - mu)/se(sd, sample_size)

#half if two tailed else not and /100
def critical_z(percentile, one_tailed):
    #calculate critical z from critical percentage
    return round(abs(scipy.stats.norm.ppf((100 - percentile)/100.)), 3) if one_tailed \
    else round(abs(scipy.stats.norm.ppf((100 - percentile)/200.)), 3)

def calc_probability_sample_mean(mu, xbar, sd, sample_size):
    #calculate percentage from z score
    return scipy.stats.norm.cdf(calc_z(mu, xbar, sd, sample_size))

def calc_t(mu, xbar, sd, sample_size):
    #calculating critical t
    return round(calc_z(mu, xbar, sd, sample_size), 2)

def calc_t_independent_sample(mu, xbar, se):
    #calculating critical t when se is given/independent samples
    return round(float(xbar-mu)/se, 2)

def t_percentile(t_val, df, one_tailed=0):
    """
    Find P-Value Given T score, DF and if its 1 Tailed or 2 Tailed
    """
    return t.sf(t_val, df) if one_tailed else t.sf(t_val, df) * 2

def critical_t(percentile, df, one_tailed):
    return round(t.isf((100-percentile)/100., df), 3) if one_tailed else round(t.isf((100-percentile)/200., df), 3)

#same as critical_t, but in proportion
def t_val_from_t_percentile(t_percentile, df, one_tailed = 0):
    """
    Find T score given T percentile, DF
    and if its 1 Tailed or 2 Tailed
    t_percentile: % in proportion
    """
    return round(t.isf(t_percentile, df), 3) if one_tailed else round(t.isf(t_percentile/2., df), 3)


def t_cmp(calculated_t, critical_t):
    """
    Given Critical T Value and T score,
    Return if its Significant or Not.
    """
    return abs(calculated_t) > abs(critical_t)


### Correlation Measures using T Stats ####
def t_r_squared(t_score, df):
    #calculated t score cousin of r squared for line fitting - correlation
    """
    R square generally tells us how much of the difference between 2 means
    in a particular hypothetisis is due to the given statistical factors
    OR R squared Error;
    Higher the Value better the Model of Assumption.
    Lower the value implies we are missing out on some stat factors and that
    its a weak hypothesis.
    """
    t= t_score**2
    return t/float(t+df)


def t_significant(x, y, percentile, one_sided):
    """
    by calculating from independent sampling from a bunch of samples
    when population params are not given
    """
    df = df_independent_sample(len(x), len(y))
    t_critical = critical_t(percentile, df, one_sided)
    xbar= mean(x)
    ybar = mean(y)
    std_dev = se_pooled_t(x, y)
    t_score = calc_t_independent_sample(ybar, xbar, std_dev)
    return t_score, t_cmp(t_score, t_critical)



def ci_t_margin_error(sd, sample_size, df, t_score = None, percentile = None, one_tailed = 0):
    critical_t_score = 0.
    if percentile:
        critical_t_score = critical_t(percentile, df, one_tailed)
    elif t_score:
        critical_t_score = t_score  #this is the Critical T value!
    #std_error = se(sd, sample_size)
    return marginal_z(critical_t_score, se(sd, sample_size))



def ci_t(xbar, sd, sample_size, df, t_score = None, percentile = None, one_tailed = 0):
    margin_error = ci_t_margin_error(sd, sample_size, df, t_score = t_score, percentile = percentile, one_tailed = one_tailed)
    return (xbar - margin_error, xbar + margin_error)

def ci_margin_error(sd, sample_size, z_score = None, percentile = None, one_tailed = 0):
    """
    z score * std_error is the margine of error i C.I. from the sample mean
    """
    critical_z_score = 0.
    if percentile: critical_z_score = critical_z(percentile, one_tailed)
    elif z_score: critical_z_score = z_score
    #std_error = se(sd, sample_size)
    return marginal_z(critical_z_score, se(sd, sample_size))


#confidence interval
def ci(xbar, sd, sample_size, z_score = None, percentile = None, one_tailed = 0):
    """
    xbar: sample mean
    percentile: in %
    """
    margin_error = ci_margin_error(sd, sample_size, z_score = z_score, percentile = percentile, one_tailed=one_tailed)
    return (xbar - margin_error, xbar + margin_error)

def ci_generic(xbar, margin_error):
    return (xbar-margin_error, xbar+margin_error)

#Decision process component
def z_cmp(calculated_z_score_proportion, criticalz_percentage_proportion):
    """
    calculated_z_score_proportion: the z score calculated from mu and xbar in proportion
    criticalz_percentage_proportion: Given Critical Value proportion for acceptance criteria
    if calculated_z_score_proportion > criticalz_percentage_proportion
    then the chance of that happening is p < criticalz_percentage_proportion
    """
    return abs(calculated_z_score_proportion) > scipy.stats.norm.cdf(scipy.stats.norm.ppf((100-criticalz_percentage_proportion)/100.))


######## ANALYSIS OF VARIANCE - ANOVA ###########
"""
F = Between Group Variability/Within Group Variability
"""
f_critical = lambda alpha, dfb, dfw: round(f.isf(alpha, dfb, dfw), 4)
#between group var
#xbar_groups is an array of x means
df_ssb = lambda xbar_groups: float(len(xbar_groups) - 1)

ssb = lambda xg, xbar_groups, each_grp_sample_size: \
sum([ni * (xbar - xg)**2 for ni, xbar in izip(each_grp_sample_size, xbar_groups)])/df_ssb(xbar_groups)

#within group var
#x_samples is an array of lists - this is like df for pooled variance
df_ssw = lambda x_samples: float(reduce(lambda x,y: x+y, imap(len, x_samples)) - len(x_samples))

#this is like pooled variance for more than 2 samples - generic
ssw = lambda x_samples: sum([var(samples) * len(samples) for samples in x_samples])/df_ssw(x_samples)

#F ratio stats
anova_ratio = lambda xg, xbar_groups, x_samples, each_grp_sample_size:ssb(xg, xbar_groups, each_grp_sample_size)/ssw(x_samples)

#if greater we reject null hypothesis else we fail to reject
f_cmp = lambda f_score, f_critical: abs(f_score) > abs(f_critical)

###### CALCULATING EFFECT MEASURES (TUKEYS HSD)- measures of significance not due to samplign error or chances #########

#this can be used when the num of samples per group are not the same N
var_num_samples = lambda x_samples: len(x_samples)/reduce(lambda n1,n2: (1/float(n1) + 1/float(n2)), imap(len, x_samples))

"""
#honesty significant differences - which samples are significantly different?
#q = df[critical_per][dfw][num_samples]
This is how you calculate Q:
Q = ml - ms / sqrt of mssw / n
k = num of groups
n = Total num of elements per group across all groups for an experiment where each group has same elements, else refer to var_num_samples(x_samples)
to calculate n;
"""

calc_tukey_hsd = lambda critical_per, dfw, ssw_score, k, n: round(q_table[critical_per][dfw][k] * sqrt(float(ssw_score)/n), 4)

def tukey_hsd_cmp(mean_diff, tukey_hsd):
    """
    Reject Null if mean diff > tukey hsd
    """
    return mean_diff > tukey_hsd

#cohens d multiple comparisons ANOVA
cohens_d_multiple = lambda x1,x2,ssw_score: float(x1-x2)/sqrt(ssw_score)

#eta squared - proportion of explained differences - sort of r sqrd
#Analogous to R squared Difference of T Test;
eta_sqrd = lambda ssb_score, ss: float(ssb_score)/ss
#ss = total or ssb+ssw = \summ((xi - xg)**2) summation of all xi minus grand xg


######## OTHER TESTS ON NON GAUSSIAN DATASETS ###############
"""
Man - Whitney U-Test
u,p = manwhitneyu(x_arr, y_arr)
Non-parametric test;
- non-normal data
- find descriptive variables.
"""
manwhitneyu_test = lambda x_arr, y_arr: manwhitneyu(x_arr, y_arr)

"""
Shapiro-Wilk Test
- Plot Histogram to check for Bell Curve
- Perform Shapiro-Wilk Test
s, p = shapiro(x_arr)
"""
shapiro_test = lambda x_arr: shapiro(x_arr)
