六侥钳、SciPy 統(tǒng)計推斷

原文：statistical-inference-scipy

譯者：飛龍

協(xié)議：CC BY-NC-SA 4.0

6.1 效應(yīng)量

署名：派生于 Allen Downey 的 CompStats呈枉。協(xié)議：Creative Commons Attribution 4.0 International。

from __future__ import print_function, division

import numpy
import scipy.stats

import matplotlib.pyplot as pyplot

from IPython.html.widgets import interact, fixed
from IPython.html import widgets

# 為隨機(jī)數(shù)生成器播種堪侯，所以我們得到相同結(jié)果
numpy.random.seed(17)

# 來自 http://colorbrewer2.org/ 的一些漂亮的顏色
COLOR1 = '#7fc97f'
COLOR2 = '#beaed4'
COLOR3 = '#fdc086'
COLOR4 = '#ffff99'
COLOR5 = '#386cb0'

%matplotlib inline

為了探索量化效應(yīng)量的統(tǒng)計量嚎尤，我們將研究男女之間的身高差異。 我使用來自行為風(fēng)險因素監(jiān)測系統(tǒng)（BRFSS）的數(shù)據(jù)伍宦，來估計美國成年女性和男性的身高的平均值和標(biāo)準(zhǔn)差（cm）芽死。

我將使用scipy.stats.norm來表示分布乏梁。結(jié)果是一個rv對象（代表隨機(jī)變量）。

mu1, sig1 = 178, 7.7
male_height = scipy.stats.norm(mu1, sig1)

mu2, sig2 = 163, 7.3
female_height = scipy.stats.norm(mu2, sig2)

以下函數(shù)在平均值的 4 個標(biāo)準(zhǔn)差內(nèi)評估正態(tài)（高斯）概率密度函數(shù)（PDF）关贵。它需要rv對象并返回一對 NumPy 數(shù)組遇骑。

def eval_pdf(rv, num=4):
    mean, std = rv.mean(), rv.std()
    xs = numpy.linspace(mean - num*std, mean + num*std, 100)
    ys = rv.pdf(xs)
    return xs, ys

這是兩個分布的樣子。

xs, ys = eval_pdf(male_height)
pyplot.plot(xs, ys, label='male', linewidth=4, color=COLOR2)

xs, ys = eval_pdf(female_height)
pyplot.plot(xs, ys, label='female', linewidth=4, color=COLOR3)
pyplot.xlabel('height (cm)')
None

png

我們現(xiàn)在假設(shè)這些是總體的真實(shí)分布揖曾。當(dāng)然落萎，在現(xiàn)實(shí)生活中，我們從未觀察到真實(shí)的總體分布炭剪。我們通常需要使用隨機(jī)樣本练链。

我將使用rvs從總體分布中生成隨機(jī)樣本。請注意奴拦，這些是完全隨機(jī)的媒鼓，完全有代表性的樣品，沒有測量誤差错妖！

male_sample = male_height.rvs(1000)

female_sample = female_height.rvs(1000)

兩個樣本都是 NumPy 數(shù)組绿鸣。 現(xiàn)在我們可以計算樣本統(tǒng)計量，如均值和標(biāo)準(zhǔn)差站玄。

mean1, std1 = male_sample.mean(), male_sample.std()
mean1, std1

# (178.16511665818112, 7.8419961712899502)

樣本均值接近總體均值枚驻，但不如預(yù)期的那樣精確。

mean2, std2 = female_sample.mean(), female_sample.std()
mean2, std2

# (163.48610226651135, 7.382384919896662)

女性樣本的結(jié)果相似株旷。

Now, there are many ways to describe the magnitude of the difference between these distributions. An obvious one is the difference in the means:

difference_in_means = male_sample.mean() - female_sample.mean()
difference_in_means # 單位為厘米

# 14.679014391669767

平均而言沼侣，男性高出 14 至 15 厘米颤诀。對于某些應(yīng)用象踊，這將是描述差異的好方法侍郭，但存在一些問題：

• 如果不了解更多關(guān)于分布的信息（比如標(biāo)準(zhǔn)差），很難解釋 15 厘米之間的差異是否很大齿尽。

• 差異的大小取決于度量單位沽损，因此很難在不同的研究中進(jìn)行比較。

有許多方法可以量化分布之間的差異循头。 一個簡單的選擇是將差異表示為平均值的百分比绵估。

# 練習(xí)：均值的相對差異，表示成百分比是什么卡骂？

relative_difference = difference_in_means / male_sample.mean()
relative_difference * 100   # 百分比
# 8.2389946286916569

但是相對差異的一個問題是你必須選擇国裳，相對于哪個均值來表達(dá)它們。

relative_difference = difference_in_means / female_sample.mean()
relative_difference * 100    # 百分比
# 8.9787536605040401

第二部分

表達(dá)分布之間差異的另一種方法是查看它們重疊的程度全跨。為了定義重疊缝左，我們在兩個均值之間選擇一個閾值。簡單的閾值是均值之間的中點(diǎn)：

simple_thresh = (mean1 + mean2) / 2
simple_thresh

# 170.82560946234622

一個更好但更復(fù)雜的閾值是 PDF 交叉的地方。

thresh = (std1 * mean2 + std2 * mean1) / (std1 + std2)
thresh

# 170.6040359174722

在此示例中渺杉，兩個閾值之間沒有太大差異蛇数。

現(xiàn)在我們可以算出有多少男性低于閾值：

male_below_thresh = sum(male_sample < thresh)
male_below_thresh

# 164

有多少女性高于它：

female_above_thresh = sum(female_sample > thresh)
female_above_thresh

# 174

“重疊”是曲線下面的總面積，最終位于閾值的右側(cè)是越。

overlap = male_below_thresh / len(male_sample) + female_above_thresh / len(female_sample)
overlap

# 0.33799999999999997

或者在更實(shí)際的術(shù)語中耳舅，如果你嘗試使用高度來猜測性別，你可能會報告錯誤分類的人數(shù)：

misclassification_rate = overlap / 2
misclassification_rate

# 0.16899999999999998

量化分布之間差異的另一種方法倚评，是所謂的“優(yōu)勢概率”挽放，這是一個有問題的術(shù)語，但在這種情況下蔓纠，隨機(jī)選擇的男性比隨機(jī)選擇的女性更高的概率。

# 練習(xí)：假設(shè)我隨機(jī)選擇男性和女性吗蚌。
# 男性更高的概率是什么腿倚？
sum(x > y for x, y in zip(male_sample, female_sample)) / len(male_sample)

# 0.91100000000000003

重疊（或錯誤分類率）和“優(yōu)勢概率”有兩個好的屬性：

• 作為概率，它們不依賴于度量單位蚯妇，因此它們在研究之間具有可比性敷燎。

• 它們以操作術(shù)語表達(dá)，因此讀者可以了解差異所產(chǎn)生的實(shí)際效果箩言。

還有另一種表達(dá)分布差異的常用方法硬贯。科恩的d是均值的差陨收，通過除以標(biāo)準(zhǔn)差來標(biāo)準(zhǔn)化饭豹。這是一個計算它的函數(shù)：

def CohenEffectSize(group1, group2):
    """計算 Cohen d。

    group1: Series 或 NumPy 數(shù)組
    group2: Series 或 NumPy 數(shù)組

    返回值: float
    """
    diff = group1.mean() - group2.mean()

    n1, n2 = len(group1), len(group2)
    var1 = group1.var()
    var2 = group2.var()

    pooled_var = (n1 * var1 + n2 * var2) / (n1 + n2)
    d = diff / numpy.sqrt(pooled_var)
    return d

計算分母有點(diǎn)復(fù)雜; 事實(shí)上务漩，人們提出了幾種方法來做到這一點(diǎn)拄衰。 該實(shí)現(xiàn)使用“池化標(biāo)準(zhǔn)差”，其是兩組標(biāo)準(zhǔn)差的加權(quán)平均值饵骨。

這是男女之間身高差異的結(jié)果翘悉。

CohenEffectSize(male_sample, female_sample)

1.9274780043619493

大多數(shù)人都不太了解d = 1.9是多么大，所以讓我們做一個可視化來校準(zhǔn)居触。

這是一個函數(shù)妖混，它封裝了我們已經(jīng)看到的用于計算重疊和優(yōu)勢概率的代碼。

def overlap_superiority(control, treatment, n=1000):
    """基于樣本估計重疊和有時轮洋。
    
    control: scipy.stats rv 對象
    treatment: scipy.stats rv 對象
    n: 樣本大小
    """
    control_sample = control.rvs(n)
    treatment_sample = treatment.rvs(n)
    thresh = (control.mean() + treatment.mean()) / 2
    
    control_above = sum(control_sample > thresh)
    treatment_below = sum(treatment_sample < thresh)
    overlap = (control_above + treatment_below) / n
    
    superiority = sum(x > y for x, y in zip(treatment_sample, control_sample)) / n
    return overlap, superiority

這是使用 Cohen 的d的函數(shù)制市，繪制具有給定效應(yīng)量的正態(tài)分布，并打印它們的重疊和優(yōu)勢砖瞧。

def plot_pdfs(cohen_d=2):
    """為標(biāo)準(zhǔn)差不同的分布繪制 PDF息堂。
    
    cohen_d: 均值之間的標(biāo)準(zhǔn)差
    """
    control = scipy.stats.norm(0, 1)
    treatment = scipy.stats.norm(cohen_d, 1)
    xs, ys = eval_pdf(control)
    pyplot.fill_between(xs, ys, label='control', color=COLOR3, alpha=0.7)

    xs, ys = eval_pdf(treatment)
    pyplot.fill_between(xs, ys, label='treatment', color=COLOR2, alpha=0.7)
    
    o, s = overlap_superiority(control, treatment)
    print('overlap', o)
    print('superiority', s)

這是一個演示函數(shù)的示例：

plot_pdfs(2)

'''
overlap 0.278
superiority 0.932
'''

png

你可以使用交互式組件來顯示d的不同值的含義：

slider = widgets.FloatSliderWidget(min=0, max=4, value=2)
interact(plot_pdfs, cohen_d=slider)
None

'''
overlap 0.305
superiority 0.931
'''

png

Cohen 的d有一些不錯的屬性：

• 因?yàn)槠骄岛蜆?biāo)準(zhǔn)差具有相同的單位，它們的比例是無量綱的，所以我們可以比較不同研究中的d荣堰。

• 在通常使用d的字段中床未，人們會進(jìn)行校準(zhǔn)，來了解哪些值應(yīng)該被認(rèn)為是大的振坚，令人驚訝的或重要的薇搁。

• 給定d（并假設(shè)分布是正態(tài)），你可以計算重疊渡八，優(yōu)勢和相關(guān)統(tǒng)計量啃洋。

總之，報告效應(yīng)量的最佳方式通常取決于受眾和你的目標(biāo)屎鳍。通常在具有良好技術(shù)屬性的摘要統(tǒng)計量宏娄，和對一般受眾有意義的統(tǒng)計量之間進(jìn)行權(quán)衡。

6.2 隨機(jī)采樣

署名：派生于 Allen Downey 的 CompStats逮壁。協(xié)議：Creative Commons Attribution 4.0 International孵坚。

from __future__ import print_function, division

import numpy
import scipy.stats

import matplotlib.pyplot as pyplot

from IPython.html.widgets import interact, fixed
from IPython.html import widgets

# 為隨機(jī)數(shù)生成器播種，所以我們得到相同結(jié)果
numpy.random.seed(18)

# 來自 http://colorbrewer2.org/ 一些漂亮的顏色
COLOR1 = '#7fc97f'
COLOR2 = '#beaed4'
COLOR3 = '#fdc086'
COLOR4 = '#ffff99'
COLOR5 = '#386cb0'

%matplotlib inline

第一部分

假設(shè)我們想估計美國男女的平均體重窥淆。我們希望量化估計的不確定性卖宠。一種方法是模擬許多實(shí)驗(yàn)，并觀察從一個實(shí)驗(yàn)到下一個實(shí)驗(yàn)的結(jié)果變化有多大忧饭。

我將從不切實(shí)際的假設(shè)開始扛伍，即我們知道總體中體重的實(shí)際分布。然后我將展示如何在沒有這個假設(shè)的情況下解決問題词裤。

根據(jù) BRFSS 的數(shù)據(jù)刺洒，我發(fā)現(xiàn)美國女性的體重（千克）分布，可以由對數(shù)正態(tài)分布模擬得很好亚斋，其參數(shù)如下：

weight = scipy.stats.lognorm(0.23, 0, 70.8)
weight.mean(), weight.std()

# (72.697645732966876, 16.944043048498038)

這是分布的樣子：

xs = numpy.linspace(20, 160, 100)
ys = weight.pdf(xs)
pyplot.plot(xs, ys, linewidth=4, color=COLOR1)
pyplot.xlabel('weight (kg)')
pyplot.ylabel('PDF')
None

png

make_sample從此分布中抽取隨機(jī)樣本作媚。 結(jié)果是 NumPy 數(shù)組。

def make_sample(n=100):
    sample = weight.rvs(n)
    return sample

這是一個n = 100的例子帅刊。 樣本的均值和標(biāo)準(zhǔn)差接近總體均值和標(biāo)準(zhǔn)差纸泡，但不準(zhǔn)確。

sample = make_sample(n=100)
sample.mean(), sample.std()

# (76.308293640077437, 19.995558735561865)

我們想估計總體中的平均體重赖瞒，因此我們將使用的“樣本統(tǒng)計量”是均值：

def sample_stat(sample):
    return sample.mean()

“實(shí)驗(yàn)”的一次迭代是收集 100 名女性的樣本并計算他們的平均體重女揭。我們可以模擬多次運(yùn)行此實(shí)驗(yàn)，并收集樣本統(tǒng)計量的列表栏饮。 結(jié)果是NumPy數(shù)組吧兔。

def compute_sample_statistics(n=100, iters=1000):
    stats = [sample_stat(make_sample(n)) for i in range(iters)]
    return numpy.array(stats)

下一行運(yùn)行模擬 1000 次并將結(jié)果放在sample_means中：

sample_means = compute_sample_statistics(n=100, iters=1000)

讓我們看一下樣本均值的分布。 此分布顯示了從一個實(shí)驗(yàn)到下一個實(shí)驗(yàn)的結(jié)果有多大差異袍嬉。

請記住境蔼，這種分布與總體中的體重分布不同灶平。這是重復(fù)虛擬實(shí)驗(yàn)的結(jié)果分布。

pyplot.hist(sample_means, color=COLOR5)
pyplot.xlabel('sample mean (n=100)')
pyplot.ylabel('count')
None

png

樣本均值接近實(shí)際總體均值箍土，這很好逢享，但實(shí)際上并不重要。

sample_means.mean()

# 72.652052080657413

樣本的標(biāo)準(zhǔn)差量化了從一個實(shí)驗(yàn)到下一個實(shí)驗(yàn)的可變性吴藻，并反映了估計的精確度瞒爬。該數(shù)量稱為“標(biāo)準(zhǔn)誤差”。

std_err = sample_means.std()
std_err

# 1.6355262477017491

我們還可以使用樣本均值的分布來計算“90% 置信區(qū)間”沟堡，其中包含 90% 的實(shí)驗(yàn)結(jié)果：

conf_int = numpy.percentile(sample_means, [5, 95])
conf_int

# array([ 69.92149384,  75.40866638])

以下函數(shù)接受一組樣本統(tǒng)計量并打印 SE 和 CI：

def summarize_sampling_distribution(sample_stats):
    print('SE', sample_stats.std())
    print('90% CI', numpy.percentile(sample_stats, [5, 95]))

這里是它的樣子：

summarize_sampling_distribution(sample_means)

'''
SE 1.6355262477
90% CI [ 69.92149384  75.40866638]
'''

現(xiàn)在我們想看看當(dāng)我們改變樣本大小n時會發(fā)生什么侧但。以下函數(shù)接受n，運(yùn)行 1000 個模擬實(shí)驗(yàn)航罗，并匯總結(jié)果禀横。

def plot_sample_stats(n, xlim=None):
    sample_stats = compute_sample_statistics(n, iters=1000)
    summarize_sampling_distribution(sample_stats)
    pyplot.hist(sample_stats, color=COLOR2)
    pyplot.xlabel('sample statistic')
    pyplot.xlim(xlim)

這是n = 100的測試：

plot_sample_stats(100)

'''
SE 1.71202891175
90% CI [ 69.96057332  75.58582662]
'''

png

現(xiàn)在我們可以使用interact和不同的n值來運(yùn)行plot_sample_stats。注意：xlim設(shè)置x軸的限制粥血。

def sample_stat(sample):
    return sample.mean()

slider = widgets.IntSliderWidget(min=10, max=1000, value=100)
interact(plot_sample_stats, n=slider, xlim=fixed([55, 95]))
None

'''
SE 1.71314776815
90% CI [ 69.99274896  75.65943185]
'''

png

此框架適用于我們想要估計的任何其他數(shù)量燕侠。通過更改sample_stat，你可以計算任何樣本統(tǒng)計量的 SE 和 CI立莉。

作為練習(xí)，請使用以下任何統(tǒng)計量填寫下面的sample_stat：

• 樣本標(biāo)準(zhǔn)差
• 變異系數(shù)七问，即樣本標(biāo)準(zhǔn)差除以樣本標(biāo)準(zhǔn)均值蜓耻。
• 最小值或最大值
• 中位數(shù)（第 50 個百分位數(shù)）
• 第 10 或 90 個百分位數(shù)
• 四分位數(shù)間距（IQR），即第 75 和第 25 百分位數(shù)之間的差械巡。

你可能會發(fā)現(xiàn)有用的 NumPy 數(shù)組方法包括std刹淌，min，max和percentile讥耗。根據(jù)結(jié)果有勾，你可能需要調(diào)整xlim。

def sample_stat(sample):
    # TODO：將下面一行替換為其它樣本統(tǒng)計量
    return sample.mean()

slider = widgets.IntSliderWidget(min=10, max=1000, value=100)
interact(plot_sample_stats, n=slider, xlim=fixed([0, 100]))
None

'''
SE 1.67195986148
90% CI [ 69.82954731  75.32184298]
'''

png

第二部分

到目前為止古程，我們已經(jīng)證明蔼卡，如果我們知道總體的實(shí)際分布，我們可以計算任何樣本統(tǒng)計量的抽樣分布挣磨，從中我們可以計算 SE 和 CI雇逞。

但在現(xiàn)實(shí)生活中，我們并不知道總體的實(shí)際分布茁裙。 如果我們這樣做塘砸，我們就不需要估計了！

在現(xiàn)實(shí)生活中晤锥，我們使用樣本建立總體分布的模型掉蔬，然后使用該模型生成抽樣分布廊宪。一種簡單而流行的方法是“重采樣”，這意味著我們將樣本本身用作總體分布的模型并從中抽取樣本女轿。

在繼續(xù)之前箭启，我想收集第一部分中的一些代碼并將其組織為一個類。 此類表示用于計算采樣分布的框架谈喳。

class Resampler(object):
    """表示計算采樣分布的框架册烈。"""
    
    def __init__(self, sample, xlim=None):
        """儲存實(shí)際樣本。"""
        self.sample = sample
        self.n = len(sample)
        self.xlim = xlim
        
    def resample(self):
        """通過帶放回地選擇原始樣本生成新樣本婿禽。
        """
        new_sample = numpy.random.choice(self.sample, self.n, replace=True)
        return new_sample
    
    def sample_stat(self, sample):
        """計算樣本統(tǒng)計量赏僧，使用原始樣本或者模擬樣本。
        """
        return sample.mean()
    
    def compute_sample_statistics(self, iters=1000):
        """模擬許多實(shí)驗(yàn)并收集所得樣本統(tǒng)計量扭倾。
        """
        stats = [self.sample_stat(self.resample()) for i in range(iters)]
        return numpy.array(stats)
    
    def plot_sample_stats(self):
        """運(yùn)行模擬的實(shí)驗(yàn)淀零，并匯總結(jié)果。
        """
        sample_stats = self.compute_sample_statistics()
        summarize_sampling_distribution(sample_stats)
        pyplot.hist(sample_stats, color=COLOR2)
        pyplot.xlabel('sample statistic')
        pyplot.xlim(self.xlim)

以下函數(shù)實(shí)例化一個Resampler并運(yùn)行它膛壹。

def plot_resampled_stats(n=100):
    sample = weight.rvs(n)
    resampler = Resampler(sample, xlim=[55, 95])
    resampler.plot_sample_stats()

這是一個n = 100的測試：

plot_resampled_stats(100)

'''
SE 1.72606450921
90% CI [ 71.35648645  76.82647135]
'''

png

現(xiàn)在我們可以在交互中使用plot_resampled_stats：

slider = widgets.IntSliderWidget(min=10, max=1000, value=100)
interact(plot_resampled_stats, n=slider, xlim=fixed([1, 15]))
None

'''
SE 1.67407589545
90% CI [ 69.60129748  75.13161693]
'''

png

練習(xí)：編寫一個名為StdResampler的新類驾中，它繼承自Resampler并覆蓋sample_stat，因此它計算重采樣數(shù)據(jù)的標(biāo)準(zhǔn)差模聋。

class StdResampler(Resampler):   
    """計算標(biāo)準(zhǔn)差的采樣分布肩民。"""
    
    def sample_stat(self, sample):
        """計算樣本統(tǒng)計量，使用原始樣本或者模擬樣本链方。
        """
        return sample.std()

使用下面的單元格測試你的代碼：

def plot_resampled_stats(n=100):
    sample = weight.rvs(n)
    resampler = StdResampler(sample, xlim=[0, 100])
    resampler.plot_sample_stats()
    
plot_resampled_stats()

'''
SE 1.30056137605
90% CI [ 13.70615766  18.05008376]
'''

png

當(dāng)你的StdResampler能用時持痰，你應(yīng)該能夠與它進(jìn)行交互：

slider = widgets.IntSliderWidget(min=10, max=1000, value=100)
interact(plot_resampled_stats, n=slider)
None

'''
SE 1.29095098626
90% CI [ 15.13442137  19.27452588]
'''

png

第三部分

我們可以擴(kuò)展這個框架來計算 SE 和 CI，來獲得均值的差祟蚀。例如工窍，男性平均比女性重。以下是女性的分布（來自 BRFSS 數(shù)據(jù)）：

female_weight = scipy.stats.lognorm(0.23, 0, 70.8)
female_weight.mean(), female_weight.std()

# (72.697645732966876, 16.944043048498038)

這是男性的分布：

male_weight = scipy.stats.lognorm(0.20, 0, 87.3)
male_weight.mean(), male_weight.std()

# (89.063576984335782, 17.992335889366288)

我將模擬 100 名男性和 100 名女性的樣本：

female_sample = female_weight.rvs(100)
male_sample = male_weight.rvs(100)

平均值的差應(yīng)為 17 千克左右前酿，但從一個隨機(jī)樣本到另一個會有所不同：

male_sample.mean() - female_sample.mean()

# 15.521337537414979

這是計算 Cohen 的d的函數(shù)：

def CohenEffectSize(group1, group2):
    """計算 Cohen d患雏。

    group1: Series 或 NumPy 數(shù)組
    group2: Series 或 NumPy 數(shù)組

    返回值: float
    """
    diff = group1.mean() - group2.mean()

    n1, n2 = len(group1), len(group2)
    var1 = group1.var()
    var2 = group2.var()

    pooled_var = (n1 * var1 + n2 * var2) / (n1 + n2)
    d = diff / numpy.sqrt(pooled_var)
    return d

男女之間的體重差異約為 1 個標(biāo)準(zhǔn)差：

CohenEffectSize(male_sample, female_sample)

# 0.9271315108152719

現(xiàn)在我們可以編寫一個版本的Resampler來計算d的采樣分布。

class CohenResampler(Resampler):
    def __init__(self, group1, group2, xlim=None):
        self.group1 = group1
        self.group2 = group2
        self.xlim = xlim
        
    def resample(self):
        group1 = numpy.random.choice(self.group1, len(self.group1), replace=True)
        group2 = numpy.random.choice(self.group2, len(self.group2), replace=True)
        return group1, group2
    
    def sample_stat(self, groups):
        group1, group2 = groups
        return CohenEffectSize(group1, group2)
    
    # 注：下面的函數(shù)和 Resampler 中相同罢维，
    # 所以我可以僅僅繼承它淹仑，但是我為了可讀性而包含它
    def compute_sample_statistics(self, iters=1000):
        stats = [self.sample_stat(self.resample()) for i in range(iters)]
        return numpy.array(stats)
    
    def plot_sample_stats(self):
        sample_stats = self.compute_sample_statistics()
        summarize_sampling_distribution(sample_stats)
        pyplot.hist(sample_stats, color=COLOR2)
        pyplot.xlabel('sample statistic')
        pyplot.xlim(self.xlim)

現(xiàn)在我們可以實(shí)例化一個CohenResampler并繪制采樣分布。

resampler = CohenResampler(male_sample, female_sample)
resampler.plot_sample_stats()

'''
SE 0.160707033098
90% CI [ 0.6808076  1.1974013]
'''

png

該示例展示了計算框架優(yōu)于數(shù)學(xué)分析的優(yōu)點(diǎn)肺孵。Cohen 的d等統(tǒng)計量是其他統(tǒng)計數(shù)據(jù)的比率攻人，相對難以分析。 但是通過計算方法悬槽，所有樣本統(tǒng)計量都同樣“容易”怀吻。

關(guān)于詞匯的一個注解：我在這里稱之為“重采樣”的東西，是一種稱為“自舉”的特定重采樣初婆。其他技術(shù)也考慮重新采樣蓬坡，包括置換檢驗(yàn)猿棉，我們將在下一節(jié)中看到，以及“jackknife”重采樣屑咳。 你可以在 http://en.wikipedia.org/wiki/Resampling_(statistics) 上閱讀更多內(nèi)容萨赁。

6.3 假設(shè)檢驗(yàn)

署名：派生于 Allen Downey 的 CompStats。協(xié)議：Creative Commons Attribution 4.0 International兆龙。

from __future__ import print_function, division

import numpy
import scipy.stats

import matplotlib.pyplot as pyplot

from IPython.html.widgets import interact, fixed
from IPython.html import widgets

import first

# 為隨機(jī)數(shù)生成器播種杖爽，所以我們得到相同結(jié)果
numpy.random.seed(19)

# 來自 http://colorbrewer2.org/ 的一些漂亮的顏色
COLOR1 = '#7fc97f'
COLOR2 = '#beaed4'
COLOR3 = '#fdc086'
COLOR4 = '#ffff99'
COLOR5 = '#386cb0'

%matplotlib inline

第一部分

作為一個例子，讓我們看看分組之間的差異紫皇。 我在 Think Stats 中使用的例子是與其他嬰兒相比的第一個嬰兒慰安。first模塊提供代碼，將數(shù)據(jù)讀入三個 pandas 數(shù)據(jù)幀聪铺。

live, firsts, others = first.MakeFrames()

我們感興趣的表觀效應(yīng)是均值的差異化焕。其他示例可能包括變量之間的相關(guān)性或線性回歸中的系數(shù)。量化效應(yīng)量的數(shù)字铃剔，無論它是什么撒桨，都是“測試統(tǒng)計量”。

def TestStatistic(data):
    group1, group2 = data
    test_stat = abs(group1.mean() - group2.mean())
    return test_stat

對于第一個例子键兜，我提取了第一個嬰兒和其他人的懷孕時間凤类。結(jié)果是 pandas Series對象。

group1 = firsts.prglngth
group2 = others.prglngth

平均值的實(shí)際差異為 0.078 周普气，僅為13小時踱蠢。

actual = TestStatistic((group1, group2))
actual

# 0.078037266777549519

零假設(shè)是組之間沒有差異。我們可以通過形成包括第一個嬰兒和其他嬰兒的合并樣本來對其進(jìn)行建模棋电。

n, m = len(group1), len(group2)
pool = numpy.hstack((group1, group2))

然后我們可以通過打亂池子，并將其分成兩組來模擬零假設(shè)苇侵，使用與實(shí)際樣本相同的大小赶盔。

def RunModel():
    numpy.random.shuffle(pool)
    data = pool[:n], pool[n:]
    return data

運(yùn)行該模型的結(jié)果是兩個 NumPy 數(shù)組，其具有打亂的孕期長度：

RunModel()

# (array([36, 40, 39, ..., 43, 42, 40]), array([43, 39, 32, ..., 37, 35, 41]))

然后我們使用模擬數(shù)據(jù)計算相同的測試統(tǒng)計量：

TestStatistic(RunModel())

# 0.081758440969863955

如果我們運(yùn)行模型 1000 次并計算測試統(tǒng)計量榆浓，我們可以看到測試統(tǒng)計量在零假設(shè)下變化了多少于未。

test_stats = numpy.array([TestStatistic(RunModel()) for i in range(1000)])
test_stats.shape

# (1000,)

這是零假設(shè)下的檢驗(yàn)統(tǒng)計量的抽樣分布，其中均值的實(shí)際差異用灰線表示陡鹃。

def VertLine(x):
    """在 x 處繪制豎直線烘浦。"""
    pyplot.plot([x, x], [0, 300], linewidth=3, color='0.8')

VertLine(actual)
pyplot.hist(test_stats, color=COLOR5)
pyplot.xlabel('difference in means')
pyplot.ylabel('count')
None

png

p 值是零假設(shè)下的檢驗(yàn)統(tǒng)計量超過實(shí)際值的概率。

pvalue = sum(test_stats >= actual) / len(test_stats)
pvalue

# 0.14999999999999999

在這種情況下萍鲸，結(jié)果約為 15%闷叉，這意味著即使兩組之間沒有差異，我們也可以看到樣本差異大到 0.078 周脊阴。

我們的結(jié)論是握侧，表觀效應(yīng)可能是偶然的蚯瞧，所以我們不相信它會出現(xiàn)在一般總體或同一總體的另一個樣本中。

第二部分

我們可以從上一節(jié)中獲取部分品擎，并將它們組織在一個表示假設(shè)檢驗(yàn)結(jié)構(gòu)的類中埋合。

class HypothesisTest(object):
    """表示假設(shè)檢驗(yàn)。"""

    def __init__(self, data):
        """初始化萄传。

        data: 任何格式相關(guān)的數(shù)據(jù)
        """
        self.data = data
        self.MakeModel()
        self.actual = self.TestStatistic(data)
        self.test_stats = None

    def PValue(self, iters=1000):
        """計算測試統(tǒng)計量的分布和 p 值甚颂。

        iters: 迭代數(shù)

        返回值: float p 值
        """
        self.test_stats = numpy.array([self.TestStatistic(self.RunModel()) 
                                       for _ in range(iters)])

        count = sum(self.test_stats >= self.actual)
        return count / iters

    def MaxTestStat(self):
        """返回模擬期間見到的最大的測試統(tǒng)計量。
        """
        return max(self.test_stats)

    def PlotHist(self, label=None):
        """使用觀測測試統(tǒng)計量處的豎直線條繪制 Cdf秀菱。
        """
        def VertLine(x):
            """在 x 處繪制豎直線條振诬。"""
            pyplot.plot([x, x], [0, max(ys)], linewidth=3, color='0.8')

        ys, xs, patches = pyplot.hist(ht.test_stats, color=COLOR4)
        VertLine(self.actual)
        pyplot.xlabel('test statistic')
        pyplot.ylabel('count')

    def TestStatistic(self, data):
        """計算測試統(tǒng)計量。

        data: 任何格式相關(guān)的數(shù)據(jù)       
        """
        raise UnimplementedMethodException()

    def MakeModel(self):
        """為零假設(shè)構(gòu)建模型
        """
        pass

    def RunModel(self):
        """運(yùn)行零假設(shè)的模型答朋。

        返回值: 模擬的數(shù)據(jù)
        """
        raise UnimplementedMethodException()

HypothesisTest是一個對模板進(jìn)行編碼的抽象父類贷揽。子類填寫缺少的方法。 例如梦碗，這是上一節(jié)的測試禽绪。

class DiffMeansPermute(HypothesisTest):
    """通過置換測試均值差異。"""

    def TestStatistic(self, data):
        """計算測試統(tǒng)計量.

        data: 任何格式相關(guān)的數(shù)據(jù)       
        """
        group1, group2 = data
        test_stat = abs(group1.mean() - group2.mean())
        return test_stat

    def MakeModel(self):
        """構(gòu)建零假設(shè)的模型洪规。
        """
        group1, group2 = self.data
        self.n, self.m = len(group1), len(group2)
        self.pool = numpy.hstack((group1, group2))

    def RunModel(self):
        """運(yùn)行零假設(shè)的模型印屁。

        返回值: 模擬的數(shù)據(jù)
        """
        numpy.random.shuffle(self.pool)
        data = self.pool[:self.n], self.pool[self.n:]
        return data

現(xiàn)在我們可以通過實(shí)例化DiffMeansPermute對象來運(yùn)行測試：

data = (firsts.prglngth, others.prglngth)
ht = DiffMeansPermute(data)
p_value = ht.PValue(iters=1000)
print('\nmeans permute pregnancy length')
print('p-value =', p_value)
print('actual =', ht.actual)
print('ts max =', ht.MaxTestStat())

means permute pregnancy length
p-value = 0.16
actual = 0.0780372667775
ts max = 0.173695697482

我們可以在零假設(shè)下繪制檢驗(yàn)統(tǒng)計量的抽樣分布。

ht.PlotHist()

png

作為練習(xí)斩例，編寫一個名為DiffStdPermute的類雄人，它擴(kuò)展了DiffMeansPermute并覆蓋TestStatistic來計算標(biāo)準(zhǔn)差的差異。標(biāo)準(zhǔn)差的差異是否具有統(tǒng)計學(xué)意義念赶？

class DiffStdPermute(DiffMeansPermute):
    """通過置換測試均值差異础钠。"""

    def TestStatistic(self, data):
        """計算測試統(tǒng)計量。

        data: 任何格式相關(guān)的數(shù)據(jù)
        """
        group1, group2 = data
        test_stat = abs(group1.std() - group2.std())
        return test_stat

data = (firsts.prglngth, others.prglngth)
ht = DiffStdPermute(data)
p_value = ht.PValue(iters=1000)
print('\nstd permute pregnancy length')
print('p-value =', p_value)
print('actual =', ht.actual)
print('ts max =', ht.MaxTestStat())

'''
std permute pregnancy length
p-value = 0.155
actual = 0.176049064229
ts max = 0.44299505029
'''

現(xiàn)在讓我們再次運(yùn)行DiffMeansPermute叉谜，看看第一胎和其他嬰兒的出生體重是否有差異旗吁。

data = (firsts.totalwgt_lb.dropna(), others.totalwgt_lb.dropna())
ht = DiffMeansPermute(data)
p_value = ht.PValue(iters=1000)
print('\nmeans permute birthweight')
print('p-value =', p_value)
print('actual =', ht.actual)
print('ts max =', ht.MaxTestStat())

'''
means permute birthweight
p-value = 0.0
actual = 0.124761184535
ts max = 0.0917504268392
'''

在這種情況下，在 1000 次嘗試之后停局，我們從未看到與觀察到的差異一樣大的樣本差異很钓，因此我們得出結(jié)論，表觀效應(yīng)不太可能在零假設(shè)下董栽。在正常情況下码倦，我們也可以推斷出表觀效應(yīng)不太可能是由隨機(jī)抽樣引起的。

最后一點(diǎn)：在這種情況下锭碳，我會報告p值小于 1/1000 或 0.001袁稽。 我不會報告p = 0，因?yàn)樵诹慵僭O(shè)下擒抛，標(biāo)貫效應(yīng)并非不可能运提，而是不太可能蝗柔。

這個筆記本由 Donne Martin 編寫。來源和協(xié)議信息在 GitHub 上民泵。