library("ggplot2")  # Пакет для красивых графиков 
library("grid")     # Пакет для субплотов
library('rio')  

# Отрегулируем размер картинок, которые будут выдаваться в нашей тетрадке
library('repr')
options(repr.plot.width=4, repr.plot.height=3)

lnL <- function(l){
    return(100*log(l) - 30*l)
  }

2*(lnL(10/3) - lnL(2))  # наблюдаемое значение
qchisq(0.95, df=1)      # критическое значение

df = import('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/flat.csv') 
head(df)

x = log(df$price)

lnL <- function(param,data) {
  n <- length(data)
  result <- -n/2*log(2*pi)-n/2*log(param[2]^2)-sum((data-param[1])^2)/(2*param[2]^2)
  return(result)
}

library('maxLik') 

res <- maxLik(lnL, start=c(0.1,0.1),data=x)
summary(res)

# вбили функцию правдоподобия с ограничением, нашли её оптимум
lnL_res <- function(sigma,data) lnL(c(5, sigma),x)

res_restr <- maxLik(lnL_res, start=0.1,data=x)
summary(res_restr)

2*(lnL(res$estimate, x) - lnL_res(res_restr$estimate, x))

qchisq(0.95, df=1)

mu0 = 4.8 
s0 = 0.3
2*(lnL(res$estimate, x) - lnL_res(c(mu0, s0), x))

qchisq(0.95, df=2)

library('dplyr')
library('rio') # тут import, который сам догадается до разделителей и даже excel прочтёт

df = import('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/verizon.csv')
# ILEC - внутренний клиент 
# CLEC - внешние клиенты 
head(df)

df %>% group_by(Group) %>% summarise(count = n(), mean = mean(Time), median = median(Time))

qplot(df$Time)

df %>% sample_n(size = 3, replace=TRUE)

get_bootstrap_samples <- function(data, B_sample){
    samples = list() # лист, куда мы будем записывать выборки
    
    for(i in 1:B_sample){
        # генерируем выборки с повторениями соответствующую по размерам оригинальной  
        # и записываем её в лист
        samples[[i]] = sample(data, size = length(data), replace=TRUE)
        }
    return(samples)
    }

# неколько сбустрапированных выборк
# get_bootstrap_samples(df$Time, 3)

stat_intervals <- function(stat, alpha){
    return(quantile(stat, c(alpha/2,(1-alpha/2)), name = FALSE))
    }

ilec = df %>% filter(Group == 'ILEC')
clec = df %>% filter(Group == 'CLEC')

# Делаем 1000 выборок для первой таблички и считаем по каждой медиану 
# функция lapply применяет функцию median к каждой выборке 
# команда as.numeric сделает листы векторами

clec_median_scores = as.numeric(lapply(get_bootstrap_samples(clec$Time, 10000), median))
ilec_median_scores = as.numeric(lapply(get_bootstrap_samples(ilec$Time, 10000), median))

cat("95% confidence interval for the ILEC median repair time:", stat_intervals(ilec_median_scores, 0.05),'\n')
cat("95% confidence interval for the CLEC median repair time:", stat_intervals(clec_median_scores, 0.05))

delta_median_scores = clec_median_scores - ilec_median_scores
cat("95% confidence interval for the diff median repair time:", stat_intervals(delta_median_scores, 0.05))

sum(delta_median_scores > 0)/length(delta_median_scores)





n_st = 72 
m_st = 18

n_lan = 49
m_lan = 11

n_gr = 41
m_gr = 12 

n_nw = 105
m_nw = 41 

alpha = 0.05

proportions_diff_zstat = function(n1, m1, n2, m2){ 
    p1 = m1/n1
    p2 = m2/n2
    P = (p1*n1 + p2*n2)/(n1 + n2)
    z_stat = (p1 - p2)/sqrt(P * (1-P) * (1/n1 + 1/n2))
    return(z_stat)
}

proportions_diff_ztest = function(z_stat, alternative = 'two-sided'){
    if(alternative == 'two-sided'){
        return(2 * (1 - pnorm(np.abs(z_stat))))
    }
    
    if(alternative == 'less'){
        return(pnorm(z_stat))
    }
    
    if(alternative == 'greater'){
        return(1 - pnorm(z_stat))
    }
}

# Для всех трёх случаев гипотеза о рпвенстве долей не отвергается. 
proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_gr, m_gr),'greater')
proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_lan, m_lan),'greater')
proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_nw, m_nw),'greater')

pvalues = rep(0,3)
pvalues[1] = proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_gr, m_gr),'greater')
pvalues[2] = proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_lan, m_lan),'greater')
pvalues[3] = proportions_diff_ztest(proportions_diff_zstat(n_st, m_st, n_nw, m_nw),'greater')
pvalues > 0.05/3 # не отвергается 

sort(pvalues)

seq(3,1,by=-1)
alpha/seq(3,1,by=-1)

sort(pvalues) > alpha/seq(3,1,by=-1) # гипотеза не отвергается 

p.adjust(c(0.01,0.02,0.07), "bonferroni")

p.adjust(c(0.01,0.02,0.07), "holm")

p.adjust(c(0.01,0.02,0.03), "hochberg")