Think Stats kitabının üçüncü baskısını Bookshop.org ve Amazon üzerinden edinebilirsiniz (bunlar satış ortaklığı bağlantılarıdır). Ücretsiz çevrimiçi sürümü beğeniyorsanız bana bir kahve ısmarlamayı düşünebilirsiniz.
14. Analitik yöntemler#
Bu kitap simülasyon ve yeniden örnekleme gibi hesaplamalı yöntemlere odaklandı; ancak çözdüğümüz bazı problemlerin çok daha hızlı hesaplanabilen analitik çözümleri vardır.
Bu bölüm bu yöntemlerden bazılarını sunup nasıl çalıştıklarını açıklar. Sonunda veri analizinde hesaplamalı ve analitik yöntemleri bütünleştirmek için öneriler veriyorum.
Bu notebook’u Colab’da çalıştırmak için buraya tıklayın.
14.1. Normal olasılık grafikleri#
Birçok analitik yöntem normal dağılımın özelliklerine dayanır: gerçek dünyadaki pek çok ölçüm normal dağılımla iyi yaklaşıklandırılır ve normal dağılımların analize yararlı matematiksel özellikleri vardır.
İlk noktayı penguen ölçümleriyle göstereceğiz, ardından matematiksel özellikleri inceleyeceğiz. İndirme yönergeleri notebook’tadır.
Aşağıdaki hücre verileri Allison Horst’un deposundan indirir.
Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R paketi 0.1.0. https://allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218.
Veriler şu araştırmada toplanmıştır: Gorman KB, Williams TD, Fraser WR (2014). Ecological sexual dimorphism and environmental variability within a community of Antarctic penguins (genus Pygoscelis). PLoS ONE 9(3):e90081. https://doi.org/10.1371/journal.pone.0090081
download(
"https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)
Verileri şöyle okuyabiliriz.
penguins = pd.read_csv("penguins_raw.csv")
penguins.shape
(344, 17)
Veri kümesi üç penguen türünün ölçümlerini içerir. Bu örnekte Adélie penguenlerini seçeceğiz.
adelie = penguins.query('Species.str.startswith("Adelie")').copy()
len(adelie)
152
Ağırlıkların normal dağılım izleyip izlemediğini görmek için verinin ampirik CDF’sini hesaplayacağız.
from empiricaldist import Cdf
weights = adelie["Body Mass (g)"].dropna()
cdf_weights = Cdf.from_seq(weights)
Aynı ortalama ve standart sapmaya sahip normal dağılımın analitik CDF’sini de hesaplayacağız.
m, s = weights.mean(), weights.std()
m, s
(np.float64(3700.662251655629), np.float64(458.5661259101348))
from scipy.stats import norm
dist = norm(m, s)
qs = np.linspace(m - 3.5 * s, m + 3.5 * s)
ps = dist.cdf(qs)
Verinin CDF’sinin normal modelle karşılaştırması şöyledir.
model_options = dict(color="gray", alpha=0.5, label="model")
plt.plot(qs, ps, **model_options)
cdf_weights.plot(label="veri")
decorate(ylabel="CDF")
Normal dağılım bu veri için yeterince iyi bir model olabilir, ancak kusursuz değildir.
Veri ve model CDF’lerini çizmek uyumu değerlendirmenin iyi yoludur; dezavantajı model parametrelerini, burada ortalama ile standart sapmayı, ne kadar iyi tahmin ettiğimize bağlı olmasıdır.
Başka bir seçenek parametre tahminine bağlı olmayan normal olasılık grafiğidir. Bu grafikte \(y\) değerleri sıralanmış ölçümlerdir.
ys = np.sort(weights)
\(x\) değerleri ise normal dağılımın karşılık gelen yüzdelikleridir; bunları ters CDF olan “percent point function”ı hesaplayan norm.ppf ile buluruz.
n = len(weights)
ps = (np.arange(n) + 0.5) / n
xs = norm.ppf(ps)
Ölçümler gerçekten normal dağılımdan geliyorsa \(x\) ve \(y\) değerleri doğru üzerinde olmalıdır. Ne kadar yaklaştıklarını görmek için linregress ile doğru uydurabiliriz.
from scipy.stats import linregress
results = linregress(xs, ys)
intercept, slope = results.intercept, results.slope
fit_xs = np.linspace(-3, 3)
fit_ys = intercept + slope * fit_xs
Aşağıdaki şekil \(x\) ve \(y\) değerleriyle uydurulan doğruyu gösterir.
plt.plot(fit_xs, fit_ys, **model_options)
plt.plot(xs, ys, label="veri")
decorate(xlabel="Standart normal", ylabel="Vücut kütlesi (g)")
Grafik tam bir doğru oluşturmaz; bu da normal dağılımın kusursuz bir model olmadığını gösterir.
Nedenlerden biri veri kümesinde ortalamaları farklı erkek ve dişi penguenlerin birlikte bulunmasıdır. Grupları ayrı çizelim. Aşağıdaki işlev normal olasılık grafiğinin adımlarını kapsar.
def normal_probability_plot(sample, **options):
"""Uyum doğrusuyla normal olasılık grafiği oluşturur."""
n = len(sample)
ps = (np.arange(n) + 0.5) / n
xs = norm.ppf(ps)
ys = np.sort(sample)
results = linregress(xs, ys)
intercept, slope = results.intercept, results.slope
fit_xs = np.linspace(-3, 3)
fit_ys = intercept + slope * fit_xs
plt.plot(fit_xs, fit_ys, color="gray", alpha=0.5)
plt.plot(xs, ys, **options)
decorate(xlabel="Standart normal")
Erkek ve dişi penguenlerin ayrı sonuçları şöyledir.
grouped = adelie.groupby("Sex")
weights_male = grouped.get_group("MALE")["Body Mass (g)"]
normal_probability_plot(weights_male, ls="--", label="Erkek")
weights_female = grouped.get_group("FEMALE")["Body Mass (g)"]
normal_probability_plot(weights_female, label="Dişi")
decorate(ylabel="Ağırlık (g)")
İki grubun grafiği de doğruya yakındır; ağırlık dağılımlarının normal olduğunu gösterir. Grupları birleştirince farklı ortalamalı iki normal dağılımın karışımını elde ederiz; böyle bir karışım her zaman normal dağılımla iyi modellenmez.
Şimdi normal dağılımları analiz için yararlı kılan matematiksel özellikleri ele alalım.
14.2. Normal dağılımlar#
Aşağıdaki sınıf normal dağılımı temsil eden nesne tanımlar. mu ve sigma2 öznitelikleri ortalama ile varyansı gösterir. sigma2 adı varyansın, genellikle sigma ile gösterilen standart sapmanın karesi olduğunu anımsatır.
class Normal:
"""Normal dağılımı temsil eder"""
def __init__(self, mu, sigma2):
"""Normal nesnesi oluşturur.
mu: ortalama
sigma2: varyans
"""
self.mu = mu
self.sigma2 = sigma2
def __repr__(self):
"""Dize gösterimini döndürür."""
return f"Normal({self.mu}, {self.sigma2})"
__str__ = __repr__
Örnek olarak erkek penguen ağırlıklarıyla aynı ortalama ve varyansa sahip bir Normal nesnesi oluşturacağız.
m, s = weights_male.mean(), weights_male.std()
dist_male = Normal(m, s**2)
dist_male
Normal(4043.4931506849316, 120278.25342465754)
Dişi penguen ağırlıklarıyla aynı ortalama ve varyansa sahip başka bir Normal nesnesi de oluşturacağız.
m, s = weights_female.mean(), weights_female.std()
dist_female = Normal(m, s**2)
dist_female
Normal(3368.8356164383563, 72565.63926940637)
Normal dağılımdan rastgele örneklem üreten yöntemi Normal sınıfına ekleyeceğiz. Var olan sınıfa yöntem eklemek için thinkstats içinde tanımlı Jupyter sihirli komutu add_method_to kullanacağız. Bu Python’ın parçası değildir; yalnızca Jupyter notebook’larında çalışır.
%%add_method_to Normal
def sample(self, n):
"""Bu dağılımdan rastgele bir örneklem üretir."""
sigma = np.sqrt(self.sigma2)
return np.random.normal(self.mu, sigma, n)
sample ile normal dağılımın ilk yararlı özelliğini göstereceğiz: iki normal dağılımdan çekilen değerlerin toplamı da normaldir.
Oluşturduğumuz Normal nesnelerinden örneklemler üretip toplayacak ve toplamların normal olasılık grafiğini oluşturacağız.
sample_sum = dist_male.sample(1000) + dist_female.sample(1000)
normal_probability_plot(sample_sum)
decorate(ylabel="Toplam ağırlık (g)")
Grafik doğruya benzer; toplamların normal dağıldığını gösterir. İki dağılımın parametreleri biliniyorsa toplam dağılımının parametreleri de hesaplanabilir. Aşağıdaki yöntem bunu gösterir.
%%add_method_to Normal
def __add__(self, other):
"""İki normal dağılımın toplamının dağılımı."""
return Normal(self.mu + other.mu, self.sigma2 + other.sigma2)
Toplam dağılımının ortalaması ortalamaların, varyansı varyansların toplamıdır. __add__ özel yöntemini tanımladığımız için iki dağılımı + ile “toplayabilir”, yani toplamlarının dağılımını hesaplayabiliriz.
dist_sum = dist_male + dist_female
dist_sum
Normal(7412.328767123288, 192843.8926940639)
Sonucu doğrulamak için normal dağılımın analitik CDF’sini çizen aşağıdaki yöntemi kullanacağız.
%%add_method_to Normal
def plot_cdf(self, n_sigmas=3.5, **options):
"""Bu dağılımın CDF'sini çizer."""
mu, sigma = self.mu, np.sqrt(self.sigma2)
low, high = mu - n_sigmas * sigma, mu + n_sigmas * sigma
xs = np.linspace(low, high, 101)
ys = norm.cdf(xs, mu, sigma)
plt.plot(xs, ys, **options)
Rastgele örneklem toplamının ampirik CDF’siyle birlikte sonuç şöyledir.
dist_sum.plot_cdf(**model_options)
Cdf.from_seq(sample_sum).plot(label="örneklem")
decorate(xlabel="Toplam ağırlık (g)", ylabel="CDF")
Hesapladığımız parametreler doğru görünüyor; bu da iki normal dağılımı ortalama ve varyanslarını toplayarak birleştirebildiğimizi doğruluyor.
Sonuç olarak normal dağılımdan n değer üretip toplarsak toplamın dağılımı da normaldir. Göstermek için erkek ağırlığı dağılımından 73 değer üretip toplayacağız. Aşağıdaki döngü bunu 1001 kez yapar; sonuç toplamların dağılımından bir örneklemdir.
n = len(weights_male)
sample_sums_male = [dist_male.sample(n).sum() for i in range(1001)]
n
73
Aşağıdaki yöntem toplamların dağılımını temsil eden Normal nesnesi oluşturur. Parametreler için ortalama ve varyansı n ile çarparız.
%%add_method_to Normal
def sum(self, n):
"""n değerin toplamının dağılımını döndürür."""
return Normal(n * self.mu, n * self.sigma2)
n ağırlığın toplam dağılımı şöyledir.
dist_sums_male = dist_male.sum(n)
Rastgele örneklemin ampirik dağılımıyla karşılaştırması şöyledir.
dist_sums_male.plot_cdf(**model_options)
Cdf.from_seq(sample_sums_male).plot(label="örneklem")
decorate(xlabel="Toplam ağırlıklar (g)", ylabel="CDF")
Analitik dağılım örnekleme uyar; sum yöntemini doğrular. Dolayısıyla n ölçümlük örneklemin toplam dağılımını hesaplayabiliriz.
14.3. Örneklem ortalamalarının dağılımı#
Örneklem toplamının dağılımından örneklem ortalamasının dağılımını da bulabiliriz. Normal dağılımın üçüncü özelliğini kullanacağız: sabitle çarpılır ya da bölünürse sonuç normaldir. Aşağıdaki yöntemler çarpım ve bölüm dağılımının parametrelerini hesaplar.
%%add_method_to Normal
def __mul__(self, factor):
"""Multiplies by a scalar."""
return Normal(factor * self.mu, factor**2 * self.sigma2)
%%add_method_to Normal
def __truediv__(self, factor):
"""Divides by a scalar."""
return self * (1 / factor)
Çarpım dağılımında ortalamayı factor, varyansı factor karesiyle çarparız. Bu özellikle örneklem ortalamalarının dağılımını hesaplayabiliriz.
dist_mean_male = dist_sums_male / n
Sonucu doğrulamak için rastgele örneklemlerin ortalamalarını da hesaplayacağız.
sample_means_male = np.array(sample_sums_male) / n
Normal modeli örneklem ortalamalarının ampirik CDF’siyle karşılaştıracağız.
dist_mean_male.plot_cdf(**model_options)
Cdf.from_seq(sample_means_male).plot(label="örneklem")
decorate(xlabel="Ortalama ağırlık (g)", ylabel="CDF")
Model ile simülasyon uyuşur; örneklem ortalamalarının dağılımını yeniden örneklemeden çok daha hızlı biçimde analitik hesaplayabildiğimizi gösterir.
Ortalamanın örnekleme dağılımından, dağılımın standart sapması olan standart hatayı hesaplayabiliriz.
standard_error = np.sqrt(dist_mean_male.sigma2)
standard_error
np.float64(40.591222045992765)
Sonuç, örnekleme dağılımını oluşturmadan standart hatayı doğrudan hesaplama kısayolunu gösterir. Varyansı n ile çarpıp n**2’ye böldük; net etki varyansı n’ye, standart sapmayı n’nin kareköküne bölmektir.
Örneklem ortalamasının standart hatasını şöyle hesaplayabiliriz.
standard_error = weights_male.std() / np.sqrt(n)
standard_error
np.float64(40.59122204599277)
Şimdi normal dağılımlarla hesaplayabileceğimiz bir sonucu daha, farkların dağılımını ele alalım.
14.4. Farkların dağılımı#
Önceki adımları birleştirerek dişi penguen ağırlıklarının örneklem ortalaması dağılımını şöyle hesaplarız.
n = len(weights_female)
dist_mean_female = dist_female.sum(n) / n
dist_mean_female
Normal(3368.835616438356, 994.0498530055667)
Erkek ve dişi ortalama ağırlıklarının örnekleme dağılımları elimizde; farkların dağılımını hesaplayalım. Aşağıdaki yöntem iki normal dağılımdan gelen değerlerin fark dağılımını hesaplar.
%%add_method_to Normal
def __sub__(self, other):
"""Bir farkın dağılımını hesaplar."""
return Normal(self.mu - other.mu, self.sigma2 + other.sigma2)
Bekleyeceğiniz gibi farkların ortalaması ortalamaların farkıdır. Belki beklemeyeceğiniz biçimde varyans, varyansların farkı değil toplamıdır. Çıkarmayı iki adımda düşünün:
İkinci dağılımın işaretini ters çevirince ortalamanın işareti değişir, varyans aynı kalır.
İlk dağılımı ekleyince toplamın varyansı varyansların toplamıdır.
İkna olmadıysanız sınayalım. Farkların analitik dağılımı şöyledir.
dist_diff_means = dist_mean_male - dist_mean_female
dist_diff_means
Normal(674.6575342465753, 2641.697160192656)
Rastgele bir fark örneklemi de şöyledir.
sample_sums_female = [dist_female.sample(n).sum() for i in range(1001)]
sample_means_female = np.array(sample_sums_female) / n
sample_diff_means = sample_means_male - sample_means_female
Aşağıdaki şekil rastgele örneklemin ampirik CDF’sini ve normal dağılımın analitik CDF’sini gösterir.
dist_diff_means.plot_cdf(**model_options)
Cdf.from_seq(sample_diff_means).plot(label="örneklem")
decorate(xlabel="Ortalama ağırlık farkı (g)", ylabel="CDF")
Uyuşmaları fark dağılımını doğru bulduğumuzu onaylar. Bu dağılımla ağırlık farkının güven aralığını hesaplayabiliriz. Ters CDF için aşağıdaki yöntemi kullanacağız.
%%add_method_to Normal
def ppf(self, xs):
sigma = np.sqrt(self.sigma2)
return norm.ppf(xs, self.mu, sigma)
ve 95. yüzdelikler %90 güven aralığını oluşturur.
ci90 = dist_diff_means.ppf([0.05, 0.95])
ci90
array([590.1162635 , 759.19880499])
Rastgele örneklemden yaklaşık aynı sonuçları alırız.
np.percentile(sample_diff_means, [5, 95])
array([595.98498856, 761.53168013])
Analitik yöntem yeniden örneklemeden hızlı ve belirlenimcidir, yani rastgele değildir.
Ancak şimdiye kadarki her şey ölçümlerin normal dağıldığı varsayımına dayanır. Bu her zaman doğru değildir; gerçek veride hiçbir zaman tam doğru değildir. Yine de çok sayıda ölçümü toplarsak, ölçümlerin dağılımı normal olmasa bile toplamlarının dağılımı çoğu zaman normale yakındır. Bu, Merkezi Limit Teoremi’nin gücüdür.
14.5. Merkezi Limit Teoremi#
Normal dağılımlardan gelen değerleri toplarsak toplam normaldir. Çoğu dağılımda bu özellik yoktur; üstel dağılımdan gelen değerlerin toplamı üstel değildir.
Ancak birçok dağılımda n değer toplanınca n arttıkça toplamın dağılımı normale yakınsar. Değer dağılımının ortalaması m, varyansı s2 ise toplam, ortalaması n*m, varyansı n*s2 olan normal dağılıma yakınsar.
Bu sonuç Merkezi Limit Teoremi’dir (MLT). İstatistiksel analizin en yararlı araçlarındandır, ancak koşulları vardır:
Değerler aynı dağılımdan gelmelidir (bu koşul gevşetilebilir).
Bağımsız çekilmelidir. Korelasyonluysa MLT uygulanmaz (korelasyon fazla güçlü değilse yine çalışabilir).
Sonlu ortalama ve varyanslı dağılımdan gelmelidir; bazı uzun kuyruklu dağılımlarda uygulanmaz.
MLT, normal dağılımların doğadaki yaygınlığını açıklar. Canlıların birçok özelliği etkileri toplanan genetik ve çevresel etmenlerden etkilenir. Ölçtüğümüz özellikler çok sayıda küçük etkinin toplamı olduğundan dağılımları normal olma eğilimindedir.
MLT’nin nasıl çalıştığını ve ne zaman çalışmadığını üstel dağılımla başlayan deneylerle görelim. Aşağıdaki döngü üstel örneklemler üretip toplar ve her n büyüklüğünü 1001 toplamlık listeye eşleyen sözlük oluşturur.
lam = 1
df_sample_expo = pd.DataFrame()
for n in [1, 10, 100]:
df_sample_expo[n] = [np.sum(np.random.exponential(lam, n)) for _ in range(1001)]
Her toplam listesinin ortalamaları şöyledir.
df_sample_expo.mean()
1 0.989885
10 9.825744
100 100.022555
dtype: float64
Dağılımın ortalaması 1’dir; 10 değer toplarsak toplam ortalaması 10’a, 100 değer toplarsak 100’e yakındır.
Bu işlev oluşturduğumuz DataFrame’i alıp her toplam listesi için normal olasılık grafiği oluşturur.
def normal_plot_samples(df_sample, ylabel=""):
"""Toplam örneklemlerinin normal olasılık grafiklerini çizer."""
plt.figure(figsize=(6.8, 2.6))
for i, n in enumerate(df_sample):
plt.subplot(1, 3, i + 1)
normal_probability_plot(df_sample[n])
decorate(
title="n=%d" % n,
xticks=[],
yticks=[],
xlabel="Standart normal",
ylabel=ylabel,
)
Aşağıdaki şekil üç toplam listesinin normal olasılık grafiklerini gösterir (normal_plot_samples tanımı notebook’tadır).
normal_plot_samples(df_sample_expo, ylabel="Üstel değerlerin toplamı")
n=1 iken toplam üstel dağılımlıdır ve grafik doğru değildir. n=10 iken yaklaşık normal, n=100 iken normalden neredeyse ayırt edilemezdir.
Üstel dağılımdan daha az çarpık dağılımlarda toplam daha küçük n değerlerinde normale yakınsar; daha çarpıklarda daha uzun sürer. Örnek olarak lognormal değerlerin toplamına bakalım.
mu, sigma = 3.0, 1.0
df_sample_lognormal = pd.DataFrame()
for n in [1, 10, 100]:
df_sample_lognormal[n] = [
np.sum(np.random.lognormal(mu, sigma, n)) for _ in range(1001)
]
Aynı örneklem büyüklüklerinin normal olasılık grafikleri şöyledir.
normal_plot_samples(df_sample_lognormal, ylabel="Lognormal değerlerin toplamı")
n=1 iken normal model uymaz; n=10 ile çok iyileşmez. n=100 iken bile kuyruklar modelden açıkça sapar.
Lognormal dağılımın ortalama ve varyansı sonludur; toplam sonunda normale yakınsar. Ancak çok çarpık dağılımlarda hiçbir uygulanabilir örneklem büyüklüğünde yakınsamayabilir; bazen hiç yakınsamaz.
14.6. Merkezi Limit Teoremi’nin sınırları#
Pareto dağılımları lognormalden daha çarpıktır. Parametrelere bağlı olarak bazılarının sonlu ortalama ve varyansı yoktur; MLT uygulanmaz.
Göstermek için sonsuz ortalama ve varyansa sahip alpha=1 parametreli Pareto dağılımından değerler üreteceğiz.
alpha = 1.0
df_sample = pd.DataFrame()
for n in [1, 10, 100]:
df_sample[n] = [np.sum(np.random.pareto(alpha, n)) for _ in range(1001)]
Çeşitli örneklem büyüklüklerinin normal olasılık grafikleri şöyledir.
normal_plot_samples(df_sample, ylabel="Pareto değerlerinin toplamı")
n=100 iken bile toplam dağılımı normale benzemez.
Değerler korelasyonluysa MLT’nin uygulanmadığını da belirtmiştim. Bunu sınamak için, seri korelasyonu —örneklemdeki ardışık öğelerin korelasyonu— verilen rho değerine eşit üstel değerler üreten generate_expo_correlated işlevini kullanacağız. İşlev notebook’ta tanımlıdır.
def generate_normal_correlated(n, rho):
"""Standart normal dağılımdan korelasyonlu değerler dizisi üretir."""
xs = np.empty(n)
xs[0] = np.random.normal(0, 1)
sigma = np.sqrt(1 - rho**2)
for i in range(1, n):
xs[i] = rho * xs[i - 1] + np.random.normal(0, sigma)
return xs
Normal dağılımdan korelasyonlu bir dizi verildiğinde aşağıdaki işlev üstel dağılımdan korelasyonlu dizi üretir.
from scipy.stats import expon
def generate_expo_correlated(n, rho):
"""Üstel dağılımdan korelasyonlu değerler dizisi üretir."""
normal = generate_normal_correlated(n, rho)
uniform = norm.cdf(normal)
expo = expon.ppf(uniform)
return expo
Korelasyonlu normal değerleri normal CDF ile 0–1 arasında tekdüze değerlere, ardından üstel ters CDF ile üstel değerlere dönüştürür.
Aşağıdaki döngü her örneklem büyüklüğü için bir sütun ve her sütunda 1001 toplam bulunan DataFrame oluşturur.
rho = 0.8
df_sample = pd.DataFrame()
for n in [1, 10, 100]:
df_sample[n] = [np.sum(generate_expo_correlated(n, rho)) for _ in range(1001)]
Toplam dağılımlarının normal olasılık grafikleri şöyledir.
normal_plot_samples(df_sample, ylabel="Korelasyonlu değerlerin toplamı")
rho=0.8 iken ardışık öğeler güçlü korelasyonludur ve toplam yavaş yakınsar. Uzak öğeler arasında da güçlü korelasyon varsa hiç yakınsamayabilir.
Önceki bölüm MLT’nin çalışmasını, bu bölüm çalışmadığında ne olduğunu gösterdi. Şimdi nasıl kullanacağımıza bakalım.
14.7. MLT’nin uygulanması#
MLT’nin yararını görmek için Bölüm 9’daki ilk bebeklerle diğerlerinin ortalama gebelik süresi farkı sınamasına dönelim. NSFG verilerini yeniden kullanacağız; indirme yönergeleri notebook’tadır.
Aşağıdaki hücre verileri indirir.
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/nb/nsfg.py")
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/2002FemPreg.dct")
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/2002FemPreg.dat.gz")
Verileri okumak ve ilk bebeklerle diğerlerini ayırmak için get_nsfg_groups kullanacağız.
from nsfg import get_nsfg_groups
live, firsts, others = get_nsfg_groups()
Gördüğümüz gibi ilk bebekler ortalama biraz daha geç doğar; gözlenen fark yaklaşık 0,078 haftadır.
delta = firsts["prglngth"].mean() - others["prglngth"].mean()
delta
np.float64(0.07803726677754952)
Farkın şans eseri olup olmadığını görmek için sıfır hipotezi olarak gebelik sürelerinin ortalama ve varyansının iki grupta aynı olduğunu varsayacak, tüm canlı doğumlardan tahmin edeceğiz.
all_lengths = live["prglngth"]
m, s2 = all_lengths.mean(), all_lengths.var()
Gebelik süresi normal dağılım izlemez; yine de ortalamanın örnekleme dağılımına yaklaşmak için normal dağılım kullanabiliriz.
Aşağıdaki işlev değer dizisini alıp verilerle aynı ortalama ve varyanslı normal dağılımdan çekilen n büyüklüğündeki örneklemin ortalama dağılımını temsil eden Normal nesnesi döndürür.
def sampling_dist_mean(data, n):
mean, var = data.mean(), data.var()
dist = Normal(mean, var)
return dist.sum(n) / n
Sıfır hipotezi altında ilk doğumların ortalama gebelik süresinin örnekleme dağılımına normal yaklaşım şöyledir.
n1 = firsts["totalwgt_lb"].count()
dist_firsts = sampling_dist_mean(all_lengths, n1)
n1
np.int64(4363)
Diğer bebeklerin örnekleme dağılımı şöyledir.
n2 = others["totalwgt_lb"].count()
dist_others = sampling_dist_mean(all_lengths, n2)
n2
np.int64(4675)
Farkın örnekleme dağılımını şöyle hesaplayabiliriz.
dist_diff = dist_firsts - dist_others
dist_diff
Normal(0.0, 0.003235837567930557)
Ortalama 0’dır; aynı dağılımdan iki örneklem çekersek ortalama farkını ortalama olarak 0 bekleriz. Örnekleme dağılımının 0,0032 varyansı, şans nedeniyle ne kadar değişim beklediğimizi gösterir.
Bu yaklaşımı doğrulamak için dağılımı yeniden örneklemeyle de tahmin edebiliriz.
sample_firsts = [np.random.choice(all_lengths, n1).mean() for i in range(1001)]
sample_others = [np.random.choice(all_lengths, n2).mean() for i in range(1001)]
sample_diffs = np.subtract(sample_firsts, sample_others)
Yeniden örneklenmiş farkların ampirik CDF’siyle normal model şöyledir. Dikey noktalı çizgiler gözlenen pozitif ve negatif farkı gösterir.
dist_diff.plot_cdf(**model_options)
Cdf.from_seq(sample_diffs).plot(label="örneklem")
plt.axvline(delta, ls=":")
plt.axvline(-delta, ls=":")
decorate(xlabel="Gebelik süresi farkı", ylabel="CDF")
Örneklem büyüklükleri fazla ve ölçümlerin çarpıklığı sınırlı olduğundan örnekleme dağılımı normalle iyi yaklaşıklandırılır. Bu nedenle p-değerini normal CDF ile hesaplayabiliriz. Aşağıdaki yöntem normal CDF’yi hesaplar.
%%add_method_to Normal
def cdf(self, xs):
sigma = np.sqrt(self.sigma2)
return norm.cdf(xs, self.mu, sigma)
Sıfır hipotezi altında delta kadar büyük fark olasılığı, örnekleme dağılımının sağ kuyruk alanıdır.
right = 1 - dist_diff.cdf(delta)
right
np.float64(0.08505405315526993)
-delta kadar negatif fark olasılığı ise sol kuyruk alanıdır.
left = dist_diff.cdf(-delta)
left
np.float64(0.08505405315526987)
Normal dağılım simetrik olduğu için left ile right aynıdır. Toplamları pozitif ya da negatif, delta kadar büyük fark olasılığıdır.
left + right
np.float64(0.1701081063105398)
p-değeri 0,170’tir; Bölüm 9’da yeniden örneklemeyle hesapladığımız tahminle tutarlıdır.
Bu hesap bağımsız örneklemler \(t\)-testine benzer. SciPy’nin ttest_ind işlevi iki örneklemin ortalama farkı p-değerini hesaplar.
from scipy.stats import ttest_ind
result = ttest_ind(firsts["prglngth"], others["prglngth"])
result.pvalue
np.float64(0.16755412639415002)
Örneklemler büyükken \(t\)-testi sonucu normal dağılımlarla hesapladığımıza yakındır. Adı, normal yerine \(t\) dağılımına dayanmasından gelir. Sonraki bölümde \(t\) dağılımını korelasyonun anlamlılığını sınamakta da kullanacağız.
14.8. Korelasyon testi#
Bölüm 9’da doğum ağırlığı ile annenin yaşı arasındaki korelasyonu permütasyonla sınayıp p-değerinin 0,001’den küçük olduğunu bulduk.
Şimdi analitik olarak yapabiliriz. Yöntem şu sonuca dayanır: Normal dağılımlardan n büyüklüğünde iki örneklem üretip Pearson korelasyonu r’yi hesaplar ve aşağıdaki işlevle dönüştürürsek:
def transform_correlation(r, n):
return r * np.sqrt((n - 2) / (1 - r**2))
Dönüştürülmüş korelasyonlar n-2 parametreli \(t\) dağılımı izler. Bunu görmek için standart normal dağılımdan korelasyonsuz örneklemler üreten aşağıdaki işlevi kullanacağız.
def generate_data(n):
"""Standart normal dağılımdan korelasyonsuz diziler üretir."""
xs = np.random.normal(0, 1, n)
ys = np.random.normal(0, 1, n)
return xs, ys
Korelasyonlarını hesaplamak için bu işlevi kullanacağız.
def correlation(data):
xs, ys = data
return np.corrcoef(xs, ys)[0, 1]
Aşağıdaki döngü çok sayıda örneklem çifti üretip korelasyonlarını bir listede toplar.
n = 100
rs = [correlation(generate_data(n)) for i in range(1001)]
Ardından dönüştürülmüş korelasyonları hesaplayacağız.
ts = transform_correlation(np.array(rs), n)
Bu ts değerlerinin \(t\) dağılımı izleyip izlemediğini denetlemek için \(t\) dağılımının CDF’sini temsil eden nesne oluşturan işlevi kullanacağız.
from scipy.stats import t as student_t
def make_student_cdf(df):
"""Student t dağılımının CDF’sini hesaplar."""
ts = np.linspace(-3, 3, 101)
ps = student_t.cdf(ts, df=df)
return Cdf(ps, ts)
\(t\) dağılımının df parametresi “degrees of freedom”, yani serbestlik derecesidir. Şekil n-2 parametreli \(t\) dağılımının CDF’sini dönüştürülmüş korelasyonların ampirik CDF’siyle gösterir.
make_student_cdf(df=n - 2).plot(**model_options)
cdf_ts = Cdf.from_seq(ts)
cdf_ts.plot(label="rastgele normal değerler")
decorate(xlabel="Dönüştürülmüş korelasyon", ylabel="CDF")
Korelasyonsuz normal örneklemlerin dönüştürülmüş korelasyonlarının \(t\) dağılımı izlediğini gösterir.
Başka dağılımlardan örneklemlerde tam olarak \(t\) dağılımı izlemez, ancak örneklem büyüdükçe yakınsar. Anne yaşı ile doğum ağırlığı korelasyonuna uygulanıp uygulanmadığına bakalım. Canlı doğumlardan geçerli verili satırları seçeceğiz.
valid = live.dropna(subset=["agepreg", "totalwgt_lb"])
n = len(valid)
n
9038
Gerçek korelasyon yaklaşık 0,07’dir.
data = valid["agepreg"].values, valid["totalwgt_lb"].values
r_actual = correlation(data)
r_actual
np.float64(0.06883397035410906)
Bölüm 9’daki gibi örneklemleri permüte ederek sıfır hipotezini simüle edebiliriz.
def permute(data):
"""x değerlerini karıştırır."""
xs, ys = data
new_xs = xs.copy()
np.random.shuffle(new_xs)
return new_xs, ys
Çok sayıda permütasyon üretip korelasyonlarını hesaplarsak sıfır hipotezi altındaki korelasyon dağılımından örneklem elde ederiz.
permuted_corrs = [correlation(permute(data)) for i in range(1001)]
Dönüştürülmüş korelasyonları şöyle hesaplarız.
ts = transform_correlation(np.array(permuted_corrs), n)
Şekil ts değerlerinin ampirik CDF’sini n-2 parametreli \(t\) dağılımı CDF’siyle gösterir.
make_student_cdf(n - 2).plot(**model_options)
Cdf.from_seq(ts).plot(label="permüte edilmiş veri")
decorate(xlabel="Dönüştürülmüş korelasyon", ylabel="CDF")
Model ampirik dağılıma iyi uyar; gözlenen korelasyonun p-değerini hesaplamakta kullanabiliriz. Önce gözlenen korelasyonu dönüştüreceğiz.
t_actual = transform_correlation(r_actual, n)
\(t\) CDF’siyle sıfır hipotezi altında t_actual kadar büyük değer olasılığını hesaplayabiliriz.
right = 1 - student_t.cdf(t_actual, df=n - 2)
right
np.float64(2.8614777214386322e-11)
-t_actual kadar negatif değer olasılığını da hesaplayabiliriz.
left = student_t.cdf(-t_actual, df=n - 2)
left
np.float64(2.8614735536574775e-11)
Toplamları, pozitif ya da negatif r_actual kadar büyük korelasyon olasılığıdır.
left + right
np.float64(5.72295127509611e-11)
SciPy aynı hesabı yapıp gözlenen korelasyonun p-değerini döndüren işlev sağlar.
from scipy.stats import pearsonr
corr, p_value = pearsonr(*data)
p_value
np.float64(5.722947107314301e-11)
Sonuçlar neredeyse aynıdır.
Yeniden örneklemeye dayanarak p-değerinin 0,001’den küçük olduğunu, ancak çok fazla yeniden örnekleme yapmadan ne kadar küçük olduğunu söyleyememiştik. Analitik yöntemlerle küçük p-değerlerini hızla hesaplarız.
Uygulamada bunun önemi olmayabilir. Genellikle p-değeri 0,001’den küçükse gözlenen etkinin şanstan kaynaklanmasının düşük olasılıklı olduğu sonucuna varabiliriz; ne kadar düşük olduğunu tam bilmek çoğu zaman önemli değildir.
14.9. Ki-kare testi#
Bölüm 9’da gözlenen sonuçlara dayanarak bir zarın hileli olup olmadığını sınadık.
from empiricaldist import Hist
qs = np.arange(1, 7)
freqs = [8, 9, 19, 5, 8, 11]
observed = Hist(freqs, qs)
observed.index.name = "outcome"
observed
| freqs | |
|---|---|
| outcome | |
| 1 | 8 |
| 2 | 9 |
| 3 | 19 |
| 4 | 5 |
| 5 | 8 |
| 6 | 11 |
Önce her sonucun beklenen sıklığını hesapladık.
num_rolls = observed.sum()
outcomes = observed.qs
expected = Hist(num_rolls / 6, outcomes)
Ardından ki-kare istatistiğini aşağıdaki işlevle hesapladık.
def chi_squared_stat(observed, expected):
diffs = (observed - expected) ** 2
ratios = diffs / expected
return np.sum(ratios.values.flatten())
observed_chi2 = chi_squared_stat(observed, expected)
Ki-kare istatistiği bu tür verilerde yaygındır; çünkü sıfır hipotezi altındaki örnekleme dağılımı verimli hesaplanabilen ki-kare dağılımına yakınsar. Görmek için adil zar atışını simüle eden işlevi kullanacağız.
def simulate_dice(observed):
n = np.sum(observed)
rolls = np.random.choice(observed.qs, n, replace=True)
hist = Hist.from_seq(rolls)
return hist
Aşağıdaki döngü simülasyonu çok kez çalıştırıp sonuçların ki-kare istatistiğini hesaplar.
simulated_chi_squared = [
chi_squared_stat(simulate_dice(observed), expected) for i in range(1001)
]
cdf_simulated = Cdf.from_seq(simulated_chi_squared)
Sonuçların ki-kare dağılımı izleyip izlemediğini denetlemek için df parametreli ki-kare dağılımının CDF’sini hesaplayan işlevi kullanacağız.
from scipy.stats import chi2 as chi2_dist
def chi_squared_cdf(df):
"""Ki-kare CDF’sine kesikli yaklaşım."""
xs = np.linspace(0, 21, 101)
ps = chi2_dist.cdf(xs, df=df)
return Cdf(ps, xs)
n olası sonuçta simüle edilmiş istatistikler n-1 parametreli ki-kare dağılımı izlemelidir.
n = len(observed)
cdf_model = chi_squared_cdf(df=n - 1)
Simüle edilmiş istatistiklerin ampirik CDF’siyle ki-kare CDF şöyledir.
cdf_model.plot(**model_options)
cdf_simulated.plot(label="simülasyon")
decorate(xlabel="Ki-kare istatistiği", ylabel="CDF")
Model simülasyon sonuçlarına iyi uyar; sıfır hipotezi altında observed_chi2 kadar büyük değer olasılığını hesaplayabiliriz.
p_value = 1 - chi2_dist.cdf(observed_chi2, df=n - 1)
p_value
np.float64(0.04069938850404997)
SciPy aynı hesabı yapan bir işlev sağlar.
from scipy.stats import chisquare
chi2_stat, p_value = chisquare(f_obs=observed, f_exp=expected)
Sonuç hesapladığımız p-değeriyle aynıdır.
p_value
np.float64(0.04069938850405)
Ki-kare istatistiğinin üstünlüğü sıfır hipotezi altındaki dağılımının verimli hesaplanmasıdır. Ancak bağlama göre gözlenen ve beklenen sonuç farkını en iyi nicelleştiren istatistik olmayabilir.
14.10. Hesaplama ve analiz#
Bu kitap yeniden örnekleme ve permütasyon gibi hesaplamalı yöntemlere odaklanır. Analitik yöntemlere göre üstünlükleri vardır:
Açıklaması ve anlaşılması kolaydır. Giriş istatistiğinin en zor konularından hipotez testinde birçok öğrenci p-değerini tam anlamaz. Bölüm 9’daki sıfır hipotezini simüle etme ve test istatistiği hesaplama yaklaşımı temel fikri daha açık kılar.
Dirençli ve esnektir. Analitik yöntemler çoğu zaman uygulamada geçersiz varsayımlara dayanır. Hesaplamalı yöntemler daha az varsayım gerektirir, daha kolay uyarlanır ve genişletilir.
Hata ayıklanabilir. Analitik yöntemler sayıları verip sonuç aldığınız kara kutular gibi olabilir. İnce hatalar yapmak kolay, sonuçlara güvenmek ve sorunu tanılamak zordur. Hesaplamalı yöntemler aşamalı geliştirme ve sınamaya uygundur; sonuçlara güveni artırır.
Dezavantajları da vardır:
Hesaplamalı yöntemler yavaş olabilir.
Yeniden örnekleme gibi rastgele yöntemler her seferinde aynı sonucu üretmez; doğruluğu denetlemek zorlaşır.
Artı ve eksileri dikkate alarak şu süreci öneriyorum:
Keşif sırasında hesaplamalı yöntem kullanın. Yeterli yanıt bulur ve çalışma süresi uygunsa durabilirsiniz.
Süre uygun değilse eniyileştirme olanaklarını arayın. Analitik yöntem kullanmak seçeneklerden biridir.
Hesaplamalı yöntemi analitik olanla değiştirmek uygunsa iki sonucu karşılıklı doğrulamak için hesaplamalı yöntemi karşılaştırma temeli olarak kullanın.
Birçok uygulamalı problemde hesaplamalı yöntemlerin süresi sorun değildir ve ilk adımın ötesine geçmek gerekmez.
14.11. Sözlük#
normal olasılık grafiği: Gözlenen değerleri normal dağılımın nicelikleriyle karşılaştırıp verinin normale ne kadar yaklaştığını gösteren grafik.
bağımsız örneklemler \(t\)-testi: İki bağımsız grubun ortalamaları arasındaki gözlenen farkın p-değerini hesaplama yöntemi.
\(t\) dağılımı: Farkın 0 olduğu sıfır hipotezi altında ortalama farkının örnekleme dağılımını ve dönüştürülmüş korelasyonların dağılımını modelleyen dağılım.
ki-kare dağılımı: Ki-kare istatistiğinin örnekleme dağılımını modelleyen dağılım.
ki-kare istatistiği: İki kesikli dağılım arasındaki farkın büyüklüğünü nicelleştiren test istatistiği.
14.12. Alıştırmalar#
14.12.1. Alıştırma 14.1#
Bu bölümde erkek ve dişi penguenlerin ağırlıklarını karşılaştırıp farkın güven aralığını hesapladık. Şimdi yüzgeç uzunluğu için aynısını yapalım. Gözlenen fark yaklaşık 4,6 mm’dir.
grouped = adelie.groupby("Sex")
lengths_male = grouped.get_group("MALE")["Flipper Length (mm)"]
lengths_female = grouped.get_group("FEMALE")["Flipper Length (mm)"]
observed_diff = lengths_male.mean() - lengths_female.mean()
observed_diff
np.float64(4.616438356164366)
Grupların aynı büyüklükte olmadığına dikkat ederek iki grubun ortalama yüzgeç uzunluğunun örnekleme dağılımlarını temsil eden Normal nesnelerini sampling_dist_mean ile oluşturun. Farkın örnekleme dağılımını ve %90 güven aralığını hesaplayın.
14.12.2. Alıştırma 14.2#
NSFG verilerinde bebeğin doğum ağırlığı ile annenin yaşı arasındaki korelasyonu hesaplayıp p-değerini \(t\) dağılımıyla bulduk. Şimdi babanın yaşını içeren hpagelb sütunuyla aynısını yapalım.
valid = live.dropna(subset=["hpagelb", "totalwgt_lb"])
n = len(valid)
n
8933
Gözlenen korelasyon yaklaşık 0,065’tir.
data = valid["hpagelb"].values, valid["totalwgt_lb"].values
r_actual = correlation(data)
r_actual
np.float64(0.06468629895432172)
Dönüştürülmüş korelasyon t_actual değerini hesaplayın. p-değeri için \(t\) CDF’sini kullanın; korelasyon istatistiksel olarak anlamlı mı? Sonucu SciPy’nin pearsonr işleviyle denetleyin.
14.12.3. Alıştırma 14.3#
Bölüm 9’daki bir alıştırmada birçok memelide cinsiyet oranının annenin yaşı, bedeni, sağlığı ve sosyal statüsü gibi “anne koşullarına” bağlı olduğunu öne süren Trivers–Willard hipotezini ele aldık. İnsan çalışmalarının sonuçları karışıktır.
Ki-kare testi alıştırması olarak bebeğin cinsiyeti ile annenin medeni durumu arasında ilişki olup olmadığına bakalım. Başlangıç yönergeleri notebook’tadır.
Önce erkek ve kız bebeklerin annelerini ayıracağız.
male = live.query("babysex == 1")
female = live.query("babysex == 2")
Şimdi her grup için bir sütunu ve medeni durumu kodlayan her fmarital değeri için bir satırı olan DataFrame oluşturacağız:
1 evli
2 dul
3 boşanmış
4 ayrı
5 hiç evlenmemiş
observed = pd.DataFrame()
observed["male"] = male["fmarital"].value_counts().sort_index()
observed["female"] = female["fmarital"].value_counts().sort_index()
observed
| male | female | |
|---|---|---|
| fmarital | ||
| 1 | 2576 | 2559 |
| 2 | 56 | 54 |
| 3 | 568 | 572 |
| 4 | 355 | 330 |
| 5 | 1086 | 985 |
Sıfır hipotezi medeni durum dağılımının iki grupta aynı olmasıdır; tüm veri kümesinden hesaplayabiliriz.
from empiricaldist import Pmf
pmf_fmarital = Pmf.from_seq(live["fmarital"])
pmf_fmarital
| probs | |
|---|---|
| fmarital | |
| 1 | 0.561653 |
| 2 | 0.012024 |
| 3 | 0.124617 |
| 4 | 0.075208 |
| 5 | 0.226498 |
Beklenen değerler için pmf_marital olasılıklarını her sütundaki toplam olgu sayısıyla çarparız.
expected = pd.DataFrame()
expected["male"] = pmf_fmarital * observed["male"].sum()
expected["female"] = pmf_fmarital * observed["female"].sum()
expected
| male | female | |
|---|---|---|
| fmarital | ||
| 1 | 2606.630739 | 2527.437691 |
| 2 | 55.805641 | 54.110188 |
| 3 | 578.349366 | 560.778312 |
| 4 | 349.038916 | 338.434631 |
| 5 | 1051.175339 | 1019.239178 |
observed ve expected ile ki-kare istatistiğini hesaplayın. Ki-kare CDF’siyle p-değerini bulun. Serbestlik derecesi, gözlenen DataFrame değer sayısı n için n-1 olmalıdır. Ardından SciPy’nin chisquare işleviyle istatistik ve p-değerini hesaplayın. İpucu: Her sütuna ayrı sınama yerine tüm DataFrame’i tek test kabul etmek için axis=None kullanın.
Bu test Trivers–Willard hipotezini destekliyor mu?
14.12.4. Alıştırma 14.4#
Bu bölümde gruplar arası farkları analiz etmekte kullandığımız yöntem, yaygın bir deney tasarımı olan “farkların farkını” incelemeye genişletilebilir. 2014 tarihli bir çalışmanın, öğrenci mühendislik ekiplerinde cinsiyet kalıp yargılarına dayalı görev dağılımını azaltmayı amaçlayan müdahalenin etkilerini araştıran verilerini kullanacağız.
Stein, L. A., Aragon, D., Moreno, D. ve Goodman, J. (Ekim 2014). Evidence for the persistent effects of an intervention to mitigate gender-stereotypical task allocation within student engineering teams. 2014 IEEE Frontiers in Education Conference (FIE) Proceedings (s. 1–9). IEEE.
http://ieeexplore.ieee.org/document/7044435/ adresinden edinilebilir.
Müdahale öncesi ve sonrasında öğrenciler, sınıf projelerinin her yönüne katkılarını 7 puanlık ölçekte değerlendirdi.
Müdahale öncesinde erkek öğrencilerin programlama puanı kadınlardan yüksekti: erkek ortalaması 3,57 ve standart hatası 0,28; kadın ortalaması 1,91 ve standart hatası 0,32 idi.
Müdahale sonrasında cinsiyet farkı küçüldü: erkek ortalaması 3,44 (SH 0,16), kadın ortalaması 3,18 (SH 0,16) idi.
Erkek ve kadın öğrenciler için müdahale öncesi ve sonrası tahmini ortalamaların örnekleme dağılımlarını temsil eden dört
Normalnesnesi oluşturun. Tahmini ortalamaların standart hataları verildiğinden örneklem büyüklüğü gerekmez.Müdahale öncesi ve sonrası cinsiyet farkının, yani ortalama farklarının örnekleme dağılımlarını hesaplayın.
Farkların farkının, yani aralık büyüklüğündeki değişimin örnekleme dağılımını hesaplayın. %95 güven aralığı ve p-değeri bulun.
Müdahaleden sonra cinsiyet farkının azaldığına dair kanıt var mı?
Think Stats: Python ile Keşifsel Veri Analizi, 3. Baskı
Telif hakkı 2024 Allen B. Downey
Kod lisansı: MIT Lisansı
Metin lisansı: Creative Commons Atıf-Gayriticari-AynıLisanslaPaylaş 4.0 Uluslararası