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.
5. Dağılımların modellenmesi#
Şimdiye kadar kullandığımız dağılımlara ampirik dağılımlar denir; çünkü bunlar ampirik gözlemlere, başka bir deyişle verilere dayanır. Gerçek dünyada karşılaştığımız birçok veri kümesine, çoğunlukla basit bir matematiksel fonksiyona dayanan kuramsal bir dağılımla yakından yaklaşılabilir. Bu bölümde bu kuramsal dağılımlardan bazılarını ve bunlarla modellenebilecek veri kümelerini ele alacağız.
Örnekler üzerinden şunları göreceğiz:
Bir skeet atıcılığı yarışmasındaki isabet ve ıskalama sayıları binom dağılımıyla iyi modellenir.
Buz hokeyi ve futbol gibi oyunlarda bir maçtaki gol sayısı Poisson dağılımını, goller arasındaki süre ise üstel dağılımı izler.
Doğum ağırlıkları Gauss dağılımı olarak da bilinen normal dağılımı, yetişkin ağırlıkları ise lognormal dağılımı izler.
Bu dağılımlara ya da bu sporlara aşina değilseniz bilmeniz gerekenleri açıklayacağım. Her örnekte basit bir modele dayalı simülasyonla başlayıp simülasyon sonuçlarının kuramsal bir dağılımı izlediğini göstereceğiz. Ardından gerçek verilerin modele ne kadar iyi uyduğunu inceleyeceğiz.
Bu notebook’u Colab’da çalıştırmak için buraya tıklayın.
5.1. Binom dağılımı#
İlk örnek olarak, yarışmacıların havaya fırlatılan kil diskleri av tüfeğiyle vurmaya çalıştığı skeet atıcılığını ele alacağız. Olimpiyatlar dâhil uluslararası yarışmalarda, her biri 25 hedef içeren beş tur ve kazananı belirlemek için gerektiğinde ek turlar yapılır.
Bir skeet atıcılığı yarışması modeli olarak her katılımcının her hedefi vurma olasılığının aynı, p, olduğunu varsayalım. Elbette bu model bir basitleştirmedir: Gerçekte bazı yarışmacıların başarı olasılığı diğerlerinden yüksektir ve tek bir yarışmacının olasılığı bile bir denemeden diğerine değişebilir. Yine de gerçekçi olmasa bile bu model, göreceğimiz üzere şaşırtıcı ölçüde doğru tahminler üretir.
Modeli simüle etmek için aşağıdaki fonksiyonu kullanacağım. Fonksiyon hedef sayısını n, her bir hedefi vurma olasılığını p olarak alır; isabetleri ve ıskalamaları göstermek üzere 1 ve 0’lardan oluşan bir dizi döndürür.
def flip(n, p):
choices = [1, 0]
probs = [p, 1 - p]
return np.random.choice(choices, n, p=probs)
Aşağıdaki örnek, her hedefi vurma olasılığının %90 olduğu 25 hedeflik bir turu simüle eder.
# Her çalıştırmada aynı sonuçları almak için rastgele sayı üretecini başlat
np.random.seed(1)
flip(25, 0.9)
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
1, 1, 1])
Daha uzun bir dizi üretip sonuçların Pmf’sini hesaplarsak 1 ve 0 oranlarının, en azından yaklaşık olarak, doğru olduğunu doğrulayabiliriz.
from empiricaldist import Pmf
seq = flip(1000, 0.9)
pmf = Pmf.from_seq(seq)
pmf
| probs | |
|---|---|
| 0 | 0.101 |
| 1 | 0.899 |
Artık flip fonksiyonunu kullanarak bir skeet atıcılığı turunu simüle edebilir ve isabet sayısını döndürebiliriz.
def simulate_round(n, p):
seq = flip(n, p)
return seq.sum()
Büyük bir yarışmada 200 yarışmacının beşer tur atış yaptığını ve hepsinin hedefi vurma olasılığının aynı, p=0.9, olduğunu varsayalım. simulate_round fonksiyonunu 1000 kez çağırarak böyle bir yarışmayı simüle edebiliriz.
n = 25
p = 0.9
results_sim = [simulate_round(n, p) for i in range(1000)]
Ortalama puan, n ile p’nin çarpımı olan 22.5 değerine yakındır.
np.mean(results_sim), n * p
(np.float64(22.522), 22.5)
Sonuçların dağılımı aşağıdaki gibi görünür.
from empiricaldist import Pmf
pmf_sim = Pmf.from_seq(results_sim, name="simülasyon sonuçları")
pmf_sim.bar()
decorate(xlabel="İsabet sayısı", ylabel="PMF")
Tepe noktası ortalamaya yakındır ve dağılım sola çarpıktır.
Simülasyon çalıştırmak yerine bu dağılımı önceden tahmin edebilirdik. Matematiksel olarak bu sonuçların dağılımı, PMF’si kolayca hesaplanabilen bir binom dağılımını izler.
from scipy.special import comb
def binomial_pmf(k, n, p):
return comb(n, k) * (p**k) * ((1 - p) ** (n - k))
SciPy, n nesneden aynı anda k tanesini seçmenin kombinasyon sayısını hesaplayan comb fonksiyonunu sağlar; bu ifade genellikle “n’in k’lı kombinasyonu” diye okunur.
binomial_pmf, p verildiğinde n denemede k isabet elde etme olasılığını hesaplar. Bu fonksiyonu bir k değerleri aralığıyla çağırırsak sonuçların dağılımını temsil eden bir Pmf oluşturabiliriz.
ks = np.arange(16, n + 1)
ps = binomial_pmf(ks, n, p)
pmf_binom = Pmf(ps, ks, name="binom modeli")
Simülasyon sonuçlarıyla karşılaştırıldığında aşağıdaki gibi görünür.
from thinkstats import two_bar_plots
two_bar_plots(pmf_sim, pmf_binom)
decorate(xlabel="İsabet sayısı", ylabel="PMF")
İki dağılım benzerdir; aralarındaki küçük farklar simülasyon sonuçlarındaki rastgele değişimden kaynaklanır. Bu uyum şaşırtıcı değildir, çünkü simülasyon ile model aynı varsayımlara, özellikle her denemenin aynı başarı olasılığına sahip olduğu varsayımına dayanır. Bir model için daha güçlü bir sınama, modelin gerçek verilerle ne kadar iyi uyuştuğudur.
2020 Yaz Olimpiyatları erkekler skeet atıcılığı yarışmasının Wikipedia sayfasından eleme turlarının sonuçlarını gösteren bir tablo çıkarabiliriz. Verileri indirme talimatları bu bölümün notebook’unda yer alır.
15 Temmuz 2024 tarihinde https://en.wikipedia.org/wiki/Shooting_at_the_2020_Summer_Olympics_–_Men’s_skeet adresinden indirilmiştir.
filename = "Shooting_at_the_2020_Summer_Olympics_Mens_skeet"
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/" + filename)
tables = pd.read_html(filename)
table = tables[6]
table.head()
| Rank | Athlete | Country | 1 | 2 | 3 | 4 | 5 | Total[3] | Shoot-off | Notes | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | Éric Delaunay | France | 25 | 25 | 25 | 24 | 25 | 124 | +6 | Q, OR |
| 1 | 2 | Tammaro Cassandro | Italy | 24 | 25 | 25 | 25 | 25 | 124 | +5 | Q, OR |
| 2 | 3 | Eetu Kallioinen | Finland | 25 | 25 | 24 | 25 | 24 | 123 | NaN | Q |
| 3 | 4 | Vincent Hancock | United States | 25 | 25 | 25 | 25 | 22 | 122 | +8 | Q |
| 4 | 5 | Abdullah Al-Rashidi | Kuwait | 25 | 25 | 24 | 25 | 23 | 122 | +7 | Q |
Tabloda her yarışmacı için bir satır ve beş turun her biri için bir sütun vardır. Bu sonuçları içeren sütunları seçip NumPy’nin flatten fonksiyonuyla tek bir diziye dönüştüreceğiz.
columns = ["1", "2", "3", "4", "5"]
results = table[columns].values.flatten()
30 yarışmacı olduğundan, her biri 25 atıştan oluşan 150 turun sonuçlarına sahibiz. Toplam 3750 denemede 3575 isabet vardır.
total_shots = 25 * len(results)
total_hits = results.sum()
n, total_shots, total_hits
(25, 3750, np.int64(3575))
Buna göre genel başarı oranı %95,3’tür.
p = total_hits / total_shots
p
np.float64(0.9533333333333334)
Şimdi n=25 ve az önce hesapladığımız p değeriyle binom dağılımını temsil eden bir Pmf hesaplayalım.
ps = binomial_pmf(ks, n, p)
pmf_binom = Pmf(ps, ks, name="binom modeli")
Bunu gerçek sonuçların Pmf’siyle karşılaştırabiliriz.
pmf_results = Pmf.from_seq(results, name="gerçek sonuçlar")
two_bar_plots(pmf_results, pmf_binom)
decorate(xlabel="İsabet sayısı", ylabel="PMF")
Binom modeli veri dağılımına iyi uyar; üstelik bütün yarışmacıların aynı ve değişmeyen yeteneğe sahip olduğu yönündeki gerçekçi olmayan varsayıma rağmen.
5.2. Poisson dağılımı#
Spor karşılaşmalarının sonuçlarının öngörülebilir örüntüler izlediği bir başka örnek olarak buz hokeyi maçlarında atılan gol sayısına bakalım.
Önce 60 dakikalık, yani 3600 saniyelik bir maçı simüle edeceğiz. Takımların maç başına ortalama toplam 6 gol attığını ve herhangi bir saniyedeki gol olasılığının, p, aynı olduğunu varsayacağız.
n = 3600
m = 6
p = m / 3600
p
0.0016666666666666668
Aşağıdaki fonksiyonla n saniyeyi simüle edip atılan toplam gol sayısını döndürebiliriz.
def simulate_goals(n, p):
return flip(n, p).sum()
Çok sayıda maç simüle edersek maç başına ortalama gol sayısının 6’ya yakın olduğunu doğrulayabiliriz.
goals = [simulate_goals(n, p) for i in range(1001)]
np.mean(goals)
np.float64(6.021978021978022)
Bu sonuçları modellemek için binom dağılımını kullanabilirdik; ancak n büyük ve p küçük olduğunda sonuçlar Poisson dağılımıyla da iyi modellenir. Bu dağılım, genellikle “lambda” diye okunan Yunanca λ harfiyle gösterilen bir değerle belirlenir. Kodda bu değeri lam değişkeniyle temsil ediyoruz (lambda, Python’da anahtar sözcük olduğu için geçerli bir değişken adı değildir). Örneğimizde lam, maç başına 6 gol olan gol atma hızını temsil eder.
Poisson dağılımının PMF’sini hesaplamak kolaydır. lam verildiğinde, bir maçta k gol görme olasılığını aşağıdaki fonksiyonla hesaplayabiliriz.
from scipy.special import factorial
def poisson_pmf(k, lam):
"""Poisson PMF'yi hesaplar.
k (int veya dizi benzeri): Gerçekleşme sayısı
lam (float): Poisson dağılımının hız parametresi (λ)
Döndürür: float veya ndarray
"""
return (lam**k) * np.exp(-lam) / factorial(k)
SciPy, 1 ile k arasındaki tam sayıların çarpımını hesaplayan factorial fonksiyonunu sağlar.
poisson_pmf fonksiyonunu bir k değerleri aralığıyla çağırırsak sonuçların dağılımını temsil eden bir Pmf oluşturabiliriz.
lam = 6
ks = np.arange(20)
ps = poisson_pmf(ks, lam)
pmf_poisson = Pmf(ps, ks, name="Poisson modeli")
Dağılımın ortalamasının 6’ya yakın olduğunu da doğrulayabiliriz.
pmf_poisson.normalize()
pmf_poisson.mean()
np.float64(5.999925498375129)
Aşağıdaki şekil simülasyon sonuçlarını aynı ortalamaya sahip Poisson dağılımıyla karşılaştırır.
pmf_sim = Pmf.from_seq(goals, name="simülasyon")
two_bar_plots(pmf_sim, pmf_poisson)
decorate(xlabel="Gol sayısı", ylabel="PMF")
Dağılımlar, rastgele değişimden kaynaklanan küçük farklar dışında benzerdir. Bu şaşırtıcı değildir; çünkü simülasyon ile Poisson modeli, maçın herhangi bir saniyesinde gol atma olasılığının aynı olduğu varsayımına dayanır. Daha güçlü bir sınama, modelin gerçek verilere ne kadar iyi uyduğuna bakmaktır.
HockeyReference sitesinden Ulusal Hokey Ligi’nin (NHL) 2023-2024 normal sezonundaki bütün maçların sonuçlarını indirdim; play-off maçları buna dâhil değildir. Uzatma ve eşitliği bozan penaltı atışları hariç, 60 dakikalık normal oyun süresinde atılan gollere ilişkin bilgileri çıkardım. Sonuçlar, her maç için bir anahtar ve maçın başlangıcından golün atıldığı ana kadar geçen saniyelerin listesini içeren bir HDF dosyasındadır. Verileri indirme talimatları bu bölümün notebook’unda yer alır.
Ham veriler 16 Temmuz 2024 tarihinde https://www.hockey-reference.com/leagues/NHL_2024_games.html adresinden indirilmiştir.
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/nhl_2023_2024.hdf")
Dosyadaki anahtarları aşağıdaki gibi okuruz.
filename = "nhl_2023_2024.hdf"
with pd.HDFStore(filename, "r") as store:
keys = store.keys()
len(keys), keys[0]
(1312, '/202310100PIT')
Normal sezonda 1312 maç oynanmıştır. Her anahtar maçın tarihini ve ev sahibi takımın üç harfli kısaltmasını içerir. Bir anahtarı arayıp gollerin atıldığı zamanların listesini almak için read_hdf fonksiyonunu kullanabiliriz.
times = pd.read_hdf(filename, key=keys[0])
times
0 424
1 1916
2 2137
3 3005
4 3329
5 3513
dtype: int64
Sezonun ilk maçında altı gol atılmıştır: İlki oyunun 424. saniyesinde, sonuncusu ise maçın bitimine yalnızca 87 saniye kala, 3513. saniyede.
3600 - times[5]
np.int64(87)
Aşağıdaki döngü bütün maçların sonuçlarını okur, her maçtaki gol sayısını hesaplar ve sonuçları bir listede saklar.
goals = []
for key in keys:
times = pd.read_hdf(filename, key=key)
n = len(times)
goals.append(n)
Maç başına ortalama gol sayısı 6’nın biraz üzerindedir.
lam = np.mean(goals)
lam
np.float64(6.0182926829268295)
Verilerle aynı ortalamaya sahip Poisson dağılımını temsil eden bir Pmf oluşturmak için poisson_pmf fonksiyonunu kullanabiliriz.
ps = poisson_pmf(ks, lam)
pmf_poisson = Pmf(ps, ks, name="Poisson modeli")
Verilerin PMF’siyle karşılaştırıldığında aşağıdaki gibi görünür.
pmf_goals = Pmf.from_seq(goals, name="atılan goller")
two_bar_plots(pmf_goals, pmf_poisson)
decorate(xlabel="Gol sayısı", ylabel="PMF")
Poisson dağılımı verilere iyi uyar; bu da hokeydeki gol atma süreci için iyi bir model olduğunu düşündürür.
5.3. Üstel dağılım#
Önceki bölümde, maçın herhangi bir saniyesinde gol atılma olasılığının aynı olduğu basit bir hokey maçı modelini simüle ettik. Aynı model altında ilk gole kadar geçen sürenin üstel dağılımı izlediği ortaya çıkar.
Bunu göstermek için takımların yine maç başına ortalama toplam 6 gol attığını varsayalım ve her saniyedeki gol olasılığını hesaplayalım.
n = 3600
m = 6
p = m / 3600
p
0.0016666666666666668
Aşağıdaki fonksiyon n saniyeyi simüle eder ve ilk golün zamanını bulmak için argmax fonksiyonunu kullanır.
def simulate_first_goal(n, p):
return flip(n, p).argmax()
Bu yöntem çalışır; çünkü flip fonksiyonunun sonucu 1 ve 0’lardan oluşan bir dizidir, dolayısıyla en büyük değer neredeyse her zaman 1’dir. Dizide en az bir gol varsa argmax ilkinin indeksini döndürür. Hiç gol yoksa 0 döndürür; ancak bu durum yeterince seyrek gerçekleştiği için onu göz ardı edeceğiz.
1001 maçı simüle etmek ve ilk gole kadar geçen sürelerin listesini oluşturmak için simulate_first_goal fonksiyonunu kullanacağız.
np.random.seed(3)
first_goal_times = [simulate_first_goal(n, p) for i in range(1001)]
mean = np.mean(first_goal_times)
mean
np.float64(597.7902097902098)
İlk gole kadar geçen ortalama süre 600 saniyeye, yani 10 dakikaya yakındır. Bu sonuç anlamlıdır: Altmış dakikalık bir maçta 6 gol bekliyorsak ortalama olarak her 10 dakikada bir gol bekleriz.
n büyük ve p küçük olduğunda, ilk gole kadar beklenen sürenin üstel dağılımı izlediği matematiksel olarak gösterilebilir.
Simülasyon çok sayıda benzersiz zaman değeri ürettiği için dağılımları PMF yerine CDF ile karşılaştıracağız. Üstel dağılımın CDF’sini hesaplamak kolaydır.
def exponential_cdf(x, lam):
"""Üstel CDF'yi hesaplar.
x: float veya float dizisi
lam: Hız parametresi
Döndürür: Birikimli olasılığın float veya NumPy dizisi
"""
return 1 - np.exp(-lam * x)
lam değeri, bir zaman birimi başına ortalama olay sayısıdır; bu örnekte saniye başına gol sayısıdır. Simülasyon sonuçlarının ortalamasını kullanarak lam değerini hesaplayabiliriz.
lam = 1 / mean
lam
np.float64(0.0016728276636563566)
Bu fonksiyonu bir zaman değerleri aralığıyla çağırırsak ilk gol zamanlarının dağılımına yaklaşabiliriz. NumPy’nin linspace fonksiyonu eşit aralıklı değerlerden oluşan bir dizi üretir; bu örnekte 0 ve 3600 dâhil olmak üzere bu aralıkta 201 değer hesaplar.
from empiricaldist import Cdf
ts = np.linspace(0, 3600, 201)
ps = exponential_cdf(ts, lam)
cdf_expo = Cdf(ps, ts, name="üstel model")
Aşağıdaki şekil simülasyon sonuçlarını az önce hesapladığımız üstel dağılımla karşılaştırır.
cdf_sim = Cdf.from_seq(first_goal_times, name="simülasyon")
cdf_expo.plot(ls=":", color="gray")
cdf_sim.plot()
decorate(xlabel="İlk golün zamanı (saniye)", ylabel="CDF")
Üstel model simülasyon sonuçlarına çok iyi uyar; ancak daha güçlü bir sınama, gerçek verilerde nasıl sonuç verdiğine bakmaktır.
Aşağıdaki döngü bütün maçların sonuçlarını okur, ilk golün zamanını alır ve sonucu bir listede saklar. Hiç gol atılmamışsa listeye nan ekler.
filename = "nhl_2023_2024.hdf"
with pd.HDFStore(filename, "r") as store:
keys = store.keys()
firsts = []
for key in keys:
times = pd.read_hdf(filename, key=key)
if len(times) > 0:
firsts.append(times[0])
else:
firsts.append(np.nan)
Gol atma hızını tahmin etmek için nan değerlerini göz ardı ederek zamanların ortalamasını hesaplayan nanmean fonksiyonunu kullanabiliriz.
lam = 1 / np.nanmean(firsts)
lam
np.float64(0.0015121567467720825)
Şimdi verilerle aynı gol atma hızına sahip üstel dağılımın CDF’sini hesaplayabiliriz.
ps = exponential_cdf(ts, lam)
cdf_expo = Cdf(ps, ts, name="üstel model")
Verilerin CDF’sini hesaplamak için nan değerlerini sona dâhil eden dropna=False argümanını kullanacağız.
cdf_firsts = Cdf.from_seq(firsts, name="veri", dropna=False)
cdf_firsts.tail()
| probs | |
|---|---|
| 3286.0 | 0.996951 |
| 3581.0 | 0.997713 |
| NaN | 1.000000 |
Aşağıdaki şekil üstel dağılımı verilerin dağılımıyla karşılaştırır.
cdf_expo.plot(ls=":", color="gray")
cdf_firsts.plot()
decorate(xlabel="İlk golün zamanı (saniye)", ylabel="CDF")
Veriler bazı yerlerde modelden sapar; ilk 1000 saniyede modelin öngördüğünden daha az gol varmış gibi görünüyor. Yine de model verilere iyi uyar.
Bu modellerin temel varsayımı – goller için Poisson modeli ve süreler için üstel model – maçın herhangi bir saniyesinde gol atılma olasılığının eşit olmasıdır. Bunun doğru olup olmadığını bir hokey taraftarına sorarsanız “hayır” der ve haklı olur; gerçek dünya bu tür varsayımları pek çok yönden ihlal eder. Buna rağmen kuramsal dağılımlar çoğu zaman gerçek verilere dikkate değer ölçüde iyi uyar.
5.4. Normal dağılım#
Gerçek dünyada ölçtüğümüz birçok şey, Gauss dağılımı veya “çan eğrisi” olarak da bilinen normal dağılımı izler. Bu dağılımların nereden geldiğini görmek için dev kabakların büyümesini açıklayan bir model ele alalım. Her gün hava kötüyse kabağın 1 libre, ılımlıysa 2 libre, iyiyse 3 libre aldığını varsayalım. Ayrıca her gün havanın kötü, ılımlı veya iyi olma olasılığının eşit olduğunu kabul edelim.
Aşağıdaki fonksiyonla bu modeli n gün boyunca simüle edip toplam ağırlık artışını döndürebiliriz.
def simulate_growth(n):
choices = [1, 2, 3]
gains = np.random.choice(choices, n)
return gains.sum()
NumPy’nin random modülü, bir değerler dizisinden n rastgele seçim üreten choice fonksiyonunu sağlar; bu örnekte dizi choices değişkenidir.
Şimdi farklı yerlerde ve farklı hava koşullarında 1001 kişinin dev kabak yetiştirdiğini varsayalım. Büyüme sürecini 100 gün boyunca simüle edersek 1001 ağırlıktan oluşan bir liste elde ederiz.
sim_weights = [simulate_growth(100) for i in range(1001)]
m, s = np.mean(sim_weights), np.std(sim_weights)
m, s
(np.float64(199.37062937062936), np.float64(8.388630840376777))
Ortalama 200 libreye, standart sapma ise 8 libreye yakındır. Ağırlıkların normal dağılımı izleyip izlemediğini görmek için aşağıdaki fonksiyonu kullanacağız. Bu fonksiyon bir örneklem alır ve örneklemle aynı ortalama ve standart sapmaya sahip normal dağılımı temsil eden bir Cdf oluşturur; dağılımı ortalamanın dört standart sapma altından dört standart sapma üstüne kadar değerlendirir.
from scipy.stats import norm
def make_normal_model(data):
m, s = np.mean(data), np.std(data)
low, high = m - 4 * s, m + 4 * s
qs = np.linspace(low, high, 201)
ps = norm.cdf(qs, m, s)
return Cdf(ps, qs, name="normal model")
Fonksiyonu aşağıdaki gibi kullanırız.
cdf_model = make_normal_model(sim_weights)
Şimdi simülasyon sonuçlarının dağılımını temsil eden bir Cdf oluşturabiliriz.
cdf_sim_weights = Cdf.from_seq(sim_weights, name="simülasyon")
Dağılımları karşılaştırmak için aşağıdaki fonksiyonu kullanacağız. cdf_model ve cdf_data birer Cdf nesnesidir. xlabel bir string, options ise cdf_data çiziminin biçimini denetleyen seçenekler sözlüğüdür.
def two_cdf_plots(cdf_model, cdf_data, xlabel="", **options):
cdf_model.plot(ls=":", color="gray")
cdf_data.plot(**options)
decorate(xlabel=xlabel, ylabel="CDF")
Sonuçlar aşağıda gösterilmiştir.
two_cdf_plots(cdf_model, cdf_sim_weights, xlabel="Ağırlık (libre)")
Normal model ağırlıkların dağılımına çok iyi uyar. Genel olarak yeterli sayıda rastgele etmeni topladığımızda toplam normal dağılımı izleme eğilimindedir. Bu, 14. Bölümde yeniden ele alacağımız Merkezi Limit Teoremi’nin bir sonucudur.
Önce normal dağılımın gerçek verilere ne kadar iyi uyduğuna bakalım. Örnek olarak Ulusal Aile Büyümesi Araştırması’ndaki (NSFG) doğum ağırlıklarının dağılımını inceleyeceğiz. Verileri read_fem_preg ile okuyup doğum ağırlıklarını libre cinsinden kaydeden totalwgt_lb sütununu seçebiliriz.
Aşağıdaki hücreler veri dosyalarını indirir ve verileri okumak için gereken statadict paketini kurar.
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")
try:
import statadict
except ImportError:
%pip install statadict
import nsfg
preg = nsfg.read_fem_preg()
birth_weights = preg["totalwgt_lb"].dropna()
Doğum ağırlıklarının ortalaması yaklaşık 7,27 libre, standart sapması ise 1,4 libredir; ancak gördüğümüz üzere bu veri kümesinde muhtemelen hatalı olan bazı aykırı değerler vardır.
m, s = np.mean(birth_weights), np.std(birth_weights)
m, s
(np.float64(7.265628457623368), np.float64(1.40821553384062))
Aykırı değerlerin tahmin edilen ortalama ve standart sapma üzerindeki etkisini azaltmak için SciPy’nin trimboth fonksiyonuyla en yüksek ve en düşük değerleri kaldıracağız.
from scipy.stats import trimboth
trimmed = trimboth(birth_weights, 0.01)
m, s = np.mean(trimmed), np.std(trimmed)
m, s
(np.float64(7.280883100022579), np.float64(1.2430657948614345))
Kırpılmış verilerde ortalama biraz, standart sapma ise belirgin ölçüde düşüktür. Normal modeli oluşturmak için kırpılmış verileri kullanacağız.
cdf_model = make_normal_model(trimmed)
Bunu verilerin Cdf’siyle karşılaştıracağız.
cdf_birth_weight = Cdf.from_seq(birth_weights, name='veri')
two_cdf_plots(cdf_model, cdf_birth_weight, xlabel="Doğum ağırlığı (libre)")
Normal model 5 librenin altı dışında verilere iyi uyar. Bu aralıkta veri dağılımı modelin solundadır; başka bir deyişle en hafif bebekler normal dağılım altında beklediğimizden daha hafiftir. Gerçek dünya çoğunlukla basit matematiksel modellerden daha karmaşıktır.
5.5. Lognormal dağılım#
Önceki bölümde kabakların hava durumuna bağlı olarak günde 1-3 libre büyüdüğü varsayımıyla büyümeyi simüle ettik. Bunun yerine büyümenin mevcut ağırlıkla orantılı olduğunu varsayalım; böylece büyük kabaklar küçüklerden günde daha fazla ağırlık kazanır. Bu muhtemelen daha gerçekçi bir modeldir.
Aşağıdaki fonksiyon bu tür orantılı büyümeyi simüle eder: Hava kötüyse kabak ağırlığının %3’ünü, ılımlıysa %5’ini, iyiyse %7’sini kazanır. Yine herhangi bir gündeki hava durumlarının eşit olasılıklı olduğunu varsayacağız.
def simulate_proportionate_growth(n):
choices = [1.03, 1.05, 1.07]
gains = np.random.choice(choices, n)
return gains.prod()
Bir kabak ağırlığının %3’ünü kazanırsa son ağırlık, ilk ağırlıkla 1,03 çarpanının çarpımıdır. Bu nedenle rastgele çarpanlar seçip bunları çarparak 100 gün sonraki ağırlığı hesaplayabiliriz.
1001 kabağı simüle etmek ve ağırlıklarını saklamak için bu fonksiyonu 1001 kez çağıracağız.
sim_weights = [simulate_proportionate_growth(100) for i in range(1001)]
np.mean(sim_weights), np.std(sim_weights)
(np.float64(130.80183363824722), np.float64(20.956047434921466))
Ortalama ağırlık yaklaşık 131 libre, standart sapma ise yaklaşık 21 libredir. Dolayısıyla bu modeldeki kabaklar önceki modeldekilerden daha küçük fakat daha değişkendir.
Bunların lognormal dağılımı izlediği matematiksel olarak gösterilebilir; yani ağırlıkların logaritmaları normal dağılımı izler. Bunu kontrol etmek için ağırlıkların logaritmalarını, bunların ortalamasını ve standart sapmasını hesaplayacağız. Herhangi bir tabandaki logaritmayı kullanabilirdik; sonuçları yorumlamayı kolaylaştırdığı için 10 tabanını kullanacağım.
log_sim_weights = np.log10(sim_weights)
m, s = np.mean(log_sim_weights), np.std(log_sim_weights)
m, s
(np.float64(2.1111299372609933), np.float64(0.06898607064749826))
Şimdi logaritmaların dağılımını aynı ortalama ve standart sapmaya sahip normal dağılımla karşılaştıralım.
cdf_model = make_normal_model(log_sim_weights)
cdf_log_sim_weights = Cdf.from_seq(log_sim_weights, name="simülasyon")
two_cdf_plots(
cdf_model, cdf_log_sim_weights, xlabel="Kabak ağırlığı (log10 libre)"
)
Model, beklediğimiz üzere simülasyon sonucuna çok iyi uyar.
İnsanlar da kabaklar gibiyse ve ağırlıktaki yıllık değişim mevcut ağırlıkla orantılıysa yetişkin ağırlıklarının lognormal dağılımı izlemesini bekleyebiliriz. Bunu sınayalım.
Ulusal Kronik Hastalıkları Önleme ve Sağlığı Geliştirme Merkezi, Davranışsal Risk Faktörü Gözetim Sistemi’nin (BRFSS) bir parçası olarak yıllık bir araştırma yürütür. Merkez 2008’de 414.509 katılımcıyla görüşerek demografik özellikleri, sağlık durumları ve sağlık riskleri hakkında sorular sormuştur. Toplanan veriler arasında 398.484 katılımcının ağırlıkları da vardır. Verileri indirme talimatları bu bölümün notebook’unda yer alır.
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/CDBRFS08.ASC.gz")
thinkstats modülü, BRFSS verilerini okuyup bir Pandas DataFrame döndüren bir fonksiyon sağlar.
from thinkstats import read_brfss
brfss = read_brfss()
Yetişkin ağırlıkları kilogram cinsinden wtkg2 sütununda kaydedilmiştir.
adult_weights = brfss["wtkg2"].dropna()
m, s = np.mean(adult_weights), np.std(adult_weights)
m, s
(np.float64(78.99245299685809), np.float64(19.546132387397257))
Ortalama yaklaşık 79 kg’dır. Logaritmaları hesaplamadan önce ağırlıkların normal dağılımı izleyip izlemediğine bakalım.
cdf_model = make_normal_model(adult_weights)
cdf_adult_weights = Cdf.from_seq(adult_weights, name="yetişkin ağırlığı")
two_cdf_plots(cdf_model, cdf_adult_weights, xlabel="Yetişkin ağırlığı (kilogram)")
Normal dağılım bazı amaçlar için bu verilere yeterince iyi bir model olabilir; ancak daha iyisini yapıp yapamayacağımıza bakalım.
Aşağıda logaritmik dönüşüm uygulanmış ağırlıkların dağılımı ile aynı ortalama ve standart sapmaya sahip normal model gösterilmiştir.
log_adult_weights = np.log10(adult_weights)
cdf_model = make_normal_model(log_adult_weights)
cdf_log_adult_weights = Cdf.from_seq(log_adult_weights, name="log yetişkin ağırlığı")
two_cdf_plots(cdf_model, cdf_log_adult_weights, xlabel="Yetişkin ağırlığı (log10 kg)")
Normal model logaritmalara, ağırlıkların kendisine uyduğundan daha iyi uyar. Bu, orantılı büyümenin ağırlık artışı için toplamsal büyümeden daha iyi bir model olduğunu düşündürür.
5.6. Neden modelleme?#
Bu bölümün başında gerçek dünyadaki birçok olgunun kuramsal dağılımlarla modellenebileceğini söyledim. Fakat bunun neden önemli olduğu açık olmayabilir.
Bütün modeller gibi kuramsal dağılımlar da soyutlamadır; yani ilgisiz kabul edilen ayrıntıları dışarıda bırakırlar. Örneğin gözlenen bir dağılım, ölçüm hataları veya örnekleme özgü tuhaflıklar içerebilir; kuramsal modeller bu kendine özgü özellikleri göz ardı eder.
Kuramsal modeller aynı zamanda bir veri sıkıştırma biçimidir. Bir model veri kümesine iyi uyduğunda az sayıda değer büyük miktarda veriyi özetleyebilir.
Doğal bir olgudan elde edilen verilerin kuramsal bir dağılıma uyması bazen şaşırtıcıdır; ancak bu gözlemler fiziksel sistemlere ilişkin içgörü sağlayabilir. Bazen gözlenen bir dağılımın neden belirli bir biçime sahip olduğunu açıklayabiliriz. Örneğin önceki bölümde yetişkin ağırlıklarının lognormal dağılımla iyi modellendiğini gördük; bu da ağırlıktaki yıllık değişimlerin mevcut ağırlıkla orantılı olabileceğini düşündürür.
Ayrıca kuramsal dağılımlar, 14. Bölümde göreceğimiz gibi matematiksel analize elverişlidir.
Ancak bütün modellerin kusurlu olduğunu unutmamak gerekir. Gerçek dünya verileri hiçbir kuramsal dağılıma kusursuz biçimde uymaz. Bazen veriler modeller tarafından üretiliyormuş gibi konuşulur; örneğin insan boylarının dağılımının normal veya gelir dağılımının lognormal olduğu söylenebilir. Sözcüğü sözcüğüne alındığında bu iddialar doğru olamaz; gerçek dünya ile matematiksel modeller arasında her zaman farklar vardır.
Modeller, gerçek dünyanın ilgili yönlerini yakalayıp gereksiz ayrıntıları dışarıda bıraktıklarında yararlıdır. Ancak neyin ilgili ya da gereksiz olduğu, modeli hangi amaçla kullanacağınıza bağlıdır.
5.7. Sözlük#
binom dağılımı: İsabet ve ıskalamalardan oluşan bir dizideki başarı veya isabet sayısını modellemek için sık kullanılan kuramsal dağılım.
Poisson dağılımı: Bir zaman aralığında gerçekleşen olayların sayısını modellemek için sık kullanılan kuramsal dağılım.
üstel dağılım: Olaylar arasındaki süreyi modellemek için sık kullanılan kuramsal dağılım.
normal dağılım: Simetrik ve çan biçimli bir eğri izleyen verileri modellemek için sık kullanılan kuramsal dağılım.
lognormal dağılım: Sağa çarpık, çan benzeri bir eğri izleyen verileri modellemek için sık kullanılan kuramsal dağılım.
5.8. Alıştırmalar#
5.8.1. Alıştırma 5.1#
NSFG katılımcı dosyasındaki numfmhh sütunu, her katılımcının hanesindeki “aile üyesi sayısını” kaydeder. Dosyayı read_fem_resp ile okuyabilir, görüşme sırasında 25 yaşında veya daha büyük olan katılımcıları query ile seçebiliriz.
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/2002FemResp.dct")
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/2002FemResp.dat.gz")
from nsfg import read_fem_resp
resp = read_fem_resp()
older = resp.query("age >= 25")
num_family = older["numfmhh"]
Bu daha yaşlı katılımcılar için numfmhh değişkeninin Pmf’sini hesaplayın ve aynı ortalamaya sahip Poisson dağılımıyla karşılaştırın. Poisson modeli verilere ne kadar iyi uyuyor?
5.8.2. Alıştırma 5.2#
Bu bölümün önceki kısımlarında bir hokey maçında ilk gole kadar geçen sürenin üstel dağılımı izlediğini gördük. Gol atma modelimiz doğruysa önceki golün üzerinden ne kadar zaman geçtiğinden bağımsız olarak her an gol atılma olasılığı eşittir. Bu doğruysa goller arasındaki sürenin de üstel dağılımı izlemesini bekleriz.
Aşağıdaki döngü hokey verilerini yeniden okur; bir maçta birden fazla gol varsa ardışık goller arasındaki süreleri hesaplar ve bunları bir listede toplar.
filename = "nhl_2023_2024.hdf"
with pd.HDFStore(filename, "r") as store:
keys = store.keys()
intervals = []
for key in keys:
times = pd.read_hdf(filename, key=key)
if len(times) > 1:
intervals.extend(times.diff().dropna())
Gözlenen aralıklarla aynı ortalamaya sahip üstel dağılımın CDF’sini exponential_cdf ile hesaplayın ve bu modeli verilerin CDF’siyle karşılaştırın.
5.8.3. Alıştırma 5.3#
İnsan boylarının dağılımı normal dağılıma mı, yoksa lognormal dağılıma mı daha çok benzer? Bunu anlamak için BRFSS’deki boy verilerini aşağıdaki gibi seçebiliriz:
adult_heights = brfss["htm3"].dropna()
m, s = np.mean(adult_heights), np.std(adult_heights)
m, s
(np.float64(168.82518961012298), np.float64(10.352640156455918))
Bu değerlerin CDF’sini hesaplayın ve aynı ortalama ile standart sapmaya sahip normal dağılımla karşılaştırın. Ardından boyların logaritmalarını hesaplayın ve logaritmaların dağılımını aynı ortalama ve standart sapmaya sahip normal dağılımla karşılaştırın. Görsel karşılaştırmaya göre hangi model verilere daha iyi uyuyor?
Think Stats: Exploratory Data Analysis in Python, 3rd Edition
Telif hakkı 2024 Allen B. Downey
Kod lisansı: MIT Lisansı
Metin lisansı: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International