# -*- coding: utf-8 -*-
"""
Created on Thu Nov 11 14:53:54 2021

@author: Sim
"""
def test_var_2samp(var1, n1, var2, n2, alternative='two-sided'):
    import scipy.stats as st
    
    df1 = n1 - 1; df2 = n2 - 1 
    F0 = var1/var2
   
    if alternative == 'greater':
        pval = 1 - st.f.cdf(F0, df1, df2)
    elif alternative == 'less':
        pval = st.f.cdf(F0,df2, df2)
    else:
        pval = st.f.cdf(F0,df1,df2)
        if pval < 0.5:
            pval = pval * 2
        else:
            pval = (1-pval)*2
    
    return(F0, pval)

def ci_var_2samp(var1, n1, var2, n2, level=95):
    import scipy.stats as st
    
    F0 = var1/var2
    df1 = n1 - 1; df2 = n2 - 1 
    
    lbd = F0/st.f.ppf(1-(100-level)/100/2, df1, df2)
    ubd = F0/st.f.ppf((100-level)/100/2, df1, df2)
    return((lbd,ubd))

F1 = test_var_2samp(15,10, 10, 10)
print(F1)

ci1 = ci_var_2samp(15, 10, 10, 10, level=95)
print(ci1)

#----------------
x = [0.005, 0.600, 1.325, 1.505, 1.520, 2.000, 2.255, 2.575, 
     3.000, 3.300, 3.500] 
y = [0.005, 0.005, 0.255, 0.505, 0.525, 1.000, 1.250, 1.505, 
     1.750, 2.750]
import numpy as np
Fval = test_var_2samp(np.var(x, ddof=1), len(x), np.var(y, ddof=1), len(y))
print(Fval)
