#!/usr/bin/env python
# coding: utf-8

# # Khoảng tin cậy (Confidence interval)

# In[1]:


# import libraries
import math
import numpy as np
import scipy.stats


# ## Sai số chuẩn (Standard error)
# - Ước lượng sai số dựa vào tỉ lệ (for proportion)
# $$\sigma_{\hat{p}} = \sqrt{\frac{p(1-p)}{n}}$$

# **Bài toán: tỉ lệ khách hàng thoả mãn với dịch vụ hiện tại**

# In[2]:


# standard error for proportion
# p hat: estimate proportion for population
p = 0.6
n = 100
std_err_prop = math.sqrt(p * (1 - p) / n) * 100
print "Sai số ước lượng:", std_err_prop
print "Tỉ lệ chặn dưới:", p * 100 - std_err_prop
print "Tỉ lệ chặn trên:", p * 100 + std_err_prop


# In[3]:


# increase sample size will decrease standard error
n = 1000
std_err_prop = math.sqrt(p * (1 - p) / n) * 100
print "Sai số ước lượng:", std_err_prop
print "Tỉ lệ chặn dưới:", p * 100 - std_err_prop
print "Tỉ lệ chặn trên:", p * 100 + std_err_prop


# - Ước lượng sai số từ giá trị trung bình (For mean)
# $$\sigma_{\bar{x}} = \frac{\sigma}{\sqrt{n}}$$

# **Bài toán: thời gian trung bình để một khách hàng uống hết ly cà phê**

# In[4]:


# lấy mẫu 5 khách hàng vào mỗi ngày thứ Hai như sau
n = 5
sample_A = np.asarray([2.1, 1.3, 0.6, 1.5, 2.4])
sample_B = np.asarray([0.7, 1.3, 1.8, 0.7, 1.9])
sample_C = np.asarray([2.1, 1.9, 2.5, 1.8, 0.9])
sample_D = np.asarray([1.6, 2.1, 1, 0.9, 1.2])
print "Sample A mean:", sample_A.mean()
print "Sample B mean:", sample_B.mean()
print "Sample C mean:", sample_C.mean()
print "Sample D mean:", sample_D.mean()

average_of_means = (sample_A.mean() + sample_B.mean() + sample_C.mean() + sample_D.mean()) / 4
list_of_means = np.asarray([sample_A.mean(), sample_B.mean(), sample_C.mean(), sample_D.mean()])
std_of_means = list_of_means.std()
print "Thời gian trung bình từ các mẫu trung bình:", average_of_means
print "Độ lệch của từ các mẫu trung bình:", std_of_means


# In[5]:


# estimate standard deviation of population
std_of_population = std_of_means * math.sqrt(n)
print "Ước lượng độ lệch chuẩn cho quần thể:", std_of_population
print "Chặn dưới cho ước lượng:", average_of_means - std_of_means
print "Chặn trên cho ước lượng:", average_of_means + std_of_means


# ## Cách tính khoảng tin cậy
# $$proportionCI = \hat{p} \pm Zscore_{\alpha/2} * \frac{\sigma}{\sqrt{n}}$$
# $$meanCI = \bar{x} \pm Zscore_{\alpha/2} * \frac{\sigma}{\sqrt{n}}$$
# $$uperLimit = sampleParam + Zscore_{\alpha/2} * standardError$$
# $$lowerLimit = sampleParam - Zscore_{\alpha/2} * standardError$$
# 
# **Bài toán: dự đoán khả năng trúng cử của ứng viên A so với ứng viên B với độ tin cậy là 95%**

# In[6]:


# lấy mẫu 100 người
n = 100

# ủng hộ ứng viên A là 55%
candidate_A_prop = 0.55

# ủng hộ ứng viên B là 45%
candidate_B_prop = 0.45

# tính Zscore cho ứng viên A với độ tin cậy là 95%
qnorm = scipy.stats.norm.ppf(0.975)
print "Z-score:", qnorm

# calculate p hat
p = candidate_A_prop
std_err_prop = math.sqrt(p * (1 - p) / n)
print "Độ lỗi chuẩn cho tỉ lệ ủng hộ ứng viên A:", std_err_prop

# calculate standard error of mean
std_of_means = std_err_prop / math.sqrt(n)
print "Độ lỗi chuẩn cho tỉ lệ thu thập được:", std_of_means

# calculate interval
upper_limit = candidate_A_prop + qnorm * std_of_means
lower_limit = candidate_A_prop - qnorm * std_of_means
print "Tỉ lệ ủng hộ chặn trên:", upper_limit
print "Tỉ lệ ủng hộ chặn dưới:", lower_limit


# **Bài toán: ước lượng thời gian trung bình (phút) một người chạy trong vòng 1km với độ tin cậy là 95%**

# In[7]:


# thống kê thời gian các lượt chạy
list_run_time = np.asarray([11.3, 10.3, 9.6, 9.1, 8.9, 8.9, 8.1, 7.5, 6.9])
sample_size = list_run_time.size
sample_mean = list_run_time.mean()
sample_std = list_run_time.std()
standard_err = sample_std / math.sqrt(sample_size)

print "Total runners:", sample_size
print "Sample mean:", sample_mean
print "Sample standard deviation:", sample_std
print "Standard/Sampling error:", standard_err

# calculate interval
qnorm = scipy.stats.norm.ppf(0.975)
print "Z-score:", qnorm
upper_limit = sample_mean + qnorm * standard_err
lower_limit = sample_mean - qnorm * standard_err
print "Thời gian trung bình chặn trên:", upper_limit
print "Thời gian trung bình chặn dưới:", lower_limit


# # Kiểm định giả thuyết (Hypothesis testing)
# Thực hiện 4 bước
# 1. Thiệt lập giả thuyết $H_0, H_a$ và chọn mức ý nghĩa/ngưỡng chấp nhận/độ không chắc chắn/độ tin cậy $\alpha$
# 2. Xác định dạng kiểm định (binomial-test, z-test, t-test, chi-square-test, f-test)
# 3. Tính p-value
# 4. Đối chiếu p-value với mức ý nghĩa đã chọn

# **Bài toán: nữ sinh bị gọi lên bảng trả bài nhiều hơn nam đúng hay sai?**

# In[8]:


# gọi tên 50 người thì thấy tỉ lệ nam nữ như sau
n_boys = 14
n_girls = 36


# 1. Thiết lập giả thuyết
#     - $H_0$ cho rằng nữ bị gọi là do tình cờ (happend by chance) nghĩa là xác suất bị gọi $p \le 0.5$
#     - $H_a$ cho rằng nữ bị gọi là do thiên vị (not by chance) nghĩa là xác suất bị gọi $p > 0.5$
#     - Chọn mức ý nghĩa $\alpha = 5\%$. Nếu tỉ lệ nữ bị gọi nhiều hơn 36 có xác suất nhỏ hơn 5% thì ta sẽ loại bỏ $H_0$.
# 2. Xác định dạng kiểm định:
#     - Đây là binomial-test với $p = 0.5$
# 3. Tính p-value (chính là Zscore)        

# In[9]:


n = n_boys + n_girls
x = n_girls
p = 0.51
alpha = 0.05
p_value = scipy.stats.binom.pmf(x, n, p)
print "p-value/Z-score:", p_value


# In[10]:


# reject H0 or not
if p_value < alpha:
    print "Reject H0 => nữ bị gọi là do thiên vị"
else:
    print "Accept H0 => nữ bị gọi là do tình cờ"


# **Bài toán: test khả năng thắng cử của ứng cử viên**
# 1. Thiết lập giả thuyết
#     - $H_0: p \le 0.5$ - ứng viên thua cử
#     - $H_a: p > 0.5$ - ứng viên thắng cử
#     - Chọn mức ý nghĩa $\alpha = 5\%$.
# 2. Xác định dạng kiểm định:
#     - Đây là z-test với $z_p = \frac{\hat{p} - p_0}{\sqrt{\frac{p_0(1 - p_0)}{n}}}$
# 3. Tính p-value (chính là Zscore)        

# In[11]:


# lấy mẫu thực nghiệm
n = 500
p_hat = yes = 0.54
p_0 = 0.5
z_p = (p_hat - p_0)/(math.sqrt(p_0 * (1 - p_0) / n))
print "Z-score:", z_p

# one-tail test
z_cdf = scipy.stats.norm.cdf(z_p)
print "Z-cdf:", z_cdf

# calculate p-value
p_value = 1 - z_cdf
print "p-value:", p_value


# In[12]:


# reject H0 or not
if p_value < alpha:
    print "Reject H0 => ứng viên thắng cử"
else:
    print "Accept H0 => ứng viên thua cử"


# **Bài toán: khách hàng không thích sản phẩm nhỏ hơn 20kg, trước khi ship nhà sản xuất cần phải test điều này**
# 1. Thiết lập giả thuyết
#     - $H_0: \mu \ge 20.15 \ kg$
#     - $H_a: \mu < 20.15 \ kg$
#     - Chọn mức ý nghĩa $\alpha = 5\%$.
# 2. Xác định dạng kiểm định:
#     - Đây là z-test với $z = \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}}$
# 3. Tính p-value (chính là Zscore)        

# In[13]:


# lấy mẫu thực nghiệm
n = 100
mu = 20.15
avg_weight = 20.1
std_weight = 0.26
z = (avg_weight - mu)/(std_weight / math.sqrt(n))
print "Z-score:", z

# one-tail test
p_value = scipy.stats.norm.cdf(z)
print "p-value:", p_value


# In[14]:


# reject H0 or not
if p_value < alpha:
    print "Reject H0 => không ship hàng"
else:
    print "Accept H0 => ship hàng"