<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 1. Assumptions of most tests
1.1 Normality
Step1:  How to test it
- QQ plot. Measure our quantiles vs hypothetical quantiles.
- Most tests are resistant to small deviations.
Step2:  1.2 Equal variance (if two or more groups)
Step3:  How to test it
- Levene test (it doesn't assume normality)
Step4:  1.3 i.i.d. (independent and identically distributed)
This comes from the data collection.
- Each observation is independent
- There are no subpopulations in the population
1.4 Uncorrelated residuals (for linear regression)
Step5:  1.4 No or little multicolinearity (for linear regression)
Your dependent variables are not measuring the same thing
2. t-test
2.1 Sample against value
Imagine you measure increase in productivity, and you want to see if the increase is significant
Step6:  2.2 Two independent samples
Imagine you measure productivity for a control group and a experimental group
Imagine it's normally distributed 
You want to see if the productivity of the groups is different
Step7:  2.3 Two paired samples
Imagine you measure productivity for a group before starting a new program and after the program (each individual has two observations)
Imagine it's normally distributed
You can test if the productivity has increased
Step8:  2.4 What happens if they don't have equal variance? (But normality)
Step9:  2.5 What happens if they are not normally distributed (but equal variance)
Use only when the number of observation in each sample is > 20 and you have 2 independent samples of ranks
Mann-Whitney-U
Step10:  2.6 What happens if they are not normally distributed (but equal variance) and paired samples?
Use only when the number of observation in each sample is > 20 and you have 2 independent samples of ranks
Wilcoxon signed ranksum
Step11:  2.6 What happens if they are not normally distributed and not equal variance?
You can only test if the distributions are different (booo)
scipy.stats.mannwhitneyu(x,y)
3. More than two samples
ANOVA (Parametric) or Kruskal-Wallis (Non parametric)
Analysis of variance (ANOVA) is a collection of statistical models used to analyze the differences between group means and their associated procedures (such as "variation" among and between groups). 
- ~ t-test
- If normality
The Kruskal–Wallis one-way analysis of variance by ranks (named after William Kruskal and W. Allen Wallis) is a non-parametric method for testing whether samples originate from the same distribution. 
- ~ Mann-Whitney-U 
- If no-normality but equal variance
3.1 Compare if the groups are different
Step12:  3.2 Problem with many groups (look-elsewhere bias)
Step13:  3.3 Compare which groups are different
3.3.1 Bonferroni correction
Corrects for the number of comparisons that you WILL make (Set p-value before the experiment. Don't cheat please)
New threhold to consider something significant is $\frac{0.05}{#\ of\ comparisons}$
Step14:  3.3.2 Tukey's HSD (honest significant difference) test. Best
Assumes normality
Much more power (less false negatives) than Bonferroni (bonferroni with many samples is hard)
Step15:  4. Fitting distributions
Why do we care about distributions?
They tell us a lot about the mechanism generating our values
Normal distribution
Step16:  5. Calculate Confidence Intervals with bootstrap methods
How does it work?
You have a sample (imagine 100 values)
Do 10000 times
Step17:  In class exercise
Step18:  PartA. Compare the number of employees of industrial vs financial companies
Step19:  PartB. Compare the market capitalization for different types of companies
<ASSISTANT_TASK:>
Python Code:
Image(url="https://kanbanize.com/blog/wp-content/uploads/2014/07/Standard_deviation_diagram.png", width=500)
Explanation: 1. Assumptions of most tests
1.1 Normality
End of explanation
def qq_plot(x):
    import scipy.stats
    (osm, osr),(slope, intercept, r)  = scipy.stats.probplot(x, dist='norm', plot=None)
    plt.plot(osm, osr, '.', osm, slope*osm + intercept)
    plt.xlabel('Quantiles',fontsize=14)
    plt.ylabel('Quantiles Obs',fontsize=14)  
x = np.random.randn(100)
plt.subplot(2,2,1)
qq_plot(x)
x = np.random.exponential(2,100)
plt.subplot(2,2,2)
qq_plot(x)
plt.tight_layout()
plt.show()
Explanation: How to test it
- QQ plot. Measure our quantiles vs hypothetical quantiles.
- Most tests are resistant to small deviations.
End of explanation
Image(url="http://goodsciencebadscience.nl/wp-content/uploads/2012/09/variances1.gif",width=500)
Explanation: 1.2 Equal variance (if two or more groups)
End of explanation
#random normal data (mean 0 std 1)
x = np.random.randn(100)
#random normal data (mean 0 std 2)
y = np.random.randn(100)*2
#Test
print(scipy.stats.levene(x,y))
Explanation: How to test it
- Levene test (it doesn't assume normality)
End of explanation
print("These are the residuals")
Image(url="https://upload.wikimedia.org/wikipedia/commons/e/ed/Residuals_for_Linear_Regression_Fit.png",width=500)
Image(url="https://i.stack.imgur.com/RU17l.png")
Explanation: 1.3 i.i.d. (independent and identically distributed)
This comes from the data collection.
- Each observation is independent
- There are no subpopulations in the population
1.4 Uncorrelated residuals (for linear regression)
End of explanation
from scipy.stats import norm
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.randn(100)+2
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=False,fit=norm)
plt.subplot(2,2,2)
qq_plot(x)
plt.show()
print("-"*20)
print("Test for mean of group different than zero")
scipy.stats.ttest_1samp(x,popmean=0)
Explanation: 1.4 No or little multicolinearity (for linear regression)
Your dependent variables are not measuring the same thing
2. t-test
2.1 Sample against value
Imagine you measure increase in productivity, and you want to see if the increase is significant
End of explanation
from scipy.stats import norm
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.randn(100)+0.1
#Our other sample, normally distributed with mean 0.5 and std 1
y = np.random.randn(100)+0.5
#Equal variance?
#We are testing if we can do the test
#If significant it means that the variances are different
print("Test for equal variance")
print(scipy.stats.levene(x,y))
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=True)
sns.distplot(y,kde=True)
plt.subplot(2,2,2)
qq_plot(x)
qq_plot(y)
plt.tight_layout()
plt.show()
print("-"*20)
print("Test for different means")
scipy.stats.ttest_ind(x,y)
Explanation: 2.2 Two independent samples
Imagine you measure productivity for a control group and a experimental group
Imagine it's normally distributed 
You want to see if the productivity of the groups is different
End of explanation
from scipy.stats import norm
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.randn(100)+0.1
#Our other sample, similar to x but adding 0.05 and some random noise
y = x+np.random.randn(100)/10 + 0.05
#Equal variance?
print("Test for equal variance")
print(scipy.stats.levene(x,y))
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=True)
sns.distplot(y,kde=True)
plt.subplot(2,2,2)
qq_plot(x)
qq_plot(y)
plt.tight_layout()
plt.show()
print("WRONG: What would happen if we don't use paired t-test")
print(scipy.stats.ttest_ind(x,y))
print("-"*20)
print("Test for different means")
scipy.stats.ttest_rel(x,y)
Explanation: 2.3 Two paired samples
Imagine you measure productivity for a group before starting a new program and after the program (each individual has two observations)
Imagine it's normally distributed
You can test if the productivity has increased
End of explanation
scipy.stats.ttest_rel?
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.randn(100)+0.1
#Our other sample, normally distributed with mean 0.2 and std 2
y = np.random.randn(50)*2+0.2
#Equal variance?
print("Test for equal variance")
print(scipy.stats.levene(x,y))
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=True)
sns.distplot(y,kde=True)
plt.subplot(2,2,2)
qq_plot(x)
qq_plot(y)
plt.tight_layout()
plt.show()
print("WRONG: What happens if we try the test assumming equal variance")
print(scipy.stats.ttest_ind(x,y))
print("-"*20)
print("Test for different means with different variance")
scipy.stats.ttest_ind(x,y,equal_var=False)
Explanation: 2.4 What happens if they don't have equal variance? (But normality)
End of explanation
#Our sample, uniformly distributed 
x = np.random.exponential(2,100)
#Our other sample, niformly distributed  + 0.1
y = np.random.exponential(2,100)+0.1
#Equal variance?
print("Test for equal variance")
print(scipy.stats.levene(x,y))
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=True)
sns.distplot(y,kde=True)
plt.subplot(2,2,2)
qq_plot(x)
qq_plot(y)
plt.show()
print("WRONG: What happens if we try t-test")
print(scipy.stats.ttest_ind(x,y))
print("-"*20)
print("Test for different means with no normal distributions, equal variance")
scipy.stats.mannwhitneyu(x,y)
Explanation: 2.5 What happens if they are not normally distributed (but equal variance)
Use only when the number of observation in each sample is > 20 and you have 2 independent samples of ranks
Mann-Whitney-U
End of explanation
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.random(100)
#Our other sample, normally distributed with mean 0.05 and std 1 and some random noise
y = x+np.random.randn(100)/10 + 0.05
#Equal variance?
print("Test for equal variance")
print(scipy.stats.levene(x,y))
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=True)
sns.distplot(y,kde=True)
plt.subplot(2,2,2)
qq_plot(x)
qq_plot(y)
plt.show()
print("WRONG: What happens if we try t-test")
print(scipy.stats.ttest_rel(x,y))
print("-"*20)
print("Test for different means with no normal distributions and equal variance")
print(scipy.stats.wilcoxon(x,y))
Explanation: 2.6 What happens if they are not normally distributed (but equal variance) and paired samples?
Use only when the number of observation in each sample is > 20 and you have 2 independent samples of ranks
Wilcoxon signed ranksum
End of explanation
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.randn(100)+0.1
#Our other sample, normally distributed with mean 0.2 and std 1
y = np.random.randn(100)+0.2
#Our other sample, normally distributed with mean 0.3 and std 1
z = np.random.randn(100)+0.3
#Equal variance?
print("Test for equal variance")
print(scipy.stats.levene(x,y,z))
#histogram
plt.subplot(2,2,1)
sns.distplot(x,kde=True)
sns.distplot(y,kde=True)
sns.distplot(z,kde=True)
plt.subplot(2,2,2)
qq_plot(x)
qq_plot(y)
qq_plot(z)
plt.show()
print("-"*20)
print("ONE-WAY ANOVA")
print(scipy.stats.f_oneway(x, y, z)    )
print("-"*20)
print("ONE-WAY KRUSKAL WALLIS") #use if not-normal and equal variance
print(scipy.stats.mstats.kruskalwallis(x,y,z))
Explanation: 2.6 What happens if they are not normally distributed and not equal variance?
You can only test if the distributions are different (booo)
scipy.stats.mannwhitneyu(x,y)
3. More than two samples
ANOVA (Parametric) or Kruskal-Wallis (Non parametric)
Analysis of variance (ANOVA) is a collection of statistical models used to analyze the differences between group means and their associated procedures (such as "variation" among and between groups). 
- ~ t-test
- If normality
The Kruskal–Wallis one-way analysis of variance by ranks (named after William Kruskal and W. Allen Wallis) is a non-parametric method for testing whether samples originate from the same distribution. 
- ~ Mann-Whitney-U 
- If no-normality but equal variance
3.1 Compare if the groups are different
End of explanation
#Let's compare two identical samples, see how often we get false results
for i in range(100):
    x = np.random.randn(100)
    y = np.random.randn(100)
    stat, pvalue = scipy.stats.ttest_ind(x,y)
    if pvalue < 0.05: 
        print(i,pvalue)
    
Explanation: 3.2 Problem with many groups (look-elsewhere bias)
End of explanation
#Our sample, normally distributed with mean 0.1 and std 1
x = np.random.randn(100)+0.1
#Our other sample, normally distributed with mean 0.2 and std 1
y = np.random.randn(100)+0.2
#Our other sample, normally distributed with mean 0.3 and std 1
z = np.random.randn(100)+0.3
#Our other sample, normally distributed with mean 0.3 and std 1
w = np.random.randn(100)+0.5
print("New threshold (instead of 0.05) = ",0.05/6)
#Compare x and y
print("x and y: ", scipy.stats.ttest_ind(x,y))
#Compare x and z
print("x and z: ", scipy.stats.ttest_ind(x,z))
#Compare x and w
print("x and w: ", scipy.stats.ttest_ind(x,w))
#Compare y and z
print("y and z: ", scipy.stats.ttest_ind(y,z))
#Compare y and w
print("y and w: ", scipy.stats.ttest_ind(y,w))
#Compare z and w
print("z and w: ", scipy.stats.ttest_ind(z,w))
Explanation: 3.3 Compare which groups are different
3.3.1 Bonferroni correction
Corrects for the number of comparisons that you WILL make (Set p-value before the experiment. Don't cheat please)
New threhold to consider something significant is $\frac{0.05}{#\ of\ comparisons}$
End of explanation
from statsmodels.stats.multicomp import pairwise_tukeyhsd
df = pd.read_csv("data/tukey_example.csv")
df.sample(10)
res2 = pairwise_tukeyhsd(df["productivity"],df["group"])
print(res2)
res2.plot_simultaneous(comparison_name=None,xlabel='diffs',ylabel='Group')
plt.show()
Explanation: 3.3.2 Tukey's HSD (honest significant difference) test. Best
Assumes normality
Much more power (less false negatives) than Bonferroni (bonferroni with many samples is hard)
End of explanation
from scipy.stats import lognorm,norm,expon,powerlaw
x = np.random.lognormal(1,0.8,3000) 
plt.subplot(2,2,1)
sns.distplot(x,fit=norm,kde=False)
plt.title("norm")
plt.subplot(2,2,2)
sns.distplot(x,fit=lognorm,kde=False)
plt.title("lognorm")
plt.subplot(2,2,3)
sns.distplot(x,fit=expon,kde=False)
plt.title("expon")
plt.subplot(2,2,4)
sns.distplot(x,fit=powerlaw,kde=False)
plt.title("powerlaw")
plt.tight_layout()
plt.show()
Explanation: 4. Fitting distributions
Why do we care about distributions?
They tell us a lot about the mechanism generating our values
Normal distribution: 
Nothing weird going on (just normal)
Lognormal distribution: 
Multiplicative process: 
The value depends on the multiplication of many normal variables.
For example, each unit of time you can grow depending on some random factors.
Many financial indicators are lognormally distributed
Looks like a normal distribution if you take the logarithm of your values (np.log(x))
Exponential distribution: 
Time between events (generated by poisson events)
If the time between events follows an exponential distribution the events are independent
It's a straight line if the y axis is in log scale
Power-law distribution: 
Rich-get-richer mechanism: Your growth depends on your value -> Inequality is unavoidable
The tails of some human and financial phenomena follow power-laws.
The tail is a straight line in log-log scale (both y and x axis)
Binomial distribution: Number of successes out of a number of trials. Each event independent
Poisson distribution: Number of successes in some time. Each event independent
End of explanation
import scikits.bootstrap as bootstrap
#Bootstrap of mean and std
x = np.random.randn(100)*3 #this could be one of our columns
CIs = bootstrap.ci(x, statfunction=np.mean,n_samples=100000)  
print('CI for mean with bootstrapping = ', CIs)
CIs = bootstrap.ci(x, statfunction=np.std,n_samples=100000)  
print('CI for std with bootstrapping = ', CIs)
x = np.random.randn(100)*3
np.mean(x),np.std(x)
CIs = bootstrap.ci(x, statfunction=np.mean,n_samples=100000)  
print('CI for mean with bootstrapping = ', CIs)
scipy.stats.ttest_1samp(x,popmean=0)
#Bootsrap of p-value
def return_pvalue(x):
    stat,pvalue = scipy.stats.ttest_1samp(x,0)
    return pvalue
CIs = bootstrap.ci(x, statfunction=return_pvalue, n_samples=10000)  
print('CI for p-value with bootstrapping = ', CIs)
Explanation: 5. Calculate Confidence Intervals with bootstrap methods
How does it work?
You have a sample (imagine 100 values)
Do 10000 times: 
Take 100 values from the sample (with replacement)
Calculate the mean (or std)
Take the percentiles 2.5 and 97.5 of those 10000 values -> Those are your confidence intervals (kind-of)
End of explanation
df = pd.read_csv("data/big3_position.csv",sep="\t")
df.head()
Explanation: In class exercise
End of explanation
industrial = df.loc[df["TypeEnt"]=="Industrial company"]
financial = df.loc[df["TypeEnt"]=="Financial company"]
#Employees (like most is lognormally distributed -> we can convert it to a normal distribution and run our tests happily)
i = np.log(industrial["Employees"].dropna())
f = np.log(financial["Employees"].dropna())
from scipy.stats import lognorm,norm
plt.figure(figsize=(10,3))
plt.subplot(1,2,1)
sns.distplot(df["Employees"].dropna(),fit=lognorm,kde=False)
plt.subplot(1,2,2)
sns.distplot(np.log(df["Employees"].dropna()),fit=norm,kde=False)
i.head()
f.head()
#Check if our distributions have the same variance
scipy.stats.levene(i,f)
sns.distplot(i)
sns.distplot(f)
#Plot histograms and qq_plots for i and f
qq_plot(i)
qq_plot(f)
#Run the correct statistical test to know if the 
#number of employees of financial and industrial companies is different
scipy.stats.man(i,f)
Explanation: PartA. Compare the number of employees of industrial vs financial companies
End of explanation
df = pd.read_csv("data/big3_position.csv",sep="\t")
df.head()
#Create a new column that is the logarithm of "MarketCap", call it "logMarketCap"
df["logMarketCap"] = np.log(df["MarketCap"])
df.head()
df.describe()
#Replace infinite values (log(0) = infinite) with nan
df = df.replace([np.inf, -np.inf], np.nan)
#Drop na values if "logMarketCap" is na 
df = df.dropna(subset=["logMarketCap"])
df.describe()
#Find the distribution of "MarketCap" and "logMarketCap", fitting the rigth distribution
from scipy.stats import norm,lognorm,expon
sns.distplot(df["MarketCap"],fit=lognorm,kde=False)
#Find the distribution of "MarketCap" and "logMarketCap", fitting the rigth distribution
from scipy.stats import norm,lognorm,expon
sns.distplot(df["logMarketCap"],fit=norm,kde=False)
#Kepe only three types of entities
df = df.loc[df["TypeEnt"].isin(["Financial company","Industrial company","Bank"])]
#Plot distributions for each type
plt.figure(figsize=(6,4))
for typeent,group in df.groupby("TypeEnt"):
    sns.distplot(group["logMarketCap"],kde=False,norm_hist=True,label=typeent)
    plt.legend()
#Run ANOVA
bank = df.loc[df["TypeEnt"] == "Bank"]
bank_values = bank["logMarketCap"]
ind = df.loc[df["TypeEnt"] == "Industrial company"]
ind_values = ind["logMarketCap"]
fin = df.loc[df["TypeEnt"] == "Financial company"]
fin_values = fin["logMarketCap"]
scipy.stats.f_oneway(bank_values,ind_values,fin_values)
#Run Tukey test
res2 = pairwise_tukeyhsd(df["logMarketCap"],df["TypeEnt"])
print(res2)
#Plot tukey test
res2.plot_simultaneous(comparison_name=None,xlabel='diffs',ylabel='Group')
plt.show()
Explanation: PartB. Compare the market capitalization for different types of companies
End of explanation
<END_TASK>