#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Blatt 3

@author: ich
"""

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

M=1000
NN = np.array([ 100])
nr_nus = 201
nu = np.linspace(0,2,nr_nus)

gauss = np.random.randn(M)
cauchy = np.random.standard_cauchy(M)
expo = np.random.standard_exponential(M)

plt.figure('Plot_distr')
plt.hist(gauss,40,normed=True,range=[-5,5],alpha=0.5, label = "Gauss")
plt.hist(cauchy,40,normed=True,range=[-5,5],alpha=0.5, label = "Cauchy")
plt.hist(expo,40,normed=True,range=[-5,5],alpha=0.5, label = "Exp")
plt.legend(loc='best')
plt.ylabel('Probability density P(x)')
plt.xlabel('x')
plt.show()

print("Exercise 1 and 2, testing")
reject_norm = np.empty([len(NN),len(nu)])
reject_cauchy = np.empty([len(NN),len(nu)])
reject_expo = np.empty([len(NN),len(nu)])

reject_norm_rank = np.zeros([len(NN),len(nu)], dtype=float)
reject_cauchy_rank = np.zeros([len(NN),len(nu)], dtype=float)
reject_expo_rank = np.zeros([len(NN),len(nu)], dtype=float)

for i_n in range(len(NN)):
    x_tmp_norm = np.random.randn(M,NN[i_n])
    x_tmp_cauchy = np.random.standard_cauchy([M,NN[i_n]])
    x_tmp_expo = np.random.standard_exponential([M,NN[i_n]])
    for i_nu in range(len(nu)):
        y_tmp_norm = np.random.randn(M,NN[i_n])+nu[i_nu]
        t_norm = stats.ttest_ind(x_tmp_norm,y_tmp_norm)        
        reject_norm[i_n,i_nu] = float(sum(t_norm.pvalue<0.05))/NN[i_n]
        
        y_tmp_cauchy = np.random.standard_cauchy([M,NN[i_n]])+nu[i_nu]      
        t_cauchy = stats.ttest_ind(x_tmp_cauchy,y_tmp_cauchy)
        reject_cauchy[i_n,i_nu] = float(sum(t_cauchy.pvalue<0.05))/NN[i_n]
              
        y_tmp_expo = np.random.exponential(1+nu[i_nu],size=[M,NN[i_n]])
        t_expo = stats.ttest_ind(x_tmp_expo,y_tmp_expo)
        reject_expo[i_n,i_nu] = float(sum(t_expo.pvalue<0.05))/NN[i_n]
        
        for grr in range(NN[i_n]):
            rank_norm = stats.wilcoxon(x_tmp_norm[:,grr],y_tmp_norm[:,grr])
            reject_norm_rank[i_n,i_nu] += float(rank_norm.pvalue<0.05)/NN[i_n]
            
            rank_cauchy = stats.wilcoxon(x_tmp_cauchy[:,grr],y_tmp_cauchy[:,grr])
            reject_cauchy_rank[i_n,i_nu] += float(rank_cauchy.pvalue<0.05)/NN[i_n]
            
            rank_expo = stats.wilcoxon(x_tmp_expo[:,grr],y_tmp_expo[:,grr])
            reject_expo_rank[i_n,i_nu] += float(rank_expo.pvalue<0.05)/NN[i_n]
        
        
plt.figure('ttest fig')
plt.plot(np.linspace(0,2,nr_nus),reject_norm.T,'r', label = "Normal t-test")
plt.plot(np.linspace(0,2,nr_nus),reject_cauchy.T,'g', label = "Cauchy t-test")
plt.plot(np.linspace(0,2,nr_nus),reject_expo.T,'b', label = "Exp t-test")
plt.plot(np.linspace(0,2,nr_nus),reject_norm_rank.T,'r--', label = "Normal Wilcoxon")
plt.plot(np.linspace(0,2,nr_nus),reject_cauchy_rank.T,'g--', label = "Cauchy Wilcoxon")
plt.plot(np.linspace(0,2,nr_nus),reject_expo_rank.T,'b--', label = "Exp Wilcoxon")
plt.xlabel("v")
plt.ylabel("power")
plt.legend(loc='best')
plt.show()