import numpy as np
import scipy.stats
一個(gè)樣本的均值
例1：某工廠的苯含量檢測
給出某鋼鐵廠20分空氣樣本的苯含量檢測數(shù)據(jù)（單位ppm），判斷該工廠空氣是否達(dá)標(biāo)（苯含量小于百萬分之1，即 1ppm）。 （數(shù)據(jù)來自 Mendenhall所著的《統(tǒng)計(jì)學(xué)》）
data = [0.21, 1.44, 2.54, 2.97, 0.00, 3.91, 2.24, 2.41, 4.50, 0.15,
0.30, 0.36, 4.50, 5.03, 0.00, 2.89, 4.71, 0.85, 2.60, 1.26]
方法一：手工計(jì)算
mean = np.mean(data)
n = len(data)
mean, n
(2.1435000000000004, 20)
t_statistics = (mean - 1) / (np.std(data, ddof=1) / np.sqrt(n))
t_statistics
2.9457560457212408
# alpha = 0.05， 單邊檢驗(yàn)
t_critical = scipy.stats.t.isf(0.05, df=n-1)
t_critical
1.7291328115213678
p_value = scipy.stats.t.sf(t_statistics, df=n-1)
p_value#為什么復(fù)制過來下面的值不一樣了
0.0041496038528359052
當(dāng)置信度水平 αα 取0.05時(shí)，因?yàn)?t_statistics > t_critical (或者 p_value < αα ) , t統(tǒng)計(jì)量落在拒絕域中，所以拒絕原假設(shè)，即該工廠的空氣中苯含量超過標(biāo)準(zhǔn)值。
方法二：使用 scipy.stats.ttest_1samp()
scipy.stats.ttest_1samp(data, 1) # 給出的是雙邊檢驗(yàn)的結(jié)果
Ttest_1sampResult(statistic=2.9457560457212408, pvalue=0.0082992077056718103)
兩個(gè)配對的樣本的均值比較
例2：水的金屬含量
飲用水中的金屬會(huì)影響水的口感，如果濃度過高甚至?xí)】诞a(chǎn)生危害。這里有10份飲用水，分別測量它們底部的水與表面水中的含鋅的濃度。判斷底部的水中鋅的濃度是否和表面水的濃度一樣？（數(shù)據(jù)來源 https://onlinecourses.science.psu.edu/stat500/node/51）
bottom = [0.430, 0.266, 0.567, 0.531, 0.707, 0.716, 0.651, 0.589, 0.469, 0.723]
surface = [0.415, 0.238, 0.390, 0.410, 0.605, 0.609, 0.632, 0.523, 0.411, 0.612]
scipy.stats.ttest_rel(bottom, surface)
Ttest_relResult(statistic=4.8638127451351831, pvalue=0.00089111545782254793)
scipy.stats.t.isf(0.05/2, df=9)
2.262157162740992
當(dāng)置信度水平 αα 取0.05時(shí)， 因?yàn)? p?value<αp?value<α ， 所以拒絕原假設(shè)，即兩種水的鋅濃度不一樣。
兩個(gè)獨(dú)立的樣本的均值比較
例3 ：自動(dòng)打包機(jī)器
某打包工廠，用機(jī)器來包裝紙箱。給出新、老機(jī)器的打包時(shí)間數(shù)據(jù)（單位：秒），判斷新機(jī)器是否比舊機(jī)器打包得更快？ （數(shù)據(jù)來源 https://onlinecourses.science.psu.edu/stat500/node/50 ）
old = [42.7, 43.8, 42.5, 43.1, 44.0, 43.6, 43.3, 43.5, 41.7, 44.1]
new = [42.1, 41.3, 42.4, 43.2, 41.8, 41.0, 41.8, 42.8, 42.3, 42.7]
scipy.stats.ttest_ind(new, old)
Ttest_indResult(statistic=-3.3972307061176026, pvalue=0.0032111425007745158)
t_statistics, p_value = scipy.stats.ttest_ind(old, new)
p_value = p_value / 2
p_value
0.0016055712503872579
-scipy.stats.t.isf(0.05, df=len(old)+len(new)-2) # 注意自由度df
-1.7340636066175359
當(dāng)置信度水平 αα 取0.05時(shí)， 因?yàn)? p?value<αp?value<α ， 所以拒絕原假設(shè)，即新機(jī)器打包更快。
計(jì)算機(jī)模擬之 bootstrap 方法
data = [0.21, 1.44, 2.54, 2.97, 0.00, 3.91, 2.24, 2.41, 4.50, 0.15,
0.30, 0.36, 4.50, 5.03, 0.00, 2.89, 4.71, 0.85, 2.60, 1.26]
def bootstrap_replicate_1d(data, func): # 進(jìn)行一次重新抽樣，并返回檢驗(yàn)統(tǒng)計(jì)量
return func(np.random.choice(data, size=len(data)))
def draw_bs_reps(data, func, size=1):
bs_replicates = np.empty(size) # 初始一個(gè)空數(shù)組
for i in range(size): # 進(jìn)行多次重新抽樣
bs_replicates[i] = bootstrap_replicate_1d(data, func)
return bs_replicates # 返回多次抽樣的檢驗(yàn)統(tǒng)計(jì)量數(shù)組
def bootstrap_pvalue_1samp(data, pop_stats, func, size=1):
sample_stats = func(data) # 計(jì)算原有樣本的檢驗(yàn)統(tǒng)計(jì)量
translated_data = data - sample_stats + pop_stats # 數(shù)據(jù)平移
bs_replicates = draw_bs_reps(translated_data, func, size) # 重新抽樣
p = np.sum( bs_replicates > sample_stats) / size # 計(jì)算抽樣統(tǒng)計(jì)量大于原有統(tǒng)計(jì)量的概率，根據(jù)實(shí)際情況來
return p
bootstrap_pvalue_1samp(data, 1, np.mean, size=10000)
0.0016999999999999999
P值小于0.05， 拒絕原假設(shè)
計(jì)算機(jī)模擬之 Permutation 方法
old = [42.7, 43.8, 42.5, 43.1, 44.0, 43.6, 43.3, 43.5, 41.7, 44.1]
new = [42.1, 41.3, 42.4, 43.2, 41.8, 41.0, 41.8, 42.8, 42.3, 42.7]
def diff_of_means(data_1, data_2):
diff = np.mean(data_1) - np.mean(data_2) # 計(jì)算兩組數(shù)據(jù)均值的差
return diff
def permutation_sample(data1, data2): # 產(chǎn)生新的分組數(shù)據(jù)
data = np.concatenate((data1, data2)) # 合并兩組數(shù)據(jù)
permuted_data = np.random.permutation(data) # 對合并后的數(shù)據(jù)進(jìn)行重新排列
perm_sample_1 = permuted_data[:len(data1)] # 分成新的兩組數(shù)據(jù)
perm_sample_2 = permuted_data[len(data1):]
return perm_sample_1, perm_sample_2
def draw_perm_reps(data_1, data_2, func, size=1): # 進(jìn)行多次重新分組的操作
perm_replicates = np.empty(size)
for i in range(size):
perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)
perm_replicates[i] = func(perm_sample_1, perm_sample_2)
return perm_replicates
def permutation_pvalue(data_1, data_2, func, size=1): # 計(jì)算P值
empirical_test_stats = func(data_1, data_2)
perm_replicates = draw_perm_reps(data_1, data_2, func, size)
p = np.sum(perm_replicates < empirical_test_stats) / len(perm_replicates) # 根據(jù)實(shí)際情況修改
return p
permutation_pvalue(new, old, diff_of_means, size=10000)#值都有點(diǎn)不一樣
0.0012999999999999999
P值小于0.05， 拒絕原假設(shè)
基本作業(yè)：
同一類動(dòng)物的肱骨大概具有相同的長寬比，考古學(xué)家根據(jù)這一性質(zhì)來鑒定物種。考古學(xué)家發(fā)掘了41塊肱骨化石，假設(shè)它們來自于同一物種，判斷它們是不是物種A（已知物種A的肱骨長寬比為8.5）。取 α=0.01α=0.01 。
數(shù)據(jù)為： [10.73, 8.89, 9.07, 9.20, 10.33, 9.98, 9.84, 9.59, 8.48, 8.71, 9.57, 9.29, 9.94, 8.07, 8.37, 6.85, 8.52, 8.87, 6.23, 9.41, 6.66, 9.35, 8.86, 9.93, 8.91, 11.77, 10.48, 10.39, 9.39, 9.17, 9.89, 8.17, 8.93, 8.80, 10.02, 8.38, 11.67, 8.30, 9.17, 12.00, 9.38]
data = [10.73, 8.89, 9.07, 9.20, 10.33, 9.98, 9.84, 9.59, 8.48, 8.71, 9.57, 9.29, 9.94, 8.07, 8.37, 6.85, 8.52, 8.87, 6.23, 9.41, 6.66, 9.35, 8.86, 9.93, 8.91, 11.77, 10.48, 10.39, 9.39, 9.17, 9.89, 8.17, 8.93, 8.80, 10.02, 8.38, 11.67, 8.30, 9.17, 12.00, 9.38]
#怎么設(shè)置全部顯示哦
#表示完全看不懂了.不知道怎么做這個(gè)，還有上面的圖片怎么弄進(jìn)去的也不知道
#markdown的格式也沒對。
