感謝Dr.fish的耐心講解和細(xì)致回答殴俱。
本次課的隨堂作業(yè)如下：
同一類(lèi)動(dòng)物的肱骨大概具有相同的長(zhǎng)寬比，考古學(xué)家根據(jù)這一性質(zhì)來(lái)鑒定物種枚抵∠哂考古學(xué)家發(fā)掘了41塊肱骨化石，假設(shè)它們來(lái)自于同一物種俄精，判斷它們是不是物種A（已知物種A的肱骨長(zhǎng)寬比為8.5）询筏。取α=0.01。
進(jìn)階作業(yè)為：
使用計(jì)算機(jī)模擬方法計(jì)算課堂作業(yè)
直接上代碼
# 導(dǎo)入分析包
import numpy as np
import pandas as pd
import scipy.stats
# 加載數(shù)據(jù)
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]
課堂作業(yè)解題思路
- 僅有一組數(shù)據(jù)竖慧，故排除配對(duì)及雙樣本檢驗(yàn)嫌套；
- 非整體數(shù)據(jù)，故排除z檢驗(yàn)圾旨；
- 最終選擇t檢驗(yàn)踱讨。
t statisic
手工計(jì)算過(guò)程：
- 設(shè)置原假設(shè)和備擇假設(shè)
原假設(shè) 及 備擇假設(shè)
檢驗(yàn)分布圖
- 計(jì)算t統(tǒng)計(jì)量
1. 隨堂作業(yè)部分
方法一：手工計(jì)算單樣本 t 檢驗(yàn)
# 計(jì)算 t 統(tǒng)計(jì)量
mean = np.mean(data)
n = len(data)
t_statistics = (mean - 8.5) / (np.std(data, ddof = 1) / np.sqrt(n))
t_statistics
# 輸出結(jié)果
4.0303238468687361
# 計(jì)算 t 臨界值
# alpha = 0.01， 雙邊檢驗(yàn)砍的，取α/2
t_critical = scipy.stats.t.isf(0.01 / 2 , df = n - 1)
t_critical
# 輸出結(jié)果
2.7044592674331502
#計(jì)算 p 值驗(yàn)證下
p_value = scipy.stats.t.sf(t_statistics, df = n - 1) * 2 # 雙邊檢驗(yàn)p值×2
p_value
# 輸出結(jié)果
0.00024267304119873163
方法二：使用 scipy.stats.ttest_1samp()計(jì)算單樣本 t 檢驗(yàn)計(jì)算
t_statistics, p_value = scipy.stats.ttest_1samp(data, 8.5)
t_statistics, p_value
# 輸出結(jié)果
(4.030323846868737, 0.00024267304119873163)
結(jié)論
當(dāng)置信度水平α取0.01時(shí)痹筛，因?yàn)?t_statistics > t_critical (或者 p_value < α) , t統(tǒng)計(jì)量落在拒絕域中，所以拒絕原假設(shè)，樣本不是物種A帚稠。
2. 進(jìn)階作業(yè) - bootstrap 方法
bootstrap解題思路
bootstrap為有放回抽樣谣旁，所以需要：
- 定義一個(gè)單次抽樣計(jì)算P值函數(shù)
- 定義一個(gè)重復(fù)抽樣并計(jì)算p值均值的函數(shù)
- 輸入抽樣次數(shù)
# 第一種方法 -- 公式手算
# 定義一個(gè)單次抽樣并計(jì)算P值的函數(shù)，并返回單次抽樣的P值
def bs_replicate_1d(data):
bs_replicate_1d = np.random.choice(data, size=len(data))
stats_mean = np.mean(bs_replicate_1d)
stats_n = len(bs_replicate_1d)
t_statistics = (stats_mean - 8.5) / (np.std(bs_replicate_1d, ddof = 1) / np.sqrt(stats_n)) #計(jì)算t統(tǒng)計(jì)量
stats_p = scipy.stats.t.sf(t_statistics, df = stats_n - 1) * 2
return stats_p
# 定義一個(gè)重復(fù)抽樣并計(jì)算p值均值的函數(shù)滋早，返回p值的均值
def draw_bs_reps(data, size=1):
bs_replicates = np.empty(size) # 初始一個(gè)空數(shù)組
# 進(jìn)行多次重新抽樣
for i in range(size):
bs_replicates[i] = bs_replicate_1d(data)
p_value = sum(bs_replicates) / size
return p_value
# 重復(fù)抽樣10000次
draw_bs_reps(data, size=1000)
# 輸出結(jié)果
0.0057400175854891624
# 第二種方法 -- 套函數(shù)
def bs_replicate_1d(data):
bs_replicate_1d = np.random.choice(data, size=len(data))
t_statistics, stats_p = scipy.stats.ttest_1samp(bs_replicate_1d, 8.5)
stats_p = stats_p
return stats_p
# 定義一個(gè)重復(fù)抽樣并計(jì)算p值均值的函數(shù)榄审，返回p值的均值
def draw_bs_reps(data, size=1):
bs_replicates = np.empty(size) # 初始一個(gè)空數(shù)組
# 進(jìn)行多次重新抽樣
for i in range(size):
bs_replicates[i] = bs_replicate_1d(data)
p_value = sum(bs_replicates) / size
return p_value
# 重復(fù)抽樣10000次
draw_bs_reps(data, size=1000)
# 輸出結(jié)果
0.0045788518825149857
不曉得為什么用計(jì)算機(jī)模擬算出來(lái)的p值會(huì)比用樣本數(shù)據(jù)算出來(lái)的大（雖然也是落在拒絕域內(nèi)吧……）
