import numpy as np
from scipy.stats import norm
 
mu = 4.55
sigma = 0.06
X_bar = 4.57
n = 10
alpha = 0.05
Z_alpha_div_2 = norm.ppf(1 - alpha / 2)
Z = (X_bar - mu) / (sigma / np.sqrt(n))
Z, Z_alpha_div_2

import matplotlib.pyplot as plt
from scipy.stats import norm
import numpy as np
 
alpha = 0.05
plt.plot(np.linspace(-3, 3, 100), norm.pdf(np.linspace(-3, 3, 100)))  # 绘制曲线
plt.fill_between(np.linspace(-3, norm.ppf(alpha / 2), 50),
                 norm.pdf(np.linspace(-3, norm.ppf(alpha / 2), 50)), 
                 alpha=.15, color="red", label="Rejection region |Z|>k")
plt.fill_between(np.linspace(norm.ppf(1 - alpha / 2), 3, 50),
                 norm.pdf(np.linspace(norm.ppf(1 - alpha / 2), 3, 50)), 
                 alpha=.15, color="red")
plt.fill_between(np.linspace(norm.ppf(alpha / 2), norm.ppf(1 - alpha / 2), 50),
                 norm.pdf(np.linspace(norm.ppf(alpha / 2), norm.ppf(1 - alpha / 2), 50)), 
                 alpha=.15, color="blue", label="Acceptance region |Z|<= k")
plt.legend()
plt.show()  # 显示图片

import pandas as pd
import pylab
import math
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
import scipy.stats
import warnings
warnings.filterwarnings('ignore')
 
df = pd.read_csv('normtemp.txt', sep='   ', names=['Temperature', 'Gender', 'Heart Rate'])
df.describe()

df.head()

observed_temperatures = df['Temperature'].sort_values()
bin_val = np.arange(start= observed_temperatures.min(), stop= observed_temperatures.max(), step = .05)
mu, std = np.mean(observed_temperatures), np.std(observed_temperatures)
p = norm.pdf(observed_temperatures, mu, std)
plt.hist(observed_temperatures,bins = bin_val, density=True, stacked=True)
plt.plot(observed_temperatures, p, color = 'red')
plt.xticks(np.arange(95.75,101.25,0.25),rotation=90)
plt.xlabel('Human Body Temperature Distributions')
plt.xlabel('human body temperature')
plt.show()
print('Average (Mu): '+ str(mu) + ' / ' 'Standard Deviation: '+str(std))

# Shapiro-Wilk Test: https://en.wikipedia.org/wiki/Shapiro%E2%80%93Wilk_test
"""
语法:scipy.stats.shapiro(x, a=None, reta=False)
x 为待检验的数据
返回(统计数, P 值)
"""
shapiro_test, shapiro_p = scipy.stats.shapiro(observed_temperatures)
print("Shapiro-Wilk Stat:",shapiro_test, "\nShapiro-Wilk p-Value:", shapiro_p)
"""
语法 scipy.stats.normaltest(a, axis=0, nan_policy='propagate')
a 是待检验数据
axis 正态分布测试将沿其计算的轴
返回：(对数据集进行假设检验的 k2 值, P 值)
"""
k2, p = scipy.stats.normaltest(observed_temperatures)
print('p:',p)

# stats.probplot(QQ 图)
scipy.stats.probplot(observed_temperatures, dist="norm", plot=pylab)
pylab.show()

def ecdf(data):
    """Compute ECDF"""
    n = len(data)
    x = np.sort(data)
    y = np.arange(1, n+1) / n
    return x, y
 
# Compute empirical mean and standard deviation 计算经验平均值和标准差
 
# Number of samples
n = len(df['Temperature']) 
 
# Sample mean 样本均值
mu = np.mean(df['Temperature']) 
 
# Sample standard deviation 样本方差
std = np.std(df['Temperature']) 
 
print('Mean temperature: ', mu, 'with standard deviation of +/-', std)
 
# Random sampling of the data based off of the mean of the data.
# 根据数据的平均值对数据进行随机抽样。
normalized_sample = np.random.normal(mu, std, size=10000)
x_temperature, y_temperature = ecdf(df['Temperature'])
normalized_x, normalized_y = ecdf(normalized_sample)
 
# Plot the ECDFs
fig = plt.figure(figsize=(8, 5))
plt.plot(normalized_x, normalized_y)
plt.plot(x_temperature, y_temperature, marker='.', linestyle='none')
plt.ylabel('ECDF')
plt.xlabel('Temperature')
plt.legend(('Normal Distribution', 'Sample data'))

from scipy import stats
 
CW_mu = 98.6
# 计算单独样本 t 检验
stats.ttest_1samp(df['Temperature'], CW_mu, axis=0)

female_temp = df.Temperature[df.Gender == 2]
male_temp = df.Temperature[df.Gender == 1]
 
# 当不确定两总体方差是否相等时, 应先利用 levene 检验, 检验两总体是否具有方差齐性
rvs1 = stats.norm.rvs(female_temp)
rvs2 = stats.norm.rvs(male_temp)
print(stats.levene(rvs1, rvs2))
 
mean_female_temp = np.mean(female_temp)
mean_male_temp = np.mean(male_temp)
print('Average female body temperature = ' + str(mean_female_temp))
print('Average male body temperature = ' + str(mean_male_temp))
 
# Compute independent t-test 
"""
计算两个独立分数样本的平均值的 T-test。
这是对 2 个独立样本具有相同平均(预期)值的零假设的检验。此测试假定默认情况下总体具有相同的方差。
"""
print(stats.ttest_ind(female_temp, male_temp, axis=0))

import pandas as pd
import numpy as np
from scipy import stats
 
data = pd.io.stata.read_stata('us_job_market_discrimination.dta')
data.head()

blacks = data[data.race == 'b']
whites = data[data.race == 'w']
blacks.call.describe()

whites.call.describe()

blacks_called = len(blacks[blacks['call'] == True])
blacks_not_called = len(blacks[blacks['call'] == False])
whites_called = len(whites[whites['call'] == True])
whites_not_called = len(whites[whites['call'] == False])
observed = pd.DataFrame({'blacks': {'called': blacks_called, 'not_called': blacks_not_called},
                         'whites': {'called' : whites_called, 'not_called' : whites_not_called}})
observed

# 统计获得职位和未获得职位的人数
num_called_back = blacks_called + whites_called
num_not_called = blacks_not_called + whites_not_called
# 得到获得职位的理论比率
rate_of_callbacks = num_called_back / (num_not_called + num_called_back)
# 获得职位和未获得职位的理论值
expected_called = len(data) * rate_of_callbacks
expected_not_called = len(data) * (1 - rate_of_callbacks)
# 理论值四格表(黑人和白人的录取率相等)
expected = pd.DataFrame({'blacks':{'called': expected_called / 2, 'not_called': expected_not_called / 2},
                        'whites': {'called': expected_called / 2, 'not_called': expected_not_called / 2}})
expected

import scipy.stats as stats
 
observed_frequencies = [blacks_not_called, whites_not_called, whites_called, blacks_called]
expected_frequencies = [expected_not_called / 2, expected_not_called / 2, expected_called / 2, expected_called / 2]
# 卡方检验
stats.chisquare(f_obs=observed_frequencies, f_exp=expected_frequencies)

import numpy as np
from scipy import stats
 
mean1 = 30.97
mean2 = 21.79
 
std1 = 26.7
std2 = 12.1
 
nobs1 = 10
nobs2 = 10
 
modified_std1 = np.sqrt(np.float32(nobs1) / np.float32(nobs1-1)) * std1
modified_std2 = np.sqrt(np.float32(nobs2) / np.float32(nobs2-1)) * std2
 
(statistic, pvalue) = stats.ttest_ind_from_stats(mean1=mean1, std1=modified_std1, nobs1=10, mean2=mean2, std2=modified_std2, nobs2=10)
 
print("t statistic is: ", statistic)
print("pvalue is: ", pvalue)

from  scipy.stats import chi2_contingency
import numpy as np
 
kf_data = np.array([[37,27], [39,21]])
kf = chi2_contingency(kf_data)
print('chisq-statistic=%.4f, p-value=%.4f, df=%i expected_frep=%s' % kf)