# -*- coding: utf-8 -*-
"""
Created on Wed Sep  1 16:30:23 2021

@author: Rene
"""
import numpy as np
from scipy.stats import norm
from scipy.stats import t
import matplotlib.pyplot as plt

def get_conf_int(x,alpha):
    n = len(x)
    x_bar = np.mean(x)
    st_dev_x_bar = np.std(x,ddof=1)/(n**0.5)
    half_width = t.ppf(1-(alpha/2),n-1)*st_dev_x_bar
    bounds = np.zeros(2)
    bounds[0] = x_bar-half_width
    bounds[1] = x_bar+half_width
    return(bounds)

n = 10
the_mu = 0
the_sigma = 1
p = 0.90
m = 50
conf_bounds = np.zeros( (m, 2) )
mu_in_conf = np.zeros(m)

i_range = np.arange(0,m)

for i in i_range:
    sample = np.random.normal(the_mu, the_sigma, size = n)
    x_bar = np.mean(sample)
    conf_bounds[i,] = get_conf_int(sample,1-p)
    mu_in_conf[i] = (conf_bounds[i,0] <= the_mu) and (conf_bounds[i,1] >= the_mu)

# Plotting the analysis results    
LB_x = 0
UB_x = m+1
off_set_x = 1

LB_y = -1.5
UB_y = 1.5
off_set_y = 0.05
    
# Setting the attributes for the Figure     
plt.rc('font', family='serif', size='10')
plt.figure(figsize=(10, 5))
Figure = plt.figure()
Figure.subplots_adjust(hspace=0.4, wspace=0.4)
x_lims = (LB_x-off_set_x,UB_x+off_set_x)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('i')
plt.ylabel('confindence interval for $\mu$')

for i in i_range:
    if mu_in_conf[i]:
        conf_color = 'green'
    else:
        conf_color = 'red'
    half_width = (conf_bounds[i,1] - conf_bounds[i,0])/2
    mid_point = (conf_bounds[i,1] + conf_bounds[i,0])/2
    plt.errorbar(i, mid_point, yerr = half_width, ecolor=conf_color, 
             elinewidth=1, capsize=0,fmt='-o',ms=3,mfc=conf_color,mec=conf_color, label = 'conf.int', c = conf_color);

not_captured = m - np.sum(mu_in_conf)

text_str = 'not captured = ' + str(not_captured)
plt.text(1,-1.5,text_str,color ='red')

text_str = str(m)+ " sampled " + f'{(100*p):2.0f}' + " % confidence intervals - variance unknown"
Figure.suptitle(text_str,size='14')
plt.savefig('Confidence_Interval_Sample_Var_Known.png', dpi=1200)