# -*- coding: utf-8 -*-
"""
Created on Wed Nov 10 16:59:41 2021

@author: Sim
"""
import statsmodels.stats.proportion as prop
import scipy.stats as st
import numpy as np

n1=100; x1=70; n2=1000; x2=540
p0 = (x1+x2)/(n1+n2); p1hat=x1/n1; p2hat = x2/n2

z0, pval = prop.test_proportions_2indep(x1, n1, x2, n2, value=0.1, 
                  compare='diff', alternative='two-sided')
print('검정값/유의확률: ', z0, pval)

cidp1 = prop.confint_proportions_2indep(x1, n1, x2, n2, method=None, 
                  compare='diff', alpha=0.05)
print(cidp1)

z1, pval1 = prop.test_proportions_2indep(x1, n1, x2, n2, value=0., 
                  compare='diff', alternative='two-sided')
print('검정값/유의확률: ', z1, pval1)

#### 직접계산
sedp0 = np.sqrt( p0*(1-p0)*(1/n1+1/n2) )
sedp1 = np.sqrt( p1hat*(1-p1hat)/n1 + p2hat*(1-p2hat)/n2 )

z2 = (p1hat - p2hat -0.1)/sedp1
pval2 = (1-st.norm.cdf(np.abs(z2)))*2
print(z2, pval2)

cidp2 = st.norm.interval(.95, loc=p1hat-p2hat, scale=sedp1)
print(cidp2)

z3 = (p1hat - p2hat)/sedp0
pval3 = (1-st.norm.cdf(np.abs(z3)))*2
print(z3, pval3)