Doç. Dr. Özgür Baştürk
Ankara Üniversitesi, Astronomi ve Uzay Bilimleri Bölümü
obasturk at ankara.edu.tr
http://ozgur.astrotux.org
Bu derste neler öğreneceksiniz?¶
İstatistiksel Yöntemler, Tahmin ve Çıkarım - Frekansçı Bakış¶
- Dağılımlar
- Bir Dağılım Foksiyonunun Elemanları
- Örnek Dağılım ve Ana Dağılım
- Dağılım Türleri
- Python'da Dağılımlar
- Merkezi Limit Teoremi
- Hipotez Testi
- t Dağılımı
- Serbestlik Derecesi
- Modelleme
- Bir Modelin Başarımına İlişkin Parametreler
- $\chi^2$ Dağılımı ve $\chi^2$ Testi
- F Dağılımı ve F Testi
- Durbin Watson İstatistiği İle Model Başarımı Testi
- Akaike Bilgi Kriteri
- Bayesian Bilgi Kriteri
- Kaynaklar
Dağılımlar¶
Olasılık teorisi terminolojisi yerleşmiş olmakla birlikte aynı şeyi ifade eden pek çok kavram bir arada, birbirlerinin yerine ve hatta bazen yanlış şekilde kullanılabilmektedir. Temelde dağılım (ing. distribution) bir veri setinin (bir örnek ya da onun seçildiği popülasyon) aldığı ya da alabileceği tüm olası değerleri ve bu değerlerin ne sıklıkta ya da hangi olasılıkla gerçekleştiğini gösteren fonskiyon, grafik, tablo ya da listedir. Aşağıda istatistiksel bir dağılımı ifade etmek üzere literatürde sıkça kullanılan kavramlara örnekler verilmiştir.
Bir Dağılım Foksiyonunun Elemanları¶
Bir dağılım fonksiyonunu yapısını anlatmak üzere kullanılan elemanlar temelde üç ölçüt grubu altında toplanır.
Merkezi Eğilim Ölçütleri
- Ortalama
- Mod
- Ortanca (medyan)
Dağılım Ölçütleri
- Standart Sapma
- Varyans
- Aralık (ing. range)
Şekil Ölçütleri
- Çarpıklık (ing. Skewness)
- Basıklık (Kurtosis)
Ortalama Değer / Beklenen Değer¶
Bir dağılımın ortalama değeri (ing. average (continuous), ing. mean (discrete)) tüm değerlerin toplamının (sürekli değişken durumunda integralinin) değer sayısına (sürekli durumda aralık uzunluğuna) bölümü ile bulunur. Ağırlıklı ortalama verilmek istenen durumda bu ağırlıkların toplamı 1 olacak şekilde düzenlenir. Ortalama değer, bu dağılımı oluşturan ölçümlerin doğruluğunu da belirler. Ölçümler ortalamadan ne kadar uzaksa onlardan o kadar şüphe duyulur.
Bir ölçüm ya da sonucun beklenen değeri (ing. expected value) ise teorik ifadelerle hesaplanan değeri ya da biliniyorsa popülasyonun ortalama değeridir. Herhangi bir ölçüm yaparken teorik olarak beklentiyle ya da poülasyonun ortalamasıyla uyum beklenmesi doğaldır. Bir örnek ölçüm için beklenen değerden farklılık gösteren ortalama değerler sistematik hatalardan kaynaklanabileceği gibi, beklenen değerin gerçek değer olmamasından (örn. teorik hesapların eksik/yanlış olmasından) ya da yeterli ölçümün yapılmamasından kaynaklanabilir.
Aritmetik ortalama;
$$ \bar{x} = \frac{1}{N} \sum\limits_{i=1}^n x_{i} $$
$w_i$ i noktasının ağırlığını göstermek üzere, ağırlıklı ortalama ise
$$ \bar{x} = \frac{1}{N} \sum\limits_{i=1}^n \frac{x_{i} / w_i^2}{1 / w_i^2} $$
il verilir.
Standart Sapma ve Varyans¶
Varyans (popülasyon: $s^2$, örnek: $\sigma^2$) ve standart sapma (popülasyon: $s$, örnek: $\sigma$), dağılım değerlerinin ortalama değer etrafında ne kadar çok saçılmış olduğunu belirler. Dolayısıyla yapılan ölçümlerin ya da elde edilen sonuçların ne kadar hassas olduğunu belirlemek ve bu değerlere dayalı tahminler yaparaken söz konusu tahminler üzerindeki belirsizliği ifade etmek üzere kullanılır.
$$ \sigma^2 = \frac{\sum_{i=1}^{N} (x_i - \mu)^2}{N - 1} $$
Varyans ve standart sapmanın bazı özellikleri:
- Varyans ve standart sapma değerleri negatif olamaz.
- Tüm ölçümler aynı değere sahip ise varyans sıfır değerini alır.
- Varyans değeri, dağılımın konumundan bağımsızdır. Tüm değerler aynı miktarda kaydırıldığında varyans değişmez.
- Varyansın birimi, ölçülen değerin biriminin karesidir. Standart sapmanın birimi ise ölçülen değerin birimidir.
Diğer Merkezi Eğilim Ölçütleri¶
Bir dağılımın modu (ing.mode) söz konusu değişkenin ölçümlerde en çok tekrar eden (en sık rastlantılan) değeridir.
Bir dağılımın medyanı (ya da ortancası) ise söz konusu değişkene ilişkim tüm ölçüm değerleri küçükten büyüğe sıralandığında tam ortada kalan değeridir. Eğer çift sayıda ölçüm yapılmışsa bu kez ortada kalan iki ölçüm değerinin ortalaması alınır.
Özellikle simetrik olmayan dağılımlar söz konusu olduğunda, tercihe göre beklenen değer olarak ortalama değer yerine mod ya da medyan değeri kullanılabilmektedir. Bu tercih gözlenen/ölçülen olgunun türüne göre yapılabilmektedir. Ortalama yerine mod ya da medyan bu asimetriye yol açan ekstrem (marjinal ya da uç) dğeerler söz konusu olduğunda merkezi dağılımı belirlemek adına tercih sebebi olabilmektedir.
Dağılıma bağlı olarak mod, medyan ve ortalama değerleri aynı olabileceği gibi farklı da olabilir.
Çarpıklık¶
Çarpıklık (ing. skewness) bir dağılımın asimetrisinin ölçütüdür. Bu ölçütün basit ve standart bir matematiksel ifadesi yoktur.
Temel olarak iki tür çarpıklık vardır.
- Negatif Çarpıklık (ing. left skewed, negatively skewed): dağılımın sol kuyruğu uzundur; dağılımın önemli bir çoğunluğu sağ tarafta toplanmıştır.
- Pozitif Çarpıklık (ing. right skewed, positively skewed): dağılımın sağ kuyruğu uzundur; dağılımın önemli bir çoğunluğu sol tarafta toplanmıştır.
Çarpıklığı bir normal dağılıma göre 0'dan farklı dağılımlar için ortalama ve ortanca birbirinden farklı olabilir.
Basıklık¶
Basıkılık bir dağılımın normal dağılıma göre ne kadar geniş ya da ‘kuyruklu’ olduğunun ölçütüdür. $kurtosis$ gibi basıklık ölçütleriyle ifade edilir.
Örnek Dağılım ve Ana Dağılım¶
Yapılan gözlemlerin sayısının arttırılması gözlemlerin oluşturduğu dağılımı, ilgili olgunun gerçek dağılımına daha fazla yaklaştıracaktır. Ancak herhangi bir değişkenin olası tüm değerlerini alabileceği çok fazla sayıda (sürekli durumda sonsuz) gözlem yapmanın mümkün olamayabilmesi sebebiyle yapılan gözlemler çoğu zaman gerçek dağılımın bir örneği niteliğini taşır.
Bir dağılımı oluşturmak için olası tüm değerlerin kullanılması kabulu, elde edilen dağılımın ana dağılım (parent distribution) olarak kabul edilmesi anlamına gelmektedir. Gerçekte herhangi bir olayın tam olarak ne tür bir dağılım gösterdiğini çoğu zaman bilemeyiz. Ancak bu dağılımı yeterli hassasiyette temsil ettiği kabul edilen matematiksel bir fonksiyonu, olgunun ana dağılımı olarak kabul edebiliriz.
Bu kabulden sonra, yapılan gözlemlerin oluşturduğu ve sonlu sayıdaki değerler ile üretilen örnek dağılımı (ing. sample distribution) kullanarak ana dağılıma ilişkin parametreleri elde edebilir ya da ana dağılımın geçerliliğini sorgulayabiliriz.
Sıklık Dağılımı¶
Sıklık ya da frekans dağılımları (frequency distribution) Bir örnek grubunda (örneğin bir dizi ölçümde) her bir sonucun (ölçümün) kaç kez gerçekleştiğini ya da tekrarlandığını gösterir. Genellikle bir tablo ya da histogramla ifade edilir.
Olasılık Dağılımları¶
Olasılık dağılımı (probability distribution), bir örnek grubunda her bir sonucun (ölçüm ya da deney sonuçları gibi) gerçekleşme olasılıklarının ifadesidir. Olasılık dağılımları değişkenin sürekli (ing. continuous) ya da süreksiz olmasına bağlı olarak iki farklı şekilde ifade edilirler.
Olasılık Dağılım Fonksiyonlarının Özellikleri¶
Rastgele Değişken (ing. random variable): Olası değerleri rastgele bir olayın sayısal sonuçları olan bir değişkendir. İki tür rastgele değişken vardır, süreksiz (kesikli) (ing. discrete) ve sürekli (ing. continuous) değişkenler.
Süreksiz bir rastgele değişken, sadece sayılabilir değerler alabilen ve böylece nicelleştirilebilen bir değişkendir. Örneğin, adil bir zarın her atılıştaki değeri rastgele bir $x$ değişkeni olarak tanımlanabilir. Bu durumda $x$, $[1,2,3,4,5,6]$ değerlerini alabilir. Bu nedenle de $x$, süreksiz rastgele bir değişkendir.
Süreksiz rastgele bir değişkenin olasılık dağılımı, olası değerlerinin her biri ile ilişkili olasılıkların bir listesidir. Buna olasılık fonksiyonu ya da daha doğru bir ifadeyle Olasılık Kütle Fonksiyonu (ing. Probability Mass Function, PMF) denir.
Rastgele bir $X$ değişkeninin k farklı değer alabileceğini varsayalım. $X = x_i$ 'ye karşılık gelen $P (X = x_i) = p_i$ olasılıkları aşağıdaki koşulları sağlar:
- $0 \lt p_i \lt 1$, $i = 1, 2, 3, ...$
- $p_1 + p_2 + ... + p_k = 1$
Sürekli rastgele değişken, sonsuz sayıda olası değeri alabilen değişkendir. Örneğin, bir sınıftaki öğrencilerin boyları rastgele bir $x$ değişkeni olarak tanımlandığında sürekli bir rastgele değişkendir Sürekli rastgele değişken bir değerler aralığı içinde tanımlandığından, o aralıkta olma olasılığı, eğrinin (veya integralin) altındaki alanla temsil edilir.
Olasılık Yoğunluk Fonksiyonu (ing. Probability Density Function, PDF olarak bilinen sürekli rasgele değişkenin olasılık dağılımı, sürekli değerleri alan fonksiyonlardır. Herhangi bir tek değeri gözlemleme olasılığı, 0'a eşittir, çünkü rastgele değişken tarafından alınabilecek değerlerin sayısı sonsuzdur ($p(x_i) = 0$. Olasılık yoğunluk fonksiyonu için $P (a \lt x \lt b) = p$ olasılıkları aşağıdaki koşulları sağlar:
Olasılık her zaman pozitiftir ($x$'in bütün değerleri için $p \gt 0$).
$p(x)$'in altında kalan 1'dir.
Ayrıca, tüm rastgele değişkenlerin (sürekli ya da değil) birikimli (kümülatif) dağılım karşılğı da vardır. Birikimli dağılım fonksiyonu, her $x$ değeri için rastgele değişken $X$'in $x$'den küçük veya ona eşit olma olasılığını veren bir fonksiyondur ($P(X \le x) = p(x))$. Kesikli rastgele bir değişken için olasılık dağılımları toplanarak birikimli (kümülatif) dağılım fonksiyonu bulunur.
Olasıık Kütle Fonksiyonu¶
Olasılık Kütle Fonksiyonu (Probability Mass Function, PMF), sadece bazı değerleri alabilen (ing. discrete) bir değişken, deney sonucu, gözlem ya da ölçümün her bir olası değerini alma olasılığını veren fonksiyondur. Normalize olasılık dağılım fonksiyonlarında her bir ölçüme ilişkin olasılık değerlerinin toplamı 1 değerine normalize edilir. Normalize olasılık kütle fonksiyonlarında tüm olası değerlerin olasılıklarının toplamı 1’i verir.
Olasılık Yoğunluğu Fonksiyonu¶
Olasılık Yoğunluğu Fonksiyonu (Probability Density Function, PDF) sürekli (ing. continuous) bir değişkenin (doğası sürekli, yani bir aralıktaki tüm değerleri alabilen, deney sonucu ya da ölçüm) herhangi bir aralıkta değer alma olasılığını veren fonksiyondur. Normalize olasılık dağılım fonksiyonlarında her bir ölçüme ilişkin olasılık değerlerinin toplamı 1 değerine normalize edilir. PDF grafik olarak ifade edildiğinde eğrinin altında kalan toplam alan normalize olasılık fonskiyonu için 1’dir. Herhangi bir değer aralığı için toplam olasılık eğrinin o aralık için altında kalana eşit olur.
Eğer x değişkeni sürekli bir aralıktan değer alabiliyorsa x’in herhangi bir değerinin olasılığı verilemez. Örneğin bir zar atışında 1 ile 6 arasındaki her bir sayının gelme olasılığı 1 / 6 iken, bir ölçümde 0 ile 1 arasındaki herhangi bir reel sayıyı alabilen x’in e, π ya da 0.1 olma olasılığı belirlenemez. Ancak x’in hangi aralık dahilinde bir değer alabileceğinin olasılığı belirlenebilir ve tüm olası değerler için bu olasılıklar Olasılık Yoğunluğu Fonksiyonu (PDF) ile ifade edilebilir.
Kümülatif ya da Birikimli Dağılım Fonksiyonu¶
Kümülatif (Birikimli) Dağılım Fonksiyonu (ing. Cumulative Distribution Function, CDF) herhangi bir olasılık dağılım fonksiyonunun sahip olduğu değerlerin toplanarak temsil edildiği fonksiyonlardır. Olasılık değerlerinin tamamı 1 ya da %100 ihtimale sahip olduğu için 1 değerine yakınsamaktadır.
Dağılım Türleri¶
Herhangi bir rastgele değişkenin alabileceği değerlerin olasılıkları farklı olasılık kütle ya da olasılık yoğunluk fonksiyonları kullanılarak belirlenebilir. Bu tamamen ele alınan rastgele değişkenin doğasına bağlıdır. Doğada her rastgele değişken en bilinen dağılım türü olan normal dağılıma uymayabilir. Ancak normal dağılım (ya da Gaussian dağılım) özellikle ölçüm rastgeleliğini ifade etmekte ve çok sayıda ölçüm olduğunda bu rastgeliliğin parametrelerini belirlemede çok kullanışlı olduğu için en sık başvurulan dağılımdır. Ancak doğadaki rastgele değişkenelerin olasılık dağılımlarını ifade etmek için başka fonksiyonlar da bulunmaktadır. Bu bölümde bu fonksiyonlardan en sık kullanılanları ele alınmıştır.
Tekdüze Dağılımlar¶
Bir değişkenin tekdüze dağılımı (ing. uniform distribution) onun bir aralık içerisindeki değerlerin tamamının aynı olasılıkla alabileceği sürekli bir dağılımdır. $x = a$ ile $x = b$ arasında tanımlı rastgele bir $x$ değişekin için tekdüze dağılımın olasılık yoğunluk fonksiyonu aşağıdaki ifade ile verilir.
$$ f(x) = \frac{1}{b - a}, a \le x \le b $$, $$ f(x) = 0, x < a~,~x > b $$
Normal Dağılım ve Gauss Fonksiyonu¶
Gauss dağılımı olarak da bilinen Normal Dağılım, deneysel bütün bilimlerde yeri bulunan sürekli bir olasılık dağılımıdır. Normal dağılım, ortalama, $\mu$ ve standart sapma, $\sigma$) ile tanımlanan çan şeklindeki yoğunluk eğrisine sahiptir. Yoğunluk eğrisi (ing. Probability Density Function, PDF), merkezi ortalamaya yakın, yayılımı standart sapması ile belirlenen; ortalamaya yakın verilerin, ortalamadan çok uzaktaki verilere göre daha sık olduğunu gösteren simetrik bir fonksiyondur. Belirli bir $x$ noktasında ortalama $\mu$ ve standart sapma $\sigma$ ile normal yoğunluk eğrisinin olasılık dağılım fonksiyonu şu şekilde verilir:
$$ p(x)~dx = \frac{1}{\sigma\sqrt{2\pi}} e^{\frac{-1}{2} (\frac{x - \mu}{\sigma})^2}~dx $$
Standart Normal Dağılım¶
Standart normal dağılımda ($\mu = 0$, $\sigma = 1$) herhangi bir $x$ değişkeninin ortalamadan z = $\pm \Delta x$ kadar sapması olasılığı scipy.stats.norm paketi kümülatif dağılım fonksiyonu (ing. cumulative distribution function) cdf kullanılarak elde edilebilir. Matematiksel olarak,
$$ p(z)~dz = \frac{1}{\sigma\sqrt{2\pi}} e^{\frac{z^2}{2}}~dz$$
z-değeri $x$ rastgele değişkeninin herhangi bir $X$ değerinin ortalamadan kaç standart sapma uzakta olduğunun ölçüsüdür. Dolayısı ile yukarıdaki ifadeyle bir normal dağılım standardize edildiği zaman, ortalama değer 0 olur (ortalamadan 0 sapma uzaklıktadır). Standart normal dağılımın standart sapması da doğal olarak $1$'e eşittir. $z = -1$ 'e karşılık gelen bir $X$ değeri, $x$ değişkeninin ortalama değerinden $1$-standart sapma küçüktür.
Poisson Dağılımı¶
Poisson rasgele değişkeni nadir gerçekleşen bir olayın bir zaman aralığında kaç kez meydana geldiğini modellemek için kullanılır. Örneğin, gökbilimdeki CCD gözlemlerinde kısa poz süreleri dahilinde dedektör üzerinde toplanan foton sayıları ya da bir web sitesinde belirli aralıklarla ziyaret edilen kullanıcı sayısı bir Poisson süreci olarak düşünülebilir.
Bir olayın bir Poisson süreci olarak tanımlanması ve bu olayın n = 1, 2, 3, ... kez gerçekleşme olasılıklarını belirlemek için Poisson dağılımının kullanılması için 3 gerekli koşul bulunmaktadır.
- Ölçümler nadir görülen olayların sayısıdır.
- Tüm ölçümler birbirinden bağımsızdır.
- Ölçümlerin görülme sıklığı ilgili zaman aralığında değişim göstermemektedir.
Poisson dağılımı olayların meydana gelme oranı ($\mu$) cinsinden tanımlanır. Bir olay bir aralıkta $n = 0, 1, 2, … $ kez oluşabilir. Bir aralıktaki ortalama olay sayısı $\lambda$ olarak belirtilir. Lambda, oran parametresi (ing. rate parameter) olarak da adlandırılan olayın gerçekleşme hızıdır. k tane olayı belirli bir aralıkta gözleme olasılığı;
$$ P(k) = e^{-\lambda}\frac{\lambda^k}{k!} $$
Aslında Poisson dağılımı, Binom dağılımının, nadir gözlenen bir olayın (örneğin foton sayımı) olasılığının çok küçük olması ($p \lt\lt 1$) nedeniyle ortalama gerçekleşme sayısının ($\mu$) olası tüm olayların sayısından ($n$) çok daha küçük olduğu ($k = \mu \lt\lt n$) için özel bir durumudur. Poisson dağılımı ($p \lt\lt 1$) için Binom dağılımının bir yaklaşımı olduğundan ortalama etrafında asimetrik dağılan bir dağılımdır ve şekli 10 kez zar atma deneyinde x kez 1 gelme olasılığının dağılımına benzer.
Poisson dağılımının ortalama değeri aslında tanımında vardır. Ancak ortalama değer, $x$ için beklenen değer ($\bar{x}$) tanımından da çıkarılabilir. Standart sapması da ona eşit olduğundan Poisson dağılımı tek bir parametreyle karakterize edilir ($\mu$).
$\lambda \rightarrow \infty$ limit durumunda Poisson dağılımı Normal dağılıma yakınsar.
Bu parametre Poisson sürecine uyan bir olayın belirli bir aralıkta kaç kez gerçekleştiğni tanımlayan oran parametresidir..
Üstel Dağılım¶
Üstel dağılım (ing. exponential distribution), bir Poisson sürecindeki (sürekli ve bağımsız ("hafızasız") bir şekilde sabit bir ortalama oranında meydana gelen) olaylar arasındaki sürenin olasılık dağılım fonksiyonu olarak tanımlanır. gamma dağılımının özel bir durumu, süreksiz bir olasılık dağılım fonksiyonu olan geometrik dağılımın sürekli durumudur. Normal, binom, Poisson, gamma gibi dağılımların üyesi olduğu üstel dağılımlar ailesinin bir üyesidir. Oran (ing. rate) parametresi olarak adlandırılan bir $\lambda$ parametresine sahiptir ve bu parametre dağılımın azalma (artma) hızını belirler. Olasılık yoğunluk fonksiyonu,
$$ p(x)~dx = \lambda~e^{-\lambda x}~dx $$
ile verilir.
Burada $\lambda$ paramatresi Poisson dağılımındaki oran parametresi ($\lambda#)'ya benzer. Örneğin, CCD dedektörün bir pikseline saniyede gelen foton sayısı, Kepler Uzay Teleskobu'nun baktığı alandaki Güneş-benzeri yıldızların sayısı, bir saat içinde bir mağazaya gelen müşteri sayısı, yıllık deprem sayısı gibi nicelikler, Poisson dağılımının parametresi olan oran parametresine karşılık gelir.
Ancak, günlük hayatta olaylar arasındaki geçen süre söz konusu olduğunda, oran yerine zaman birimiyle konuşmyı yeğleriz. Örneğin, 10000 ADU sayım üretecek kadar fotonun algılanması için 20 saniye poz süresinin gerekmesi (20 / 10000 = 0.002 bir oran yerine), bir makinenin hata vermeden çalışacağı yıl sayısının 10 yıl olması (0.1 başarısızlık / yıl demek yerine) gibi. Örnekler çoğaltılabilir, bu nedenle üstel dağılımın “ortalama” terminolojisi Poisson süreçlerdeki oran parametresi $\lambda$ yerine, bunun tersi olan zaman parametresi $1 / \lambda$ üzerine kuruludur. Örneğimizden hareketle 10000 ADU 20 saniyede algılanıyorsa, "Poisson oranı" 0.002'dir, 10000 ADU algılamak için 20 saniye poz süresi vermek gerekir. Dolayısı ile $X \sim exp(0.002)$ verildiğinde kastedilen bu Poisson oranıdır.
Üstel dağılım, bir Poisson olayının gerçekleştiği iki durum arasındaki süre üzeirnden tanımlandığından bu iki olay arasındaki bekleme süresinde herhangi bir olay gerçekleşmez. Bu durum aşağıdaki şekilde ifade edilebilir.
$$ P(x = k) = \frac{\lambda^k~e^{-\lambda}}{k!} $$ $$ P(x = 0) = \frac{\lambda^0~e^{-\lambda}}{0!} = e^{-\lambda} $$
Tek bir zaman biriminde değil de bir $t$ zaman aralığında hiçbir olayın gerçekleşmeme olasılığı bilinmek isteniyorsa P(x=0, 1/t aralıkta) P(x=0, 1/t aralıkta) ... * P(x=0, 1/t aralıkta) çarpımı t kez yapılmalıdır. Bu durumda
$$ P(x = 0, t) = e^{-\lambda~t} $$
elde edilir. Bu şekilde elde ettiğimiz olasılık "bir olay gözlendikten sonra bir sonraki olaya kadar olan bekleme süresinin t zaman biriminden büyük olma olasılığı" olup $P(t > T)$ birikimli olasılığıdır. Birikimli olasılık $P(t < T$) ile ifade olunduğundan, bu nicelik $1 - CDF$'e eşittir.
$$ P(t \le T) = 1 - P(t \gt T) = CDF = 1 - e^{-\lambda~t} $$
Birikimli yoğunluk fonksiyonunun türevi olasılık yoğunluk fonksiyonunu verecektir.
$$ PDF = \frac{d~CDF}{dt} = \frac{d(1 - e^{-\lambda~t})}{dt} \rightarrow p(t)~dt = \lambda~e^{-\lambda~t}~dt $$
olasılık yğunluk fonksiyonu elde edilmiş olur.
Sonuç olarak, birim zaman başına bir Poisson sürecine uyan olay sayısı bir Poisson dağılımını izliyorsa, olaylar arasındaki süre üstel dağılımı takip eder.
Kuvvet Yasası Dağılımları¶
Kuvvet yasası (zaman zaman "ölçekleme yasası" olarak da adlandırılır), bir nicelikteki göreli bir değişikliğin, ona kuvvet yasasıyla bağlı diğer bir nicelikte orantılı bir değişme yol açtığını belirtir. Yasanın en basit örneği örneklerinden biri ters-kare yasasıdır. İki cisim arasındaki mesafe arttıkça aralarındaki çekim kuvveti mesafenin karesiyle azalır ya da iki özdeş yıldızdan yakın olanın parlaklığı, uzak olanın parlaklığından aralarındaki mesafenin karesi kadar fazladır. Yasaya uyan bir dağılımın olasılık yoğunluk fonksiyonu,
$$ p(x)~dx = a~x^{a - 1} $$
ile verilir.
Örneğin $1 M_{\odot}$ ile $10 M_{\odot}$ arasındaki anakol yıldızlarının kütlelerinin olasılık dağılımı $f(M)~\alpha~M^{-3.5}$ ile bağıntılıdır, zira bir yıldızın anakol yaşı kütlesinin -3.5. kuvvetiyle orantılıdır. Yaşamakta olduğumuz CoVid19 hastalığına sebep olan SARS-CoV-2 virüsünün başlangıçtaki yayılma hızının üstel olduğu, sonrasında kuvvet yasasına uyduğu ve yine üstel olarak azaldığı tespit edilmiştir.
Binom Dağılımı¶
Başarı ya da başarısızlık, kazanç ya da kayıp, "yazı" ya da "tura" gibi sadece iki sonucun mümkün olduğu ve tüm denemeler için başarı (p) ve başarısızlık olasılığının (1-p) toplamının 1 olduğu süreksizdağılıma Binom Dağılımı denir. Deneyin tekrar sayısının 1 olması durumunda elde edilen özel dağılımı Bernoulli Dağılımı adı verilir. Bununla birlikte, sonuçların eşit derecede olası olması gerekmez ve her deneme birbirinden bağımsızdır. Binom dağılımının parametreleri n ve p'dir; burada n, toplam deneme sayısıdır; p, her denemede başarı olasılığıdır. Olasılık dağılım işlevi şu şekilde verilir:
$$ f(k, n, p) = Pr (k, n, p) = {n \choose k} p^k (1 - p)^{n-k} $$
Burada $n \choose k$, n'in k'lı kombinasyonlarının sayısıdır ve
$$ {n \choose k} = \frac{n!}{k!~(n-k)!} $$
ile verilir.
Bernoulli Dağılımı¶
Bernoulli dağılımı, binom dağılımının tek bir denemenin yapıldığı (n = 1) özel bir hali olup, olasılık kütle fonksiyonu (ing. Probability Mass Function, PMF) şu şekilde verilir:
$$ f(k, p) = p^k q^{(1-k)} = p^k (1-p)^{(1-k)} $$
Python'da Dağılımlar¶
Python'da tüm olasılık dağılım fonksiyonlarına ilişkin fonksiyonların yer aldığı modüller bulunmaktadır. Bu modüllerden en sık kullanılanları scipy.stats ve np.random paketleridir. Bu modüllerdeki fonksiyonlarla ilgilenilen dağılımın yapısına uygun örnekler ya da dağılım fonksiyonları oluşturularak incelenebilir. Çizim için matplotlib.pyplot fonksiyonları hist, bar ve plot en sık kullanılanlardır. Ayrıca seaborn modülü fonksiyonları da kolay kullanılabilmeleri nedeniyle önemli bir alternatif oluşturmaktadır.
scipy.stats paketindeki neredeyse tüm dağılım fonksiyonlarında bulunan bazı argüman ve parametreler, Normal Dağılım örneğinde (stats.norm) aşağıda listelenmeye çalışılmıştır.
rvs(loc=0, scale=1, size=1, random_state=None): Söz konusu dağılımdan bir örnek oluşturmak için kullanılır.
pdf(x, loc=0, scale=1): Olasılık yoğunluk fonksiyonu oluşturmak için kullanılır.
logpdf(x, loc=0, scale=1): Olasılk yoğunluk fonksiyonun logaritması (doğal logaritma, ln)'nı verir.
cdf(x, loc=0, scale=1): İlgilenilen dağılım türü temelinde bir birikimli (kümülatif) yoğunluk fonksiyonu oluşturmak için kullanılır.
logcdf(x, loc=0, scale=1): Birikimli yoğunluk fonksiyonunun logaritmasını (ln) verir.
sf(x, loc=0, scale=1): Kümülatif fonksiyon $P(X \le x)$ olasılığını verirken $sf$ (ing. survival function) $P(X \gt x) = 1 - (P(X \le x)$ olasılığını verir.
logsf(x, loc=0, scale=1): sf fonksiyonunun logaritmasıdır.
ppf(q, loc=0, scale=1): (ing. percent point function), birikimli yoğunluk fonksiyonunda elde edilen olasılığın karşılık geldiği $X = x$ değeridir. Bu anlamda birikimli yoğunluk fonksiyonunun tersidir. Birikimli yoğunluk fonksiyonu $X = x$ değeri argüman olarak sağlandığında $P(X \ge x)$ olasılığını verirken, $ppf$ verilen $P(X \ge x)$ olasılık değeri için $X = x$ değerini verir.
isf(q, loc=0, scale=1): (ing. inverse survival function), ppf'e benzer şekilde $1 - cdf$ olasılığının karşılığı olan $X = x$ değerini verir.
moment(n, loc=0, scale=1): Dağılımın n. momentini verir.
stats(loc=0, scale=1, moments=’mv’): Dağılıma ilişkin momentleri verir. İstenen momentler
momentsparametresine baş harfiyle sağlanır. (Ortalama: Mean(‘m’), Varyans: variance(‘v’), Çarpıklık: skew(‘s’), ve/veya kurtosis(‘k’)).entropy(loc=0, scale=1): Dağılımın entropisini verir.
fit(data, loc=0, scale=1):
dataparametresi üzerinden verilen bir veriye ilgili dağılımı uyumlar (fit eder) ve bu uyumalamanın dağılımın türüne göre ilgili parametrelerini döndürür.expect(func, args=(), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds): Tek değişkenli bir fonksiyonun bir dağılıma göre beklenen değerini döndüren fonksiyondur.
median(loc=0, scale=1): Dağılımın ortancasını (medyan değerini) döndürür.
mean(loc=0, scale=1): Dağılımın ortalama değerini döndürür.
var(loc=0, scale=1): Dağılımın varyansını ($\sigma^2$) döndürür.
std(loc=0, scale=1): Dağılımın standart sapmasını ($\sigma$) döndürür.
interval(alpha, loc=0, scale=1): Tüm dağılımın $\alpha$ yüzdelik bölümüne denk gelen aralığının uç noktalarını verir.
Python'da Dağılım Fonksiyonlarına Örnekler¶
Python'da herhangi bir dağılım fonksiyonuna uyan örnekler oluşturmak için genellikle numpy.random paketi fonksiyonları, bu fonksiyonları kullanarak analizler ve hipotez testleri yapmak için ise scipy.stats paketi fonksiyonları kullanılır. Başka pek çok dış paket de mevcut olmakla birlikte temel analizler için bu modüller yeterlidir. Şimdi yukarıda ortak elemanları verilen bu fonksiyonlara bazı örnekler verelim.
Tekdüze Dağılım Örneği¶
numpy.random fonksiyonlarından uniform ile tekdüze bir dağılımdan giderek daha fazla sayıda örnek seçilerek tekdüze bir dağılımın nasıl oluştuğunu örnekleyen bir dizi kod aşağıda verilmiştir. Ayrıca tekdüze dağılımın her aşamada bir normal fonksiyona göre ne kadar daha çarpık (ing. skewness) ve basık olduğuna ilişkin değerler (dağılımın $skew$ ve $kurtosis$ parametreleriyle) de verilmiştir. Tekdüze bir dağılımın simetrik olması beklendiği için örnek sayısı arttırıldıkça çarpıklığının 0'a, ortancaya (medyan), basıklığı gösteren Kurtosis değerinin ise -1.2'ye yaklaşması beklenmelidir.
import numpy as np
from matplotlib import pyplot as plt
from scipy import stats as st
%matplotlib inline
x = np.random.uniform(-5,5,10)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=20)
plt.show()
x = np.random.uniform(-5,5,100)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=20)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
print('Normal bir dagilima gore basiklik: ', st.kurtosis(x))
print('Normal bir dagilima gore carpiklik: ',st.skew(x))
x = np.random.uniform(-5,5,1000)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=20)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
print('Normal bir dagilima gore basiklik: ', st.kurtosis(x))
print('Normal bir dagilima gore carpiklik: ',st.skew(x))
N = 100000
x = np.random.uniform(-5,5,N)
print(x.mean(), x.std(), np.median(x))
plt.xlabel('x')
plt.ylabel('PDF')
plt.hist(x, bins=20, density=True)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
Standart Normal Dağılım Örneği¶
stats.norm.cdf fonksiyonu $x$'in $-\infty$'dan $\Delta x = +1.65\sigma$ 'ya kadarki aralıkta olma olasılığını verir (standart normal dağılımın bu değere kadar sol tarafının altında kalan alan).
import scipy.stats as st
z = 1.645
pz = st.norm.cdf(z)
print("Standart normal dagilimda x = +{:.2f} sigma sapma olasiligi : {:.2f}". \
format(z,pz))
Eğer x'in $\Delta x = \pm1.65\sigma$ aralığında olma olasılığı isteniyorsa bu kez bu değerlerin fonksiyonlara ayrı ayrı sağlanması ve aradaki farkın alınması gerekir.
z1 = -1.645
z2 = 1.645
pz_z1_z2 = st.norm.cdf(z2) - st.norm.cdf(z1)
print("Standart normal dagilimda x'in = +/-{:.2f} sigma dahilinde olma olasiligi : {:.2f}". \
format(z2,pz_z1_z2))
Dolayısıyla normal dağılımı temsil eden eğrinin en yüksek eğime sahip tanjantlarının eğriyi ve x eksenini kestiği noktalar; sırasıyla 1 ve 2 $\sigma$ verilerek öğrenilebilir.
z1 = -1.0
z2 = 1.0
pz_z1_z2 = st.norm.cdf(z2) - st.norm.cdf(z1)
print("Standart normal dagilimda x'in = +/-{:.2f} sigma dahilinde olma olasiligi : {:.2f}". \
format(z2,pz_z1_z2))
z1 = -2.0
z2 = 2.0
pz_z1_z2 = st.norm.cdf(z2) - st.norm.cdf(z1)
print("Standart normal dagilimda x'in = +/-{:.2f} sigma dahilinde olma olasiligi : {:.2f}". \
format(z2,pz_z1_z2))
Tersten giderek verilen bir olasılığın merkezden kaç standart sapma uzaklığa denk geldiği (yani $z$ değeri) de bulunabilir. Bunun için ppf (percent point function) fonkiyonu kullanılır. Fonksiyonun yine $+\Delta x$'e kadar olasılık verdiği (eğrinin sol tarafının altında kalan) akılda tutulmalıdır.
pz = 0.95
z = st.norm.ppf(pz)
print("Standart normal dagilimda {:.2f} olasiligin karsilik geldigi sapma degeri : {:.2f}". \
format(pz,z))
İki taraflı olasılık isteniyorsa, $px$ argümanı buna uygun olarak verilmelidir.
pz = 0.05
z = st.norm.ppf(pz)
print("Standart normal dagilimda +\-{:.2f} olasiligin karsilik geldigi sapma degeri : +/-{:.2f}". \
format(pz,z))
Aynı amaçla daha kullanışlı olan interval fonksiyonu da kullanılabilir. norm.interval fonksiyonu argüman verilen bir olasılığın $x$'in hangi z = $\pm \Delta x$ değerleri arasında olma olasılığı olduğunu verir. Bir başka deyişle, verilen olasılığa karşılık gelen alanın uç noktaları olan $z$ değerlerini verir.
pz = 0.90
z = st.norm.interval(pz)
print(z)
print("Standart normal dagilimda +\-{:.2f} olasiligin karsilik geldigi sapma degeri : +/-{:.2f}". \
format(pz,z[1]))
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Polygon
fig, ax = plt.subplots()
x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 1000)
ax.plot(x, st.norm.pdf(x), 'r-', lw=2, alpha=0.6, label='norm pdf')
xi = st.norm.ppf(0.001)
xf = st.norm.ppf(0.05)
ix = np.linspace(xi,xf)
iy = st.norm.pdf(ix)
verts = [(xi, 0), *zip(ix, iy), (xf, 0)]
poly = Polygon(verts, facecolor='b')
ax.add_patch(poly)
ax.set_xlabel("z")
ax.set_ylabel("pdf")
plt.arrow(-2.75, 0.10, 0.5, -0.05, capstyle="round",
head_width=0.02, head_length=0.05)
plt.text(-2.6,0.115,"p(z)")
plt.show()
Bu alanın değerini birikmli (kümülatif) dağılım fonksiyonu cdf üzerinden görebiliriz. Bir başka deyişle, kümülatif dağılım fonksiyonu, normal dağılım fonksiyonunun altında kalan alanın z değeri ile değişimini verir. $z \rightarrow \infty$ iken $cdf \rightarrow 1.0$ 'dır.
from matplotlib.patches import Polygon
fig, ax = plt.subplots()
x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 1000)
plt.plot(x, st.norm.cdf(x), 'r-', lw=2, alpha=0.6, label='norm pdf')
plt.axvline(x=st.norm.ppf(0.05))
plt.axhline(y=st.norm.cdf(st.norm.ppf(0.05)))
x_degeri = "z = {:.3f}".format(st.norm.ppf(0.05))
y_degeri = "p(z) = {:.3f}".format(st.norm.cdf(st.norm.ppf(0.05)))
plt.text(1, 0.1, y_degeri)
plt.text(-2.8,0.15, x_degeri)
plt.xlabel('z')
plt.ylabel('cdf')
plt.show()
Poisson Dağılımı Örneği¶
Diyelim ki bir yıldız gözleminde CCD'nin herhangi bir pikseline belirli bir dalgaboyu aralığında yıldızdan saniyede ortalam 5 foton geliyor olsun. Dedektörün bu dalgaboyu aralığındaki kuantum etkinliği de %80 olsun. Yani gelen 5 fotonun 4'ü dedektör tarafından algılansın ($\mu = 4$). Gece boyunca alınan 1'er saniyelik poz süreleriyle alınan 1000 görüntüde algılanan foton sayılarının değişimi bir Poisson dağılımı sergiler.
from scipy.stats import poisson
ornek_poisson = poisson.rvs(mu=4, size=1000)
fig,ax = plt.subplots()
ax.hist(ornek_poisson)
ax.set(xlabel='Poisson Dagilimi', ylabel='Frekans')
plt.show()
Foton sayıları için gözlenme sayıları (frekans) yerine, normalize yoğunluk (o aralıktaki ölçüm sayısının toplama bölümü) verilmek istenirse pyplot.hist fonksiyonunda density parametresi $True$ değerine ayarlanır. Poisson dağılımının bir süreksiz fonksiyon olduğunu düşünerek olasılık yoğunluk fonksiyonu (PDF) yerine, istenen her değer için Poisson olasılık kütle fonksiyonunu (PMF) örnek dağılımın üzerine çizdirmek bu şekilde mümkündür. Burada PMF, fonksiyon süreksiz olduğunu vurgulamak üzere, özellikle az sayıda (10) nokta için, 'o' sembolüyle verilmiştir.
from scipy.stats import poisson
ort = 4. # oran parametresi
ornek_poisson = poisson.rvs(mu=ort, size=1000)
fig,ax = plt.subplots()
ax.hist(ornek_poisson, density = True)
# %1 ile %99 olasiliklar arasina karsilik gelen x degerleri
# arasindaki 10 nokta icin olasilik kutle fonksiyonu (pmf)
x = np.linspace(poisson.ppf(0.01, mu=ort), poisson.ppf(0.99, mu=ort),10)
pmf = poisson.pmf(x, mu=ort)
ax.plot(x, pmf, 'ro')
ax.set(xlabel='Poisson Dagilimi', ylabel='PMF')
plt.show()
Diyelim ki poz süresini uzatıp 5 saniyeye çıkarıyoruz. Ancak birim zamanda algılanan foton sayısı da değişmiyor. Buna karşın artık $\mu = 5~x~4 = 20$ foton algılarız. Doğal olarak gece boyunca alacağımız görüntü sayısı düşecektir. Poz süresini 5 kat arttırdığımız için toplam görüntü sayısı da 5 kat azalmış ve 200'e düşmüş olsun.
from scipy.stats import poisson
ornek_poisson = poisson.rvs(mu=20, size=200)
fig,ax = plt.subplots()
ax.hist(ornek_poisson)
ax.set(xlabel='Poisson Dagilimi', ylabel='Frekans')
plt.show()
Görüldüğü gibi dağılım normal dağılıma daha fazla yaklaşmıştır. Astronomide gözlemlerinde zaman çözünürlüğünden ve toplam görüntü sayısından ödün verilebilecek durumlarda poz süresinin arttırılmasıyla Poisson dağılımı gösteren foton istatistiğinin Gaussyen dağılımı yakınsaması amaçlanır. Zira görüntü işleme rutinleri, sayımların Normal dağıldığı varsayımı üzerine kuruludur.
Merkezi Limit Teoremi¶
Olasılık teorisinde, Merkezi Limit Teoremi (ing. Central Limit Theorem, CLT), herhangi bir rastgele değişenin kendisi normal dağılıyor olması dahi, o rastgele değişenin olası tüm değerlerinden oluşturulan çok sayıdıa örnekleminin ortalamalarının (bu örneklemlerin kendisi de normal dağılmak zorunda değildir) normal dağılacağını söyler. Teorem olasılık teorisinde büyük önem taşır. Çünkü teoremden normal dağılımlar için geçerli olasılıksal ve istatistiksel yöntemlerin diğer dağılım türlerini de içeren birçok probleme uygulanabileceği sonucu çıkarılabilir. Diyelim ki seçilmiş 1000 Sefeid türü değişenin zonklama dönemleri normal dağılmıyor olabilir. Ancak çok sayıda (örneğin 100'er tane) Sefeid içeren çok sayıda örneklemin (örneğin 100 küresel kümenin) ortalama zonklama dönemlerinin dağılımına bakıldığında normal dağılım gösterdikleri görülür. Normal dağılım ve örnek için bkz.
Bir örnekle merkezi limit teoremini görselleştirerek anlamaya çalışalım. Pek çok sınıfa aynı anda uygulanan bir sınav düşünelim. Sınavda alınan notların sınıflar içinde ve sınıflarının ortalamasının da genelde nasıl dağıldığına bakalım. Öncelikle 50 öğrenciden oluşan bir sınıf için numpy.randint() fonksiyonuyla rastgele belirlenmiş sınav notlarından oluşan bir örneklem oluşturalım.
import numpy as np
# rastgele sinav notunu ayni degerleri verecek sekilde
# ayarlayalim (seed edelim)
np.random.seed(1)
# notlarin tamsayi oldugunu varsayarak 50 rastgele
# notu 40 ile 70 arasinda olusturalim.
notlar = np.random.randint(40, 70, 50)
print(notlar)
print('Sinif ortalamasi: {:.2f}'.format(notlar.mean()))
print('Standart sapmasi: {:.2f}'.format(notlar.std()))
Sınıf içindeki not dağılımını yukarıda gördük; bir de histogramına bakalım.
from matplotlib import pyplot as plt
# sinif ici not dagilimi
plt.hist(notlar)
plt.show()
Hiç de normal dağılıyormuş gibi görünmüyor. Nitekim numpy.random.randint() sınırları verilen (örneğimizde 40 ve 70) bir tekdüze dağılımdan istenen büyüklükte (örneğimizde 50) bir rastgele örneklem oluşturur. Şimdi örneklem sayısını 1 sınıftan 1000 sınıfa çıkaralım ve sınıfların ortalamasının nasıl dağıldığna bakalım.
import matplotlib.pyplot as plt
# yine ayni sekilde istikrarli bir orneklem icin
# rastgele tam sayi ureticisini seed edelim.
np.random.seed(1)
# Yine 40 ile 70 arasindaki notlar icin
# bu kez 1000 siniftan olusan orneklemimizdeki
# her bir sinifin ortalamasini alip sinif_ortalamalari
# degiskenine alalim
sinif_ortalamalari = [np.mean(np.random.randint(40, 70, 50)) for i in range(1000)]
# dagilimi bir histogramla gosterelim.
plt.hist(sinif_ortalamalari)
plt.show()
print('Sinif ortalamlarinin ortalamasi {:.2f}'.format(np.mean(sinif_ortalamalari)))
print('Sinif ortalamlarinin standart sapmasi {:.2f}'.format(np.std(sinif_ortalamalari)))
Her ne kadar tek bir sınıf içindeki not dağılımı normal değil tekdüze (uniform) olsa da 1000 sınıfın ortalamaları normal dağılmaktadır. Bir sınıf içinde 8.65 olan standart sapma değerinin de 1000 örenklemin ortalamaları söz konusu olduğunda 1.21'e kadar düştüğünü görüyoruz.
Örneklem Dağılımı (ing. Sampling Distribution): Bu şekilde normal dağılmıyor olsa dahi örneklemlerin ortalamalarının dağılımı örneklem dağılımı adını alır.
Ortalamanın Beklenen Değeri (ing. Expected Value of the Mean): Örneklem dağılımının ortalama değeri popülasyonun ortalaması için beklenen değerdir. Bir başka deyişle, bir popülasyondan rastgele seçilen örneklemlerin ortalamalarının dağılımının ortalaması o popülasyonun ortalamasının beklenen değeridir.
Standart Hata (ing. Standard Error): Örneklem dağılımının ortalama değerinin standart sapmasıdır. Dolayısı ile popülasyon ortalamasının standart sapmasının da beklenen değeridir. Tıpkı standart sapmanın herhangi bir ölçümün ortalamadan sapmasının bir ölçütü (tüm sapmaların ortalaması) olduğu gibi, standart hata da bireysel örneklemlerin ortalamalarının, örnek dağılımının ortalamasından sapmasının ölçütüdür. Bu nedenle de standart hatadır çünkü ortalamaların dağıımının (örneklem dağılımının) ortalamasının popülasyonun ortalamasını ne kadar temsil ettiğine olan güvenimizin bir ölçütüdür.
Örneklemlerin ortalamasının bir dağılımı olduğu gibi diğer istatistiklerin de (standart sapma, korelasyon katsayısı, $\chi^2$, ...) birer dağılımı vardır. Aynı şekilde ortalamaların ortalamasının bir standart hatası olduğu gibi diğer istatistiklerin de standart hataları vardır.
Hipotez Testi¶
Peki Merkezi Limit Teoremi bize nasıl yardımcı olacak? Teorem, ele alınan bir örneğin karşılaştırma yapılan bir popülasyonu temsil edip etmediği hipotezini test etmemizi sağlar. Bunun için örneğin ortalama değerini belirleyip, bu örneğin karşılaştırma yapılan popülasyondan gelip gelmediğini tahmin etmek için örneklem dağılımı ile karşılaştırabiliriz. Hipotez testinin belirli kuralları olup, öncelikle testin temellerini anlamaya ihtiyaç vardır.
Herhangi bir bilimsel problemin çözümünde öncelikle bir hipotez geliştirilir. Hipotez, bir problemin çözümüne elde az sayıda kanıt ve çalışma varken yapılan ve daha ileri araştırma için bir başlangıç noktası oluşturan öneridir.
Sıfır hipotezi ($H_0$), bir öneriyi test etmek üzere, gözlenen bir olayın sadece rastgele süreçler ile oluştuğunu, örnekleme ya da deneysel hatalardan kaynaklandığını öne süren hipotezdir. Gözleme ilişkin alternatif hipotezin ($H_1$) yani gözlenen olgunun rastgele süreçlerle oluşmadığının ve olgunun parametreleri arasında bir korelasyon bulunduğu fikrinin geçerliliğinin test edilmesi için kullanılır.
Sıfır hipotezinin yanlışlanamaması durumunda gözlenen olayın sadece rastgele süreçlerden kaynaklandığı (ya da sıfır hipotezinin geçerli olduğu) kabul edilir. Sıfır hipotezinin yanlışlanması durumunda ise alternatif hipotezin doğruluğu kanıtlanmış olmaz. Ancak gözlenen olgunun sadece rastgele süreçlerden (ya da sıfır hipotezinin öngördüğü süreçlerden) kaynaklanmadığı sonucuna varılır.
Sıfır hipotezinin yanlışlanması için, alternatif hipotezin gerçekleşebilme olasılığının, rastgele süreçler ile gerçekleşme ya da sıfır hipotezinin ilgilendiği dağılıma göre benzer bir olayın gözlenebilme olalsılığından anlamlı bir mertebede yüksek olması gerekir. Örneğin; bir bitkiyi sulamanın, bitkinin büyüme hızıyla ilişkili olduğunu söyleyen bir alternatif hipotezin kabul edilebilmesi için, bitkinin büyüme hızıyla sulama arasında bir ilişki olmadığını söyleyen sıfır hipotezinin yanlışlanması gerekir. Daha radikal bir örnek olarak; piramitleri uzaylıların inşa ettiğini iddia eden alternatif hipotezin kabul edilebilmesi için öncelikle, piramitleri uzaylıların inşa etmediğini iddia eden sıfır hipotezinin yanlışlanması gerekir. Tabi ki bu durumda piramitleri uzaylılıarın inşa ettiğine ilişkin doğrudan bir kanıtın da bulunması icap eder.
Örnek: CoVid 19 aşısı¶
Diyelim ki CoVid-19 hastalığı için bir aşı önerisi var. Aşının başarısını test etmek için $n$ bireyi rastgele seçerek bir örneklem oluşturuyor ve bu bireyleri aşıyla enjekte ediyoruz. Diyelim ki bu $n$ bireyin $k$ tanesinin bir süre sonra (örneğin 1 yıl) yapılan testleri (ve diğer klinikleri) negatif sonuç veriyor ve bu $k$ bireyi "sağlıklı" olarak değerlendiriyoruz. Şimdi aşının işe yarayıp yaramadığına yönelik hipotez testinin nasıl yapılacağına bakalım.
$H_0$ (Sıfır hipotezi): Aşının hiçbir etkisi yoktur ve sağlıklı olan $k$ birey hiçbir şey yapılmasa da hastalığı belirlenen sürede (1 yıl) kapmayacak bireylerdir.
$H_1$ (Alternatif hipotez): Aşı bireylerin hastalıktan etkilenmesini engellemektedir.
Test edilmesi gereken hipotez $H_0$ hipotezidir. Eğer bu hipotez doğruysa, aşılanan $n$ bireydeki $k$ negatif sayısının (örneklem büyüklüğüne oranının), tüm toplumdaki $K$ negatif sayısına (popülasyona oranına) benzer olması gerekir.
Diyelim ki hastalığı kapma olasılığı aşı vs. gibi bir engelleyici olmadan %25 olsun (p = 0.25). Bu durumda da kapmama olasılığı (q = 1 - p = 0.75) %75 olur ve bu deney iki olasılığın olduğu bir Bernoulli deneyidir ve dağılım da Bernoulli Dağılımı'na uyar.
Diyelim ki örneklemimiz 10000 hasta ile sınırlı (n = 10000). Bu durumda sağlıklı birey $k$ sayısı 1 yıl sonra 10000 x 0.75 = 7500 çıkarsa (Bernoulli dağılımının ortalama değeri) bunun tamamen rastgele olduğunu, aşının bir etkisi olmadığını söyleyebiliriz. Zira hiçbir şey yapılmasa da hastalığı kapmama olasılığı %75'tir. Ancak bu durumda hangi sayıda sağlıklı birey bulursak $H_0$ hipotezini reddetmeliyiz (aşının etkisi olmadığı sonucuna varmalyızı) gibi bir seçimle karşılaşırız. Zira 7500 değil de 7440 sağlıklı birey (%74.4) ya da 7590 (%75.9) sağlıklı birey bulmamız sadece örneklemimizin küçüklüğünden kaynaklanıyor olabilir mi?
Burada 7500'ün her iki tarafından belirli bir standart sapma ($\sigma$) değerini belirleyip, eğer deneydeki sağlıklı birey sayısı ortalamanın (7500) örneğin $2~\sigma$ ya da $2.5~\sigma$ içinde kalıyorsa $H_0$'ı kabul etmek aksi takdirde reddetmek iyi bir fikir olabilir. İşte burada merkezi limit teoreminin temel yaklaşımını hatırlayabiliriz. Her ne kadar ilgilendiğimiz dağılım bir Bernoulli dağılımıysa da yeterince büyük örneklemlerle çalışıldığı vakit bu dağılımın normal dağılıma yakınsayacağını biliyoruz. Dolayısı ile normal dağılımın parametreleri ve yaklaşımıyla hareket edebiliriz!
Bir Bernoulli dağılımının standart sapması $\sigma = \sqrt{n~p~q}$ ile verilir. Bu durumda örneğimiz için standart sapma değeri $\sigma = \sqrt{10000~0.25~0.75} = 43.30$ olarak bulunur. Burada iki konuya değinmek gerekir:
1) Çalışmamızın daha "güvenilir" olması için ortalamadan kaç standart sapma ($\sigma$) uzaklığı $H_0$'ı red ya da kabul etmek için kriter yapmalıyız?
2) "Yeterince büyük örneklem" nedir, bir ölçüsü var mıdır?
Önceliklle ikinci sorunun cevabının popülasyonun büyüklüğüne bağlı olduğunu, örneklemin bazı problemlerde sadece niceliğinin (büyüklüğünün) değil de aynı zamanda niteliğinin (içerdiği örneklerin popülasyon niteliğini ne kadar yansıtıyor olduğunun) de dikkate alınması gerektiğini söylemek gerekir.
Birinci sorunun cevabı ise "Güven Aralığı" (ing. Confidence Interval) ya da "Güven Düzeyi" (ing. Confidence Level) kavramlarıyla verildiğini belirterek bu kavramları açmak gerekir.
Güven Aralığı ve Güven Düzeyi¶
Güven aralığı, bir örnek dağılımın ortalama değerinin, ana dağılımın ortalama değerinden hangi olasılık mertebesi dahilinde uzak olduğunu gösteren aralıktır. Doğa bilimlerinde belirsizlik ölçütü olarak genellikle standart sapmayı kullanmak yeterlidir. Ekonomi ve siyaset bilimi gibi bazı başka disiplinlerde standart sapmayı vermek yerine güven aralığı (ing. confidence interval) ya da genellikle yüzde biriminde ifade edilen güven düzeyi (ing. confidence level) vermek tercih edilir.
Örneğin bir parametre için yapılan bir dizi ölçümün sonucu $101.25 \pm 0.01$ olarak verilmiş olsun. Aksi söylenmediği (örneğin bu hatanın ölçüm aletinin limiti olduğu gibi) sürece bir doğa bilimci bu belirsilziiği bir standart sapma olarak yorumlar ($\sigma = 0.01$). Ancak 2 standart sapma ($2~\sigma$) kullanarak ölçümün $\% 95$ güven düzeyinde $101.25$ sonucunu verdiğini söylemek de aynı içeriği taşır. Zira normal olduğu kabul edilen bir dağılımda %95 güven düzeyi de yaklaşık olarak $\pm 2\sigma$’ya (daha yaklaşık bir değer olarak $1.96~\sigma$) karşılık gelir. Benzer şekilde ölçüm sonucunun $(101.23 , 101.27)$ %95’lik güven aralığında olduğu da söylenebilir.
%95’lik güven aralığı, deney/gözlem tekrarlanmaya devam edilirse tekrarların %95’inde, gözlenen değerlerin bulunan güven aralığında çıkacağı anlamına gelir.
Örnek: CoVid-19 aşısı¶
CoVid-19 aşı önerisi örneğinde standart sapmayı Bernoulli Dağılımı'ndan hareketle belirlemiştik. Aşının hiçbir etki yapmadığı, aşı uygulamasının 1. yılı sonunda sağlıklı bireylerin tamamen rastlantısal olarak sağlıklı kaldıklarını ifade eden $H_0$ sıfır hipotezini reddetmek için $\%95$'lik bir güven düzeyi belirlemiş olalım. Bu durumda
$$ \mu \pm 2~\sigma = 7500 \pm 2 x 43.30 = 7500 \pm 86.60 $$
dahilinde bir güven aralığı oluşturabiliriz. 1 yıl sonunda sağlıklı birey sayısını 7587'den büyük bulursak $H_0$'ı reddedebiliriz. Yani "aşının hiçbir etkisi yoktur" önerisini reddetmiş oluruz! Bu alternatif hipotezi kabul ettiğimiz (yani aşının etkili olduğu) anlamına gelmez ama onu destekleyici bir kanıt oluşturur.
Örnek: Gaia Paralaksları¶
Gaia uzay teleskobu gökcisimlerinin hassas konum ölçümlerini yapmak üzere uzaya gönderilmiştir. Teleskop, galaksimiz içindeki yıldızların da paralaks ölçümlerini defalarca yaparak uzaklıklarının hassas bir şekilde belirlenmesini sağlamaktadır.
Bir yıldızın paralaks ölçümü 1000 defa yapılıyor ve uzaklğı parsek cinsinden hesaplandığında elde edilen ortalama uzaklık değeri ($\mu_d$) ve belirsizliği ($\sigma_d$) $d = 210 \pm 5$ pc olarak veriliyor. %95 güven seviyesine karşılık gelen aralığı bulunuz. Problemi hem temel istatistik bilgimizden hareketle, hem de genelleştirmek üzere Python dili ve ekstra paketlerin olanaklarıyla çözelim.
Çözüm: %95’lik güvenilirlik seviyesine karşılık gelen $z$ değeri için standart normal dağılım tablolarına ya da grafiklerine bakılır. Bu değer scipy.stats.norm.interval fonksiyonuyla doğrudan elde edilebilir ya da $p = 0.95$ için aşağıdaki formülden yararlanılarak hesaplanan $z$ değeri tablodan alınabilir. Daha sonra,
$$ z = \frac{|x - \mu|}{\sigma} $$
ifadesinden hareketle $x$ değişkeninin %95 olasılıkla içinde bulunduğu güven aralığı tesis edilir.
from scipy import stats as st
p = 0.95
z = st.norm.interval(p)
print("%{:.2f} olasiliga karsilik gelen z-değeri: {:.2f} 'dir.".
format(p*100, z[1]))
Görüldüğü üzere $z$-değeri %95'lik güven aralığı için 1.96 olarak verilmektedir. Bu durumda güven aralığı aşağıdaki şekilde tesis edilebilir.
$$ z = \frac{|x - \mu|}{\sigma} \Rightarrow x = \mu \pm z~\sigma \Rightarrow x = 210 \pm 1.96~x~5.0 = 210 \pm 9.8~pc$$
Bir başka deyişle, 1000 kez paralaks ölçümü yapılan bu yıldızın uzaklığı %95 güven düzeyinde $(200.2 , 219.8)$ pc aralığındadır.
Aynı sonuca scipy.stats.norm.interval fonksiyonunu aşağıdaki şekilde kullanarak da ulaşılabilir.
from scipy import stats as st
p = 0.95
mu = 210.0
sigma = 5.0
aralik = st.norm.interval(p,loc=mu, scale=sigma)
print("Yildizin uzakligi icin %{:.2f}'lik guven araligi ({:.2f}, {:.2f}) pc'tir.".\
format(p*100, aralik[0], aralik[1]))
Görüldüğü üzere yeterince fazla ölçüm yapıldığında ölçüm sonuçlarının normal dağıldığını varsayarak, normal dağılımın parametreleriyle (ortalama ve standart sapma) bir güven aralığı oluşturulmuştur. Oysa ki tek bir ölçüm (ya da az sayıda ölçüm) söz konusu olduğunda ölçüm hataları normal (ya da Gaussyen) dağılmaz, Poisson Dağılımı'na (ya da genelleştirilmiş hali olan Drichlet Dağılımı'na) uyar. Ancak Merkezi Limit Teoremi gereğince yeterince büyük bir örneklem söz konusu olduğunda, rastgele bir değişkenin değerinin dağılımı Normal Dağılım'a yakınsar. Bu durumda da normal dağılım yaklaşımından hareketle analiz yapılabilir.
Diyelim ki bu yıldızın uzaklığı bir başka yöntemle (örneğin çift yıldız gözlemleriyle oluşturulan kalibrasyondan hareketle ya da HIPPARCOS uzay teleskobu ölçümleriyle) belirlenmiş olsun. Bu ölçümlerin hangi güvenilirlik düzeyinde eşit kabul edileceğine karar verildikten sonra sıfır hipotezi ($H_0$) ve alternatif hipotez ($H_1$) aşağıdaki şekilde oluşturulabilir.
$H_0$: Yildizin Gaia uzakligi cift yildiz gozlemlerinden belilrlenen kalibrasyonla uyumlu degildir.
$H_1$: Yildizin Gaia uzakligi cift yildiz gozlemlerinden belilrlenen kalibrasyonla uyumludur.
Bu tür bir karşılaştırmanın nasıl yapılacağını biraz sonraya bırakarak hipotez testlerinde yapılabilecek olası hataları görelim.
Hipotez Testinde Tip I ve Tip II Hataları¶
Tip-I Hata (ing. Type-I Error): Sıfır hipotezinin ($H_0$) doğru olmasına karşın reddedilmesiyle oluşan hatadır. Özellikle örnek sayısının azlığından kaynaklanabilir.
Tip-II Hata (ing. Type-II Error): Sıfır hipotezini ($H_0$) yanlışlığına rağmen kabul edilmesiyle oluşan hatadır.
İstatistiksel Anlamlılık Seviyesi¶
İstatistiksel anlamlılık, bir örneklemden türetilen bir istatistiğin o örneklemin alındığı popülasyondaki bir olguyu ne düzeyde (hangi olasılıkla) temsil ettiğinin ölçüsüdür. Örneğin "Bir popülasyondan alınan bir örneğin ortalamasının, o örneğin popülasyonu temsil etmediği ve şans eseri oluştuğu (örneklem hatası) kanısına bu iki ortalamada hangi düzeyde bir fark olursa varılmalıdır?" sorusunun cevabı araştırmacı tarafından belirlenen istatistiki anlamlılık seviyesiyle karşılaştırma yapılarak verilir. Eğer ulaşılan sonucun (örneğin iki ortalama arası farkın) olasılığı seçilen anlamlılık seviyesinden küçükse bu sonucun istatistiksel olarak anlamlı olduğu ve şans eseri (örneklem hatası nedeniyle) oluşma ihtimalinin istatistiksel anlamlılık seviyesiyle izin verilen daha küçük olduğu sonucu çıkar.
Gözlenen olayın rastgele süreçlerden meydana gelip gelmediği test edilmek istendiği için, olayın gerçekleşme ihtimali bir normal dağılım kullanılarak incelenir. Eğer olaya ilişkin dağılımın rastgele dağılım olmadığı biliniyorsa ilgili dağılım kullanlabilir. Ancak yeterli sayıda örnek dağılım elemanı bulunuyorsa merkezi limit teoremi gereğince normal dağılım varsayımı, bulunmuyorsa normal dağılım yerine, düşük serbestlik derecelerine daha duyarlı olan t-dağılımı varsayımı yapılabilir.
Istatistiksel anlamlılığın belirlenebilmesi için kullanılan kritik değer $\alpha$, bir olasılık değeri anlamındadır. Farklı bilim dallarında farklı değerler kabul görebilmektedir. Gözlenen olgunun rastgeleliğinin büyük oranda beklenir olduğu durumda daha küçük olasılık değerleri verilmesi uygun olur.
İstatistiki anlamlılık seviyesi (ing. statistical significance), toplam olasılık olan $1$’den güven seviyesinin farkıdır.
Sıfır hipotezi ($H_0$) doğru ise ve gözlemsel verinin ya da ilgili parametrenin normal dağılım sergilemesi bekleniyorsa, normal dağılımca beklenen en olası değer ortalama değerdir. Bu durumda kritik olasılık değeri ($\alpha$) belirlendiği takdirde, kritik değere göre ortalamaya daha yakın olan toplam olasılık değerleri sıfır hipotezinin reddedilememesi anlamına gelir. Bu durumda sıfır hipotezi $1- \alpha$ güven aralığında reddedilememiş olur.
[CoVid-19 aşısı](#Örnek:-CoVid-19-aşısı) örneğinde güven düzeyinin %95 olarak belirlenmesi durumunda istatistik anlamlılık seviyesi $\alpha = 1 - 0.95 = 0.05$ 'tir. 1 yıl süren aşılama sürecinin sonunda bu düzeyin üzerinde bir iyileşme (ki bunun 7587 hastaya denk geldiğini görmüştük) istatistiki olarak anlamlı bir sonuç bulunduğu değerlendirmesi yapılabilir. Bu da alternatif hipotez $H_1$ için istatistiki olarak anlamlı bir kanıta ulaşıldığı anlamına gelir. Ancak bu $H_1$'in kabul edildiği anlamına da gelmez.
[Gaia paralaksları](#Örnek:-Gaia-Paralaksları) örneğinde;
$H_0$: Gözlenmiş çift yıldızlar üzerinden belirlediğiniz bir korelasyon %95 güven düzeyinde Gaia ölçümüyle ($\mu = 210~pc$) uyumludur.
şeklinde ifade edilen bir sıfır hipotezini yanlışlamanız için kalibrasyonunuzun $\alpha = 1 - 0.95 = 0.05$ istatistiki anlamlılık seviyesinin karşılık geldiği $9.8~pc$ değerinin ötesinde bir değer (200.2 pc'ten küçük ya da 219.8 pc'ten büyük) bir uzaklık vermesi gerekir. Bu durumda kalibrasyonunuzun Gaia ölçümleriyle uyumlu olmadığını ifade eden $H_1$ alternatif hipotezine ilişkin bir kanıta ulaşmış olursunuz.
p Değeri¶
p-değeri (ing. probability value: p-value), sıfır hipotezinin ($H_0$) doğru kabul edilmesi durumunda rastgele seçilen bir örneklemin ortalamasının belirli bir değerden büyük ya da ona eşit olması olasılığıdır. Sıfır hipotezinin ($H_0$) doğru olması durumundaki olasılık belirlendiği için koşullu bir olasılıktır. Özünde bu değer bir hipotezin doğru olduğu kabulü altında herhangi bir istatistiğin tamamen şans eseri olarak elde edilmesi olasılığını verir. Burada şanstan kasıt, rastgele seçimle örneklem oluşturulurken bu sonucu verecek örneklerin (verilerin) seçilmiş olmasıdır (örneklem hatası, ing. sampling error).
p-değeri belirlenen bir istatistiki anlamlılık derecesinden (ing. statistical significance: $\alpha$) daha küçükse sıfır hipotezi ($H_0$) reddedilir. Çünkü elde edilen sonuç belirlenen düzeyde istatistiksel olarak anlamlıdır. Sıfır hipotezi doğru iken bu sonucun elde edilmesi tamamen şans eseri de (örnekleme hatasıyla) olasıdır. Ancak bu olasılık istatistiksel anlamlılıkla izin verilenden daha azdır. Bu durumda alternatif hipotez ($H_1$) için bir kanıt elde edilmiş olunur. Eğer p-değeri belirlenen istatistiki anlamlılılk derecesinden büyükse $H_0$ reddedilemez (bu kabul edildiği anlamına gelmez!). Çünkü sonuç belirlenen düzeyde istatistiksel olarak anlamlı değildir. Bu sonucun sıfır hipotezinin doğru olması durumunda şans eseri elde edilmiş olma ihtimali belirlenen istatistiksel anlamlılık seviyesinden büyüktür ki zaten sıfır hipotezi rastgele seçilen bir örneğin doğrulaması üzerine kuruludur.
Örneğin p-değeri $0.03$ olarak hesaplanıyor ancak belirlenen istatistiksel anlamlılık değeri $\alpha = 0.05$ ise, bu $H_0$'ın doğru olması durumunda elde edilen sonucun tamamen şans eseri elde edilme olasılığının $\%5$’ten küçük bir olasılığa sahip olduğu ($\%3$) anlamına gelir. Bu nedenle $H_0$ reddedilir. Bir başka deyişle, hesaplanan değer gözlenen olgunun örneklem hatasıyla oluşma ihtimalini $\%3$ olarak belirlediği için $H_0$ reddedilir. Ancak p-değeri (örneğin biraz sonra göreceğimiz t-testi ile) 0.1 olarak hesaplanırsa bu böyle bir örneklemin $H_0$’ın doğru olması durumunda $\%10$ ($> \alpha = 0.05$) ihtimalle rastgele herhangi bir örneklem için elde edilebileceğni gösterir ki bu durumda $H_0$ reddedilemez.
Örnek: p Değeri¶
Diyelim ki açık yıldız kümelerinin içerdiği yıldızların %25’inin G tayf türünden daha sıcak yıldızları içerdiği gibi bir iddia var ve biz de bir açık yıldız kümesi gözleyerek bu iddiayı test etmek istiyoruz.
Bu durumda sıfır hipotezimiz
- $H_0$: $\mu = 0.25$,
alternatif hipotezimiz ise
- $H_1$: $\mu \neq 0.25$
olarak belirlenebilir.
Şimdi yapmamız gereken gözlediğimiz kümedeki G tayf türünden daha sıcak yıldızları saymak (diyelim ki 120 yıldızdan 40’ı böyle yıldızlar olsun) ve $H_0$ ‘ın doğru olduğu kabulü altında bu gözlem sonucunun tamamen şans eseri elde edilip, edilemeyeceğini bulmak!
Bu deney yine bir Bernoulli deneyidir. Zira kümdede ele alınan bir yıldız ya G tayf türünden daha sıcaktır. ya da değildir. Bu durumda standart sapma $\sigma = \sqrt{n p q} = \sqrt{120~0.25~(1-0.25)} = 4.74$ olarak bulunur. %25 G tayf türünden daha sıcak yıldız bulma olasılığı ise $\mu = 0.25~x~120 = 30$ yıldıza karşılık gelir.
$$ z = \frac{x - \mu}{\sigma} = \frac{40 - 30}{120~0.25~0.75} = 2.108 \sim 2.11$$
Şimdi iş bu z-değerinin karşılık geldiği p-değerini hesaplamaktır. Bu değer bir standart normal dağılımı tablosundan (z-tablosu) alınabileceği gibi scipy.stats.norm fonksiyonlarıdan cdf() da bu amçla kullanılabilir. Ancak cdf() (1 - p) değerini vereceği için, doğrudan p-değerini veren sf() fonksiyonu (ing. survival function) da tercih edilebilir.
from scipy import stats as st
z = 2.108
p = st.norm.sf(z)
print("z = {:.2f} degeri p-degeri = {:.4f}'e karsilik gelir".\
format(z, p))
olasilik = st.norm.cdf(z)
print("z = {:.2f} degerinin karsilik geldigi olasilik = %{:.2f}'tir.".\
format(z, 100*olasilik))
$z = 2.11$’e karşılık gelen p-değeri $p = 0.0175$ bulunur. Yani bu örneğin bu popülasyondan örnekleme hatasıyla türetilmiş olma olasılığı %1.75 olup istatistiksel anlamlılığın $\alpha = 0.05$ seçilmesi durumunda $p = 0.0175 < \alpha$ olduğu için $H_0$ reddedilir. Yani, ele alınan örnek iddianın doğru kabul edilmesi durumunda istatistiksel anlamlılık seviyesi dahilinde tüm açık yıldız kümeleri popülasyonlarından rastgele bir seçimle türetilmiş olabilir ama bu ihtimal "göreli olarak (istatistiksel olarak anlamlı sayılacak düzeye göre) düşüktür". Bu nedenle de elimizdeki yıldız grubu (örneklem) tüm açık yıldız kümelerindeki G-tayf türünden yıldızların (popülasyon) oranının $\%25$ olduğu iddiasıyla çelişmektedir! Ama biz evrendeki tüm açık yıldız kümelerini gözleyemeyiz ve hiçbir zaman da gözleyemeyeceğiz! O yüzden, $H_1$ kabul edilmiş de değildir ama o yönde bir kanıt elde edilmiştir. Bu sonuç daha fazla örnekle ya da istatistiksel anlamlılık seviyesi değiştirilerek, daha yüksek bir güven düzeyinde konuşmak istendiğinde yanlışlanabilir. Dolayısı ile bu güven düzeyinde $H_0$'ı reddetmiş olmakla şöyle bir yanlış yapıyor olabiliriz: Çok küçük bir ihtimalle de olsa bu veri $H_0$ (tüm yıldızların %25'inin G-tayf türünden sıcak yıldızlar olduğu) doğru olduğu halde örnekleme hatasıyla ("talihsizlik!") oluşmuş da olabilir! Yani 120 yıldızlık örneğimizde 40 yıldızın G-tayf türünden olması bu önerme doğru olduğu halde örneklem hatası (şans kaynaklı) nedeniyle düşük ihtimaldir! Bu durumda Tip-I Hatası yapılmış olur.
Ancak %99'luk güvenlik düzeyine denk gelen $\alpha = 0.01$ seçilseydi $p > \alpha$ olduğu için $H_0$ reddedilemezdi. Bir başka deyişle, elimizdeki örnek açık yıldız kümelerindeki yıldızların %25’inin G tayf türünden daha sıcak yıldızları barındırdığı iddiasını $\alpha = 0.01$ istatistiksel anlamlılık ya da $\%99$ güven seviyesinde reddetmez; çünkü bu örneklem, $H_0$'ın doğru olması durumunda, istediğimiz ya da izin verdiğimiz anlamlılık seviyesinden (%1) daha yüksek bir olasılığa (%1.75) sahiptir. Yine benzer şekilde, bu önermeyi reddetmemekle bir hata yapıyor olabiliriz. Yani $H_0$ hipotezi reddedilmediği halde, onun örneklemimiz için öngördüğü sayıya ($30$) istatistiksel anlamlılığın izin verdiği kadar yakın sayıda sıcak yıldız içeren bir açık kümeyi tamamen şans eseri (örneklem hatası sonucu) seçmiş de olabiliriz. Yani $H_0$ doğru olmayabilir ama talihsizlik eseri seçtiğimiz örnek onun öngördüğüne yakın bir sonuç vermiştir. Bu durumda da Tip-II Hatası yapılmış olur.
Şimdi seçtiğimiz istatistiksel anlamlılık seviyesinin ($\alpha = \%1 = 0.01$)'in karşılık geldiği yıldız sayılarına bakalım. Elimizdeki örnekte kaç yıldızdan daha az ya da daha çok sayıda G-tayf türünden daha sıcak yıldız olması $H_0$'ı reddetmemize yol açar, çünkü belirlediğimiz istatistiksel anlamlılık seviyesinde $H_0$'ın öngördüğünden daha az ya da daha çok G-tayt türünden sıcak yıldız içerir? Bu iki yönlü bir testtir. İstatistiksel anlamlılık nihayetinde bir olasılık olasılık olduğu ve belirleyeceğimiz yıldız sayılarından az ve çok yıldız içerme olasılıkları toplanacağı için, bu olasılık, eğrinin altında ortalamadan ($30$) her iki yönde belirlemek istediğimiz sayıların dışında kalan alandır. Bu sayıları scipy.stats.norm.interval() fonksiyonu ile belirleyebiliriz. Ancak bu fonksiyonun eğrinin $H_0$'ın kabul edildiği alanı veren uç değerleri döndürdüğüne dikkat edilmelidir. Bu nedenle istenen $\alpha$ değerinin karşılık geldiği uç noktalar için, o olasılığın 1'den farkı alınarak interval() fonksiyonuna sağlanması gerekir. Bu değer aynı zamanda güven düzeyini ifade eder.
### Secilen istatistiksel anlam seviyesinin karsilik geldigi
# yildiz sayilari
alpha = 0.01
olasilik = 1-alpha
sayi = st.norm.interval(olasilik,loc=30,scale=4.74)
print("p = {:.2f} degerinin karsilik geldigi yildiz sayilari {:.2f} ve {:.2f}'dir.".\
format(olasilik, sayi[0], sayi[1]))
Görüldüğü gibi, açık yıldız kümemizin, $H_0$ doğru iken belirlediğimiz istatistiksel anlamlılık seviyesinde (%1) 42.21'den fazla, 17.79'dan az sayıda G-tayf türünden sıcak yıldızı şans eseri (örnekleme hatasıyla) içerme olasılığı %1'in altında olmalıdır. Eğer sözgelimi bizim incelediğimiz kümede 45 böyle yıldız olsaydı, "rastgele bir seçimde bu kadar yıldız içeren bir kümeye denk gelme olasılığımız %1'in altında o nedenle $H_0$ doğru olamaz!" derdik. Oysa ki bizim kümemizde 40 tane bu tür yıldız var. O nedenle $H_0$ bu anlamlılık seviyesinde reddedilemez!
Oysa ki $H_0$'ı reddettiğimiz ilk anlamlılık seviyesi olan $\alpha = \%5 = 0.05$ için bu sayıları hesaplasak;
### Secilen istatistiksel anlam seviyesinin karsilik geldigi
# yildiz sayilari
alpha = 0.05
olasilik = 1-alpha
sayi = st.norm.interval(olasilik,loc=30,scale=4.74)
print("al = {:.2f} degerinin karsilik geldigi yildiz sayilari {:.2f} ve {:.2f}'dir.".\
format(olasilik, sayi[0], sayi[1]))
elde ederdik. Görüldüğü gibi bu anlamlılık seviyesinde $39.29$'dan daha fazla sıcak yıldız içeren bir açık kümeyi tamamen şans eseri oluşturma olasılığı $\%5$. Yani $H_0$ doğruyken rastgele seçtiğimiz bir kümenin 40 yıldız içeren bir küme olma olasılığı bundan daha az.
t Dağılımı¶
Tüm parametrelerin değerleri Gaussian olasılık dağılımına uygun dağılmayabilir. Bazı parametreler, Gaussian dağılımına uyuyor gibi görünse de ortalamadan uzak değerlerinin olaslıkları Gaussian dağılımla öngörülenden yüksek bulunablir. Önemli bir başka durum yeterince örnek sayısının olmadığı durumdur ki bu durumda ortalama ve standart sapma çok iyi belirlenmememiş olabilir. Bu durumlarda parametre değerlerinin olasılık dağılımı için t-dağılımını (ing. Student's t-distribution) kullanmak daha iyi bir seçenektir.
$t$-istatistiği, $x$ değerinin, örneğin ortalama değeri olan $\bar{x}$ değerinden kaç standart sapma ($s$) uzaklıkta olduğunu belirler ve aşağıdaki şekilde ifade edilir.
$$ t = \frac{|x - \bar{x}|}{s} $$
Bu durumda $t$-istatistiği için olasılık dağılımı
$$ P(t,\nu) = \frac{1}{\sqrt{\nu~\pi}} \frac{\Gamma[(\nu + 1) / 2]}{\Gamma(\nu / 2)} (1 + \frac{t^2}{\nu})^{-(\nu + 1)/2} $$
dir. Burada $\nu$ serbestlik derecesini gösterir ve $\bar{x}$ $N$ bağımsız ölçümün ortalamasıysa $\nu = N - 1$ ile verilir.
Aşağıdaki görselde, serbestlik derecesinin değişmesi ile t-dağılımının değişimi görülmektedir. Mavi ile gösterilen eğri standart normal dağılımdır. 1 serbestlik derecesine sahip olan t-dağılımında görüleceği üzere, $t$-dağılımı, daha geniş kanatlara sahiptir. Bu normal dağılımın ortalama değerinden daha uzakta ortalama değerlerin bulunma olasılığının normal dağılıma göre daha fazla olduğunu göstermektedir; ki az sayıda ölçüm söz konusu olduğunda bu yaklaşım daha gerçekçidir. Her bir serbestlik derecesi için verilen dağılımda önceki serbestlik dereceleri için t-dağılımı yeşille, verilen serbestlik derecesi için kırmızıyla gösterilmektedir.
Serbestlik derecesi arttıkça ya da başka bir deyişle örneklem büyüdükçe, t-dağılımı da normal dağılıma yakınsamaktadır. Bunun sonucu olarak, örneklemdeki eleman sayısının az olması durumunda t-dağılımının kullanılması, normal dağılımın kullanılmasından daha uygun olmaktadır.
Örnek: Gaia Paralaksları t dağılımı Yaklaşımı¶
Henüz Gaia görevinin ancak belirli bir bölümünü (DR-2) tamamlamış durumdadır ve aynı yıldızın paralaks ölçümünü belirli kez yapabilmiştir. Daha önceki örnekte gördüğümüz yıldızın 9 paralaks ölçümü yapılmış ve uzaklğı parsek cinsinden hesaplandığında elde edilen ortalama uzaklık değeri ($\bar{x}$) ve belirsizliği ($s$) $210 \pm 5 ~pc$ olarak veriliyor olsun. %95 güven seviyesine karşılık gelen aralığı bu kez $t- dağılımı$ kullanarak bulalım.
Bu kez ölçüm sayısı çok az olduğu ($N = 9$) için serbestlik derecesini dikkate alan t-dağılımı üzerinden bir güven aralığı tesis etmek istiyoruz. t-dağılımı için verilen kümülatif olasılık tablosundan %95 güven seviyesine $ν = 9 – 1 = 9 – 1 = 8$ serbestlik derecesi için karşılık gelen $t-istatistiği$ değeri $t = 2.3$’tür. Bu durumda;
$$ t = \frac{|x - \bar{x}|}{s} = 2.3 $$ $$ t = \frac{|x - \bar{x}|}{s} = 2.3 \Rightarrow x = \bar{x} \pm t~s \Rightarrow x = 210 \pm 2.3~x~5 = 210 \pm 11.5~pc$$
Bu durumda %95 güven seviyesine karşılık gelen aralık:
$$ x = (210 – 11.5, 210 + 11.5) = (198.5, 221.5) $$
aralığıdır.
Aynı güven seviyesi için çok daha geniş olarak bulunan güven aralığı, t-dağılımının, Gaussyen (normal) dağılımın azımsadığı (ing. underestimate) uç değerlerin olasılıklarını daha büyük bulmasından kaynaklanmaktadır. Örnek sayısı az olduğu vakit, uç değerlerin olasılıklarının daha yüksek olması beklenebilir. Örnek sayısı artıp, serbestlik derecesi buna bağlı olarak yükseldikçe dağılım, Gaussyen’e yaklaşacaktır. Örneğin $\nu = 50$ serbestlik derecesi için $t \sim 2.0$ iken normal dağılımda güven aralığının hesaplanması için kullanılan $z$ parametresi de $z \sim 1.95$ değerini almaktadır. Oysa $N = 9$ olan örneğimizde $t = 2.3$ iken aynı güven aralığında $z = 1.96$’dır.
Aynı örneği Python'da scipy paketi fonksiyonlarından t-dağılımına ilişkin araçlar sunan scipy.stats.t() ile çözmeye çalışalım. Fonksiyon scipy kütüphanesinde yer alan, Ders-2b Dağılım Fonksiyonları'nda gördüğünüz diğer dağılım fonksiyonlarına benzer şekilde çalışır ve benzer metotlarla parametrelere sahiptir. t-dağılımı sürekli bir fonksiyondur bu nedenle süreksiz fonksiyonlardaki Olasılık Kütle Fonkskiyonu'nun (ing. Probability Mass Function, PMF) yerini Olasılık Yoğunluk Fonksiyonu (ing. Probability Density Function, PDF) alır. Aynı şekilde kümülatif olasılıklar da yoğunluk fonksiyonu (ing. Cumulative Density Function, CDF) ile ifade edilir.
Öncelikle problemimizin parametreleriyle %95'lik bir güven düzeyi için güven aralığı oluşturalım. Bu amaçla güven düzeyini (0.95) alpha parametresine (bu parametreinin ismi istatistiksel anlamlılığı simgeleyen $\alpha$ ile aynı olsa da güvenlik düzeyini doğrudan vermek gerektiğine dikkat ediniz), serbestlik derecesini ($\nu = N - 1 = 9 - 1 = 8$) stats.t fonksiyonunun df parametresine, ölçümlerin ortalamasını ($\mu = 210~pc$) loc parametresine, standart sapmasını ise scale parametresine vermeliyiz.
from scipy import stats as st
N = 9
xbar = 210 # pc
s = 5 # pc
gd = 0.95 # %95
guven_araligi = st.t.interval(alpha=gd, df=N-1, loc=xbar, scale=s)
print("%{:.1f}'lik guven araligi {:.1f} ile {:.1f} arasidir.".\
format(100*gd, guven_araligi[0], guven_araligi[1]))
t-dağılımına uyan bireysel ölçümleri (örneklemi) oluşturabilmek için rvs() fonksiyonunu kullanalım. Bunun için örneklem büyüklüğümüzü size parametresi ile $N = 9$'a, serbestlik derecesini df parametresi ile $\bar{x} = N - 1 = 8$'e, ortalama değeri loc parametresinden $\mu = 210~pc$ ve standart sapmayı da scale parametresinden $s = 5~pc$'e ayarlayalım.
from scipy import stats as sta
import numpy as np
# Ornegi tekrar ayni degerlerle uretebilmek icin
# seed edelim.
np.random.seed(53)
N = 9
xbar = 210 # pc
s = 5 # pc
paralakslar = st.t.rvs(df=N-1, loc=xbar, scale=s, size=N)
print(paralakslar)
Bir t-dağılımından verilen bu parametrelerle üretilen paralaks ölçümleri için doğal olarak ortalama ve standart sapma değerleri verilenlerle birebir aynı olmayacaktır. Zira işin içine rastgelelik katılmıştır. Ancak bu örnekte bireysel Gaia ölçümleri ürettiğimiz veri setimizi ($paralakslar$) kullanabiliriz. İstediğimiz parametrelerle bir t-dağılımından rastgele örneklem seçerek oluşturduğumuz bu t-dağılımına en iyi uyan olasılık yoğunluğu fonksiyonunu stats.t.fit fonksiyonuyla bulmaya çalışalım.
t_uyumlama = st.t.fit(paralakslar)
print("Ort: {:.2f}, St.Sapma: {:.2f}".\
format(t_uyumlama[1], t_uyumlama[2]))
İstenirse bu uyumlama (fit) fonksiyonunun parametrelerine karşılık gelen bir olasılık yoğunluk fonksiyonu (PDF) stats.t.pdf() fonksiyonu ile çizdirilebilir (sürekli mavi eğri). Bu fonksiyonunun yanı sıra karşılaştırma için başlangıçtaki parametrelerimize ($\mu = 210, \sigma = 5$) karşılık gelen bir t-dağılımı çizdirelim (süreksiz mavi eğri). Ayrıca uyumlama fonksiyonu ile elde ettiğimiz parametrelere karşılık gelen bir normal dağılım için de bir olasılık yoğunluk fonksiyonu çizdirelim (turuncu sürekli eğri).
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
xbar_fit = t_uyumlama[1]
s_fit = t_uyumlama[2]
# Daha cok sayida nokta iceren bir x-ekseni tanimlayalim
x = np.linspace(200,225,50)
pdf_fit = st.t.pdf(x, df=8, loc=xbar_fit, scale=s_fit)
pdf_orj = st.t.pdf(x, df=8, loc=xbar, scale=s)
pdf_gauss = st.norm.pdf(x, loc=xbar_fit, scale=s_fit)
plt.hist(paralakslar, density=True, histtype='stepfilled', alpha=0.2)
plt.plot(x, pdf_fit, 'b-', label="Uyumlama")
plt.plot(x, pdf_orj, 'b--', label="Orjinal")
plt.plot(x, pdf_gauss, c="orange", ls="-", label="Normal")
plt.legend(loc="best")
plt.show()
Serbestlik Derecesi¶
Şu ana kadar pek çok kez kullandığımız ve daha önce de tanımladığımız serbestlik derecesi kavramını bu noktada hatırlamakta yarar var. Serbestlik derecesi değiştikçe t-dağılımının şekli (ve doğal olarak olasılıklar değiştiği) için t-tablolarına (ya da t-dağılımına dayanan scipy.stats.t fonksiyonlarına) başvurulurken serbestlik derecesi girdi parametresidir. Serbestlik derecesi (ing. Degree of Freedom, df ya da dof) temel olarak, formal tanımıyla) dinamik bir sistemin, üzerinde herhangi bir kısıtlamayı ihlal etmeden hareket edebileceği bağımsız yolların sayısı olarak tanımlanır ve sistemin eksiksiz olarak tanımlanabileceği bağımsız koordinatların sayısıdır. Bu bağlamda, serbestlik derecesi özet olarak bir sistemde bağımsız şekilde değişebilen değerlerin ya da niceliklerinin sayısıdır. Serbestlik derecesi, öncelikli olarak örnek dağılımın eleman sayısına (örnek büyüklüğüne) bağlıdır. Eleman sayısının fazla olması serbest olarak değişebilecek değerlerin sayısının da fazla olması anlamına gelir. Bir örnek dağılım için serbestlik derecesi; $dof = N -1$ ile veriliir. Burada serbestlik derecesi $dof$, örneklem büyüklüğü ise $N$’dir. Örnek dağılım, bir ana dağılımdan üretilmekte olduğu için, örnek dağılımın ortalama değeri, ana dağılımın ortalamasını temsil etmelidir. Bu durum serbestlik derecesinin 1 azalması anlamına gelir; bu nedenle serbestlik derecesi $N - 1$ olmuştur.
Aynı şekilde bir sistemin varyansı söz konusu olduğunda da serbestlik derecesi $dof = N - 1$'dir zira varyans, ortalamaya bağlı olarak hesaplanır.
Eğri uyumlamasında ise; $dof = N – m$ ile verilir. Burada $m$, uyumlamada kullanılan parametre sayısıdır. Uyumlamada kullanılan parametre sayısının fazla olması serbest olarak değişebilen değerlerin sayısının parametre sayısı kadar azalması anlamına gelir.
Modelleme¶
Bir bilimsel model gerçek bir olgunun fiziksel, kavramsal ya da matemetiksel olarak temsilidir. Bir modelin amacı bir olguyu daha kolay anlaşılır, yorumlanır, tanımlanır, nitel hale getirilebilir, görselleştirilebilir ve simüle edilebilir hale getirmektedir. Böylece bir model üzerinden geleceğe ilişkin tahminler de yapılabilir (ekstrapolasyon). Geçen ders gördüğümüz eğri uyumlama da gözlemsel (ya da deneysel) bir veri setine yapılan analitik bir modeldir.
Bir uyumlama işlemi sonunda, gözlemsel veri ile uyumlama eğrisinin değerleri arasındaki farklara artık (ing. residual) denir. Artıkların dağılımı ve gösterecekleri olası trendler, uyumlamada kullanılan yöntemin veya modelin ne kadar başarılı veya kabul edilebilir olduğunun bir ölçütü olarak kullanılabilir.
Aşağıdaki grafiğin üst panelinde HW Virgo yıldızının minimum zamanlarının (veri noktaları) değişimine yapılan bir parabol uyumlaması (kesikli eğri) görülmektedir. CCD gözlemleriyle elde edilen minimum zamanları daha büyük simgelerle, birinci minimumlar içi dolu, ikinci minimumlar içi boş çemberlerle gösterilmiştir. Alt panelde parabolik bu modelden artıklar görünmektedir. Bu panelde $y = 0$ doğrusu (sürekli doğru) parabole karşılık gelmektedir. Uyumlanan parabolün iki ucundan tutulup yatay bir doğru haline getirildiğini; bu sırada veri noktalarının da bu parabol modele olan uzaklıklarının bozulmadan taşındığnı düşünebilirsiniz. Veri noktalarının alt panelde modele göre bu şekildeki konumları, bir başka deyişle modelden farkları, artıkları verir. Yakın zamanda yapılan CCD gözlemleriyle elde edilen minimum zamanlara yapılan parabol modelden artıkların sinüsoidal bir değişim gösterdiği düşünülmüş, artık (Kepleryan) ikinci bir modelle uyumlanmıştır (kesikli eğri).
Bir Modelin Başarımına İlişkin Parametreler¶
Bir modelin başarımı veri noktaları ile o noktalara modelde karşılık gelen değerler arasındaki fark, daha doğru ifadeyle bu farkların karelerin toplamı, üzerinden tanımlanır. Bu fark ne kadar küçükse model o kadar başarılı adledilir. Ancak bu Artık Kareler Toplamı (ing. Residual Sum of Squares) olarak tanımlanan ve aşağıdaki şekilde ifade edilen bu toplam doğrudan bir model başarımı göstergesi değildir.
$$ S_r = \sum_{i=1}^{n} (y_{i} - \hat{y_i})^2 $$
Kök Ortalama Kare Hatası/Sapması (ing. Root Mean Square Error/Deviation) ya da daha iyi bilinen adıyla "RMS Hatası" artıkların uyumlama etrafında ne kadar saçıldığını gösteren bir parametredir ve yine bir model başarımı göstergesi değildir.
$$ RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n} (y_{i} - \hat{y_i})^2} = \sqrt{\frac{S_r}{n}} $$
RMSE yerine, ona çok benzer şekilde ancak fark kareler toplamının $n$ nokta üzerinden ortalaması yerine varyansa (ya da ölçüm hatalarının karesine) bölünmesiyle tanımlanan $\chi^2$ ve onun serbestlik derecesine normalize edilmesiyle tanımlanan indirgenmiş $\chi^2_{\nu}$ parametresi daha sık kullanılır. Bu parametre adıyla bilinen $\chi^2$-dağılımını göstermesi nedeniyle artıkların analizinde önemli bir avantaj sunar.
$$ \chi^2 = \sum_{i=1}^{n} \frac{(y - y_i)^2}{e_i^2} $$
$$ \chi^2_\nu = \frac{\chi^2}{\nu} $$
Burada sistemin $\nu$ serbestlik derecesi (dof, ing. degree of freedom), veri noktası sayısı ile model paramatreleri sayısı arasındaki farktır ($\nu = n - m$).
Ayrıca korelasyon katsayısı (ing. correlation coefficient, $r^2$) da model başarımını ifade etmek için zaman zaman kullanılsa da daha çok lineer ve polinom modeller için kullanılmaktadır. Hatırlatmak gerekirse; korelasyon katsayısı bir modelin ortalamaya oranla başarısını gösterir ve aşağıdaki ifade ile verilir. $\bar{y}$ ölçümlerin aritmetik ortalaması ve $\hat{y_i}$ modelde $i.$ noktaya karşılık gelen değer olmak üzere;
$$ S_t = \sum_{i=1}^{n} (y_{i} - \bar{y})^2 $$
$$ S_r = \sum_{i=1}^{n} (y_{i} - \hat{y_i})^2 $$
$$ r^2 = \frac{|S_r - S_t|}{S_t} $$
Doğrusal ya da parabolik modelden daha karmaşık modelleri ortalamayla karşılaştırmak iyi bir fikir olmayacağı için bu ölçütün de kullanım alanı bu tür modellerle sınırlıdır.
$\chi^2$ Dağılımı ve $\chi^2$ Testi¶
Serbestlik derecesi ($\nu$) kadar bağımsız standart normal rastgele değişkenin toplamının dağılımıdır. Gamma dağılımının özel bir halidir. Hipotez testi, güven aralığı oluşturma (ing. confidence interval), uyum başarısı (ing. goodness of fit) gibi bir çok istatistiksel çıkarımda kullanılmaktadır. Nadiren Helmert dağılımı olarak da isimlendirilir.
Bir x değişkenin ölçüm değerleri $x_i$ ve üzerindeki belirsizlikler $\sigma_i$ olmak üzere ağırlıklı ortalama:
$$ \bar{x} = \frac{\sum_{j = 1}^{n} (x_j / \sigma_j^2)}{\sum_{j = 1}^{n} (1 / \sigma_j^2)} $$
Tüm ölçümler üzerinde aynı belirsizlik bulunuyorsa ($\sigma = \sigma_i$) bu durumda ağırlıklandırmaya gerek olmadğından ortalama değer doğrudan elde edilir.
$$ \bar{x} = \frac{(1 / \sigma^2) \sum_{j = 1}^{n} x_j}{(1 / \sigma^2) \sum_{j = 1}^{n} 1} = \frac{1}{N} \sum_{j = 1}^{n} x_j $$
Ortalamanın Standart Sapması (ya da Varyansı) her bir veri noktasının hatası ($\sigma_i$) ya da tüm verinin standart sapmasından ($\sigma$) farklıdır:
$$ \sigma_{\nu}^2 = \frac{1}{\sum_{j = 1}^{n} (1 / \sigma_j^2)} $$
Tüm ölçümler üzerinde aynı belirsizlik bulunuyor ya da verinin tamamının standart sapması biliniyor ise ($\sigma = \sigma_i$):
$$ \sigma_{\nu}^2 = \frac{\sigma^2}{N} $$
Ölçümün (ya da ölçüm aletinin) belirsizliği:
$$ \sigma^2 \sim s^2 = \frac {1}{N-1} \Sigma_{j = 1}^n (x_j - \bar{x})^2 $$
Sınırlı sayıda ölçüm alabilmekten kaynaklanan istatistiksel dalgalanma (Poisson):
$$ \sigma^2 = \nu \sim \bar{x} $$
olmak üzere;
$\chi^2$-testi: $x_j$ olası değerleri için gözlenen bir frekans dağılımının $h(x_j)$, varsayılan bir olasılık dağılımına göre belirlenmiş olasılıklarını veren $N P(x_j)$ dağılımı ile karşılaştırılmasına dayanan testtir. Bu karşılaştırma için kullanılan $\chi^2$-istatistiği:
$$ \chi^2 = \sum_{j=1}^{n} \frac{[h(x_j) - N~P(x_j)]^2}{\sigma_j(h)^2} $$
Serbestlik derecesi ($\nu$): Veri sayısı (N) ile verilerden birbirlerinden bağımsız olarak belirlenen parametres sayısı (m) arasındaki farktır (N – m).
İndirgenmiş $\chi^2_{\nu} = \chi^2 / \nu$: $1$’e yakınsa $h(x_j)$ dağılımı $N~P(x_j)$ dağılımı ile uyumludur. Serbestlik derecesine karşılık verilen indirgenmiş $\chi^2_{\nu}$ tablolarından rastgele seçilen bir veri setinin, ana dağılımdan türetilmiş olup olmadğının olasılığı bulunur.
Testin Uygulanışı¶
$\chi^2$-testi genellikle gözlemsel verinin, gözlenen olaya ilişkin teori ile tutarlılığının sınanması için kullanılır. Genellikle testte sıfır hipotezi ($H_0$) “teorik yaklaşım gözlem verisi ile tutarlıdır.” şeklinde seçilir.
Sıfır hipotezi ($H_0$) doğru ise uyumlama sonrası artıkları, rastgele hatalardan kaynaklanır ve bu nedenle normal dağılım gösterir. Bu durumda her bir gözlem verisinin uyumlanan eğriden (ya da modelden) farklarının toplamı bir ki-kare dağılımı vereceğinden ki-kare testi uygulanabilir.
Test adımları şu şekilde uygulanır:
$\chi^2$-değeri (Pearson test istatistiği, $\chi^2$) aşağıda verilen en genel ifadeden (ya da daha önce verilen formüller kullanılarak) hesaplanır. Burada $O_i$: gözlenen değeri, $C_i$ modelden hesaplanan değeri, $\sigma_i$ gözlemsel verinin hatasını (ölçüm hatası ya da standart sapmasını) göstermektedir: $ \chi^2 = \sum \frac{(O_i - C_i)^2}{\sigma_i^2} $
Serbestlik derecesi ($\nu = N - m$) hesaplanır. Burada $N$ ölçüm sayısını, $m$ modelin parametre sayısını göstermektedir.
Sıfır hipotezi ($H_0$) kurgulanır. Sıfır hipotezi, model uyumlamalarında seçilen modelin veriyi uyumladığını ifade ederken, alternatif hipotez ise bunun tersini ifade eder. $\chi^2$ testinde test edilen $H_0$ sıfır hipotezidir.
Bir güven düzeyi seçilir.
Ki-kare ($\chi^2$) tablolarından ilgili serbestlik derecesi ve güven düzeyine karşılık gelen kritik $\chi^2$ değeri, hesaplanan $\chi^2$ ile karşılaştırılır.
Hesaplanan $\chi^2$ değeri, tablodaki kritik değerden küçük ise sıfır hipotezi (gözlem verisinin teori ya da model ile tutarlı olduğu) reddedilemez. Eğer büyük ise sıfır hipotezi reddedilir (teori ya da model gözlem ile tutarlı değildir görüşü için bir kanıt elde edilir).
$\chi^2$ testini verilen güven düzeyleri için kritik $\chi^2$ değerlerini tablolardan okuyup hesaplanan $\chi^2$ değeri ile karşılaştırarak yapmak yerine elde edilen $\chi^2$ değeri için, $H_0$'ın doğru olması durumunda o değerden daha büyük bir $\chi^2$ elde etme olasılığını veren p-değerini hesaplayıp ($p = P(\chi^2 > \chi^2_{hesap})| H_0$), bu değeri istatistiki anlamlılık seviyesi ile ($\alpha$ = 1 - güven düzeyi) karşılaştırmak mümkündür. Eğer hesaplanan bu p-değeri ile istatistiksel anlamlılık seviyesinden ($\alpha$) büyükse $H_0$ reddedilemez. Aksi durumda $H_0$ reddedilir. Bu durum $H_1$ için bir kanıt teşkil eder.
$\chi^2$ değeri aynı parametre sayısına sahip iki modelin doğrudan karşılaştırılması için bir gösterge olarak da kullanılabilir.
Model uyumlamalarında serbestlik derecesi başına $\chi^2$ değeri ("indirgenmiş ki-kare", $\chi^2_\nu = \chi^2 / \nu$) kullanılıyorsa modelin başarımı $\chi^2_\nu$ değeri 1 ile karşılaştırılarak belirlenir.
$\chi^2_{\nu} = 1$ ise uyumlama gözlem verisinin üzerindeki hata mertebesinde iyi sonuç vermektedir. $\chi^2_{\nu} > 1$ ise uyumlama gözlem verisinin üzerindeki belirsizliği tam olarak modelleyememektedir. $\chi^2_{\nu} << 1$ ise uyumlama çok zayıftır. $\chi^2_{\nu} < 1$ gözlem verisi üzerindeki hata büyüktür sonuçlarına varılır.
Örnek: Doğru Uyumlama¶
$\chi^2$-testinin nasıl uygulandığını görmek için iyi ve basit bir uygulama doğrusal bir modeli gözlemsel veri üzerine uyumlamaktaır. Bu dersle birlikte verilen gözlemsel veriye bir doğru uyumlayalım ve Pearson $\chi^2$-testiyle uyumun başarısını (ing. goodness of fit) görelim.
Öncelikle veri setimizi inceleyelim.
from astropy.io import ascii
veri = ascii.read("veri/model_karsilastirma_chi2.dat")
veri
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
x = np.array(veri['col1'])
y = np.array(veri['col2'])
yhata = np.array(veri['col3'])
plt.errorbar(x,y,yerr=yhata,fmt=".",c="k")
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Öncelikle basit bir model için doğrusal en küçük kareler yöntemini seçerek basit bir Python fonksiyonuyla doğrusal uyumlama yapalım. numpy.polyfit() fonksiyonu deg=1 parametresiyle çalıştırıldığında bu amaca bizi ulaştıracaktır.
katsayilar = np.polyfit(x, y, deg = 1) # Sondaki 1 dogru uyumlamasini gosterir
m = katsayilar[0]
n = katsayilar[1]
print("Verilen veriye uyumlanan en iyi dogru: y = {:.4f}x + {:.4f}".format(m,n))
ymodel = m*x + n
plt.errorbar(x,y,yerr=yhata,fmt=".",c="k", label="veri")
plt.plot(x,ymodel,'r-', label="model")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(loc="upper left")
plt.show()
Şimdi $\chi^2$ değerini kendi yazacağımız bir fonksiyonla hesaplayalım. Fonksiyonumuz $y$ ve $y_{hata}$ dizilerinin yanı sıra modelden hesaplanan $y_{model}$ değerlerini de gönderelim.
import numpy as np
kikare = lambda y,ye,ymod : ((y - ymod)**2 / ye**2).sum()
chi2 = kikare(y,yhata,ymodel)
dof = len(y) - len(katsayilar) # serbestlik derecesi
chi2_nu = chi2 / dof
print("Serbestlik derecesi: ", dof)
print("Chi^2 = {:.2f}, Indirgenmis Chi^2 = {:.2f}".format(chi2,chi2_nu))
Gayet basit bir fonksiyonla hem $\chi^2$, hem de $\chi^2_{\nu}$ istatistiklerini hesapladık. Şimdi bu istatistikleri kullanarak bir hipotez testi yapalım ve modelimizin başarısını bir hipotez testiyle test edelim.
- $H_0$: Lineer model gözlemsel veri ile tutarlıdır.
- $H_1$: Lineer model gözlemsel veri ile tutarlı değildir.
Bunun için öncelikle bir istatistiksel anlam seviyesi seçmeliyiz. $\chi^2$ testlerinde güven düzeyi yerine istatistiksel anlam ($\alpha$) tercih edilir. Ancak biz istatistiksel anlamın güven düzeyinin 1'den farkı olduğunu biliyoruz. Bu örnek için güven düzeyini %90 ($0.90$) olarak belirlemiş olalım. Bu güven düzeyine karşılık gelen istatistiksel anlam seviyesi ($\alpha = 0.10$) ve $\chi^2$ değerini kullanarak $H_0$'ı test edelim. Bu amaçla $\chi^2$ dağılım tablolarını kullanabilirsiniz. Böyle bir tablo kullanılırsa 0.10 istatistiksel anlam seviyesi için kritik $\chi^2$ değerinin $13.36$ olduğu ve hesaplanan $\chi^2 = 4.28 < 13.36$ olduğu görülebilir. Bu durumda lineer modelin gözlemsel veri ile tutarlı olduğuna ilişkin $H_0$ hipotezinin %90 güven düzeyinde reddedilmemesi gerekir. Bu durumda $H_0$ reddedilmez. Yani "lineer modelin gözlemsel veri ile tutarlı olması durumunda bu veri setini şans eseri elde etme olasılığı belirlenen istatistiksel anlamlılık seviyesinden büyüktür!". Aksi halde reddedilir ve bu da lineer modelin gözlemsel veriyle tutarlı olmadığını ifade eden $H_1$ için bir kanıt teşkil ederdi.
Python'da scipy.stats paketinin bu tür testler için geliştirilmiş fonksiyonlara sahip olduğunu daha önceki dağılımlar için yaptığımız uygulamalarda görmüştük. Yukarıda tablolarla yapılan test için scipy.stats.chi2.ppf fonksiyonunu da (ing. point percent function) kullanabiliriz. Ancak ppf fonksiyonuna verilen olasılık değerinin her zaman kritik değere kadarki ($\chi^2 < \chi^2_{kritik}$) olasılık olduğunu görmüştük. Dolayısı ile ppf() fonksiyonunun doğru kritik değeri verebilmesi için ona $\alpha$ istatistiksel anlamlılık seviyesi yerine güven düzeyi sağlanmalıdır ($1 - \alpha$).
from scipy import stats as st
alpha = 0.10
chi2_kritik = st.chi2.ppf(1-alpha, df=dof)
print("Kritik Chi^2 degeri: {:.4f}".format(chi2_kritik))
Görüldüğü gibi kritik $\chi^2$ değeri 13.3616'dır; Bu değer hesaplanan $\chi^2$ değeri olan 4.28'den büyüktür. Bir başka deyişle ($\chi^2 < \chi^2_{kritik}$). Bu durumda $H_0$ reddedilemez. Yani lineer modelin veriyle uyumlu olduğuna ilişkin bir kanıt elde edilmiş olur. Belirlediğimiz istatisiksel anlam seviyesi olan %10, lineer modelin veriyle uyumlu olması durumunda rastgele seçilecek 100 veri setinin 10 tanesinde uyumladığımız doğrunun $13.3616$'dan büyük bir $chi^2$ değeri ($\chi^2 > 13.3616$) verirken 90'ında bundan daha küçük bir $\chi^2$ değeri elde edileceğini ortaya koymaktadır.
$\chi^2$ testini verilen güven düzeyleri için kritik $\chi^2$ değerlerini tablolardan okuyup (ya da scipy.stats.chi2.ppf ile bulup) hesaplanan $\chi^2$ değeri ile karşılaştırarak yapmak yerine elde edilen $\chi^2$ değeri için, o değerden daha büyük bir $\chi^2$ elde etme olasılığını veren p-değerini hesaplayıp, bu değeri istatistiki anlamlılık seviyesi ile de karşılaştırabiliriz. Hesaplayacağımız p-değeri, lineer modelimiz doğruysa ($H_0$) bulduğumuz $\chi^2$ değeri'nden büyük bir artık değerini rastlantısal olarak elde etme ihtimalini verir $p = P(\chi^2 \ge 4.28 | H_0$). Eğer bulduğumuz p-değeri seçtiğimiz güven düzeyiyle belirlediğimiz istatistiksel anlamlılıktan ($\alpha = 1 - 0.90 = 0.1$) küçük çıkarsa ($p < \alpha$) $H_0$'ı reddederiz. Aksi halde $H_0$'ı reddedemeyiz.
scipy.stats.chi2.cdf() fonksiyonu bize $\nu = 8$ için $\chi^2$'in 4.28'den küçük olma olasılığını verecektir. Eğer bunun yerine $\chi^2$ ile karşılaştırmak istediğimiz p-değerini elde etmek istiyorsak ($P)\chi^2 > 4.28~|~H_0$) bu durumda scipy.stats.chi2.sf() fonksiyonu kullanılmalıdır.
olasilik = st.chi2.sf(chi2, df=dof)
print("{:d} serbestlik derecesi icin P(chi^2 > {:.2f}) = {:.4f}".\
format(dof, chi2, olasilik))
Sonuç olarak $\chi^2 > 4.28$ olma olasılığnı %83.1 olarak bulmuş olduk. $p = 0.831$ olan bu değeri $\alpha = 0.10$ ile karşılaştırdığımızda $p > \alpha$ olduğunu görüyoruz. Bu durum $H_0$'ı reddetmememizi gerektirir (tablo çözümüyle aynı sonuca ulaştık).
p-değeri ($p = 0.831$) bize lineer model doğruysa ($H_0$) elde edilen $\chi^2 = 4.28$) değerinden daha büyük bir $\chi^2$ değeri verebilecek bir veri setini şans eseri türetme olasılığını verir.
Peki bu bulduğumuz olasılık değeri (p = %83.10) ne anlama gelmektedir?
Belirlenen parametreleriyle bu lineer modelin doğru olması durumunda modelden en az bu kadar artık verebilecek (modelden en az bu kadar uzakta) bir gözlemsel (deneysel) veri seti elde etme (ölçme) olasılığımız %83.10'dur!
Bir başka deyişle, lineer modelin doğru olması durumunda bu modelden en fazla bu kadar artık verebilecek bir veri setini tamamen şans eseri üretme olasılığı %16.90'dır.
Diğer taraftan bulduğumuz bu p-değeri şu anlamlara gelmez:
- Lineer modelin doğru olma olasılığı %83.1'dir!
- Lineer modelin yanlış olma olasılığı %16.9'dur!
- Lineer modeli rededersek %16.9 ihtimalle yanılıyoruz demektir!
Zira $\chi^2$ testi aslında hiçbir modelin ne kadar olası olduğunu söylemez! Verilen güven düzeyinde ilgili modelle gözlenen artıklar kadar artık verebilecek bir veri setinin elde edilme olasılığını söyler!
Örneğimizde lineer modelin çok başarılı olmadığı da ayrıca görülmektedir. Yukarıda kendi fonksiyonumuzu yazarak elde ettiğimiz bu istatistiki parametreleri LMFIT paketi fonksiyonlarıyla da elde edebiliriz. Bunun için öncelikle doğru modelinin parametrelerini tanımlamamız gerekir.
import lmfit as lm
dogru = lambda x,egim,kesme : egim*x + kesme
dogru_modeli = lm.Model(dogru)
print(dogru_modeli.param_names)
print(dogru_modeli.independent_vars)
Daha sonra uyumlama için makul birer başlangıç değeri vermemiz gerekir. Bu değerleri verinin basit bir grafiğinden kabaca belirleyebiliriz. Ayrıca weights parametresini kullanarak ölçüm hatalarına göre bir ağırlıklandırma da sağlayabiliriz.
sonuc = dogru_modeli.fit(y, x=x, egim=0.25, kesme=1.00, weights=1/yhata)
print(sonuc.fit_report())
Gördüğünüz gibi $\chi^2$ ve $\chi^2_{\nu}$ değerleri bizim basit $kikare$ fonksiyonumuzla hesapladıklarımızla aynıdır. $\chi^2 = 4.28$ değerini ve $dof = 10 - 2 = 8$ değerini scipy.stats.chi2 fonksiyonlarında kullanarak bir $\chi^2$ testi yapabileceğimiz gibi $\chi^2_{\nu} \sim 0.54 < 1$ değerinden hareketle uyumlamanın çok başarısız olmamakla birlikte gözlemsel verideki hataların (ya da saçılmanın) modelin başarımının etkileyecek düzeyde büyük olduğu sonucuna varabiliriz.
F Dağılımı ve F Testi¶
F-Testi, bir gözlem verisine yapılan iki farklı model uyumlamasının karşılaştırılıp hangisinin istatistiksel olarak daha başarılı bir uyumlama olduğunun belirlenmesi için kullanılan sınamalardan birisidir.
$F$ adı, varyans analizinin temellerini atan istatistikçi ve biyolog Ronald A. Fisher’a ithafen verilmiştir. Model seçimi için yapılan F-testinde takip edilmesi gereken temel adımları şöyledir:
Öncelikle bir sıfır hipotezi ($H_0$) belirlenir. $H_0$, genellikle aralarında karşılaştırma yapılan iki varyansın (modelden ya da ortalamadan fark kareler toplamının serbestlik derecesine bölümü) eşit olduğunu ifade eder ($F = s_{1} / s_{2} = 1$). Alternatif hipotez ($H_1$) ise bunun tam tersidir.
Karşılaştırılması istenen iki uyumlamanın artık kareler toplamı (ing. residual sum of squares, $S_r$) hesaplanır.
Rakip uyumlamalar için serbestlik dereceleri hesaplanır.
İstatistiksel anlamlılık derecesi $\alpha$ belirlenir.
$F$ değeri hesaplanır.
$F$ ile $\alpha$ karşılaştırılır ve $H_0$ bağlamında yorumlanır.
$F$-değeri aşağıdaki şekilde hesap edilir. Burada $S_{r_i}$, i. uyumlamadan (modelden) artık kareler toplamını; $m_i$, i. modelde kullanılan parametre sayısını, $n$ ise gözlemsel veri sayısını göstermektedir.
$$ F = \frac{(\frac{S_{r_1} - S_{r_2}}{m_2 - m_1})}{(\frac{S_{r_2}}{n - m_2})} $$
Eğer modellerin parametre sayıları aynı ($m_1 = m_2$) ise F değeri aşağıdaki şekilde hesaplanabilir.
$$ F = \frac{S_{r_1}}{S_{r_2}} $$
Artık kareler toplamlarının serbestlik derecelerine oranı bir $\chi^2$-dağılımı, $\chi^2$-dağılımlarının oranları ise bir F-dağılımı verir.
Böylece hesaplanan F değeri, bir F dağılımı üzerinde istatiksel anlamlılık derecesi ($\alpha$) ile ifade edilen kritik bir olasılık değeri ile karşılaştırılabilir bir değer olur.
F-testi için, karşılaştırılan artık ya da varyansların normal dağıldığı ve örneklemlerin birbirlerinden bağımsız olduğu varsayılır.
F dağılımı, bir olayın modellenmesinde kullanılan rakip modelin hangi güvenilirlik aralığı limitleri dahilinde birbirleri ile aynı başarımda olduğunun hesaplanmasında kullanılabilidiği gibi iki dağılımın karşılaştırılması temelindeki F-testinde, yani aynı ana dağılımdan seçilen iki örnek dağılımın standart sapmalarının belirli bir güven aralığı dahilinde karşılaştırılmasında da kullanılır. Bu testlerde kullanılan F-test istatistikleri F-dağılımına sahiptir. $U_i$ değerleri birer $\chi^2$-dağılımı ve $df_i$ değeleri bu dağılımların serbestlik dereceleri, $s_i^2$ değerleri bir modelden ya da ortalamadan fark kare toplam değerlerinin serbestlik derecesine bölümü, $\sigma_i^2$ ise ilgili normal süreçlerin standart sapmalarıd olmak üzere u test istatistikleri;
$$ X = \frac{U_1 / df_1}{U_2 / df_2} $$
$$ X = \frac{s_1^2 / \sigma_1^2}{s_2^2 / \sigma_2^2} $$
$$ s_1^2 = \frac{s_1^2}{df_1}, s_2^2 = \frac{s_2^2}{df_2} $$
$$ F(x; df_1, df_2) = \frac{\sqrt{\frac{(df_1~x)^{df_2}~df_2^{df_2}}{(df_1~x + df_2)^{df_1 + df_2}}}}{x~B(\frac{df_1}{2},\frac{df_2}{2})} = \frac{1}{B(\frac{df_1}{2},\frac{df_2}{2})} (\frac{df_1}{df_2})^{\frac{df_1}{2}} x^{\frac{df_1}{2}-1}(1 + \frac{df_1}{df_2}~x)^{-\frac{df_1 + df_2}{2}}$$
Örnek: F Testi¶
Aşağıdaki tabloda birinci sütundaki zaman değerlerine karşılık ikinci sütunda bir niceliğin ölçüm değerleri görülmektedir.
| $t_i$ | $y_i$ |
|---|---|
| 2.5 | 10.07 |
| 2.8 | 11.07 |
| 5.4 | 14.59 |
| 6.5 | 20.75 |
| 9.2 | 27.94 |
| 9.5 | 31.50 |
| 11.0 | 38.04 |
| 13.3 | 49.90 |
| 14.6 | 60.72 |
| 16.4 | 75.57 |
Bu veri setini biri üstel ($y = A~e^{\alpha~t}$), diğeri kuvvet yasasına uyan $y = B~t^{\beta}$ iki ayrı modelle modelleyelim. Bu modellerden hangisinin veriyi daha iyi temsil edebildiğini doğrudan uyumlama eğrilerinden ya da artıklardan anlamamız her zaman mümkün olmaz. Uyumlamanın iyiliğini ölçebilecek istatistiksel yöntemlerin bazıları anlamlı bir seçim yapmaya yeterli olmayabilir. Örneğin en küçük kareler yöntemi ile bu veri setine yapılacak iki uyumlama için lineer korelasyon katsayılarını elde edelim.
Uyumlamak istediğimiz her iki fonksiyon da parametreleri bağlamında lineerleştirilebilir olduğu için lineer en küçük kareler yöntemlerini kullanabiliriz. Öncelikle fonksiyonları doğrusal hale getirelim.
$$ y_1 = A~e^{\alpha~t} \Rightarrow ln(y_1) = ln(A) + \alpha~t $$ $$ y_2 = B~t^{\beta} \Rightarrow ln(y_2) = ln(B) + \beta~ln(t) $$
Yukarıdaki fonksiyonları veriye uyumlamak için basit olması açısından numpy.polyfit() fonksiyonunu deg = 1 parametresiyle doğrusal bir model yapmak için kullanabiliriz. Ancak veri setimizi yukarıdaki lineerleştirmeye uygun olacak şekilde düzenlemeliyiz.
import numpy as np
t = np.array([2.5, 2.8, 5.4, 6.5, 9.2, 9.5, 11.0, 13.3, 14.6, 16.4])
y = np.array([10.07,11.07,14.59,20.70,27.94,31.50,38.04,49.90,60.72,75.57])
lnt = np.log(t)
lny = np.log(y)
ks_y1 = np.polyfit(t,lny,deg=1)
ks_y2 = np.polyfit(lnt,lny,deg=1)
# Katsayilari donusturmeliyiz
A = np.exp(ks_y1[1])
alpha = ks_y1[0]
B = np.exp(ks_y2[1])
beta = ks_y2[0]
# Fit fonksiyonlarini ekrana yazdiralim
print("y1 = {:.2f} e^({:.2f}t)".format(A,alpha))
print("y2 = {:.2f} t^{:.2f}".format(B,beta))
Fit fonksiyonlarını elde ettikten sonra bunları verimiz üzerine çizdirelim.
from matplotlib import pyplot as plt
%matplotlib inline
# Fit fonksiyonlarini biraz fazla nokta
# kullanarak cizdirelim
tyeni = np.linspace(t[0],t[-1])
y1 = A*np.exp(alpha*tyeni)
y2 = B*tyeni**(beta)
plt.plot(t,y,'ro',label="veri")
plt.plot(tyeni,y1,'b-',label="Ustel model")
plt.plot(tyeni,y2,'g-',label="Kuvvet yasasi")
plt.xlabel("t")
plt.ylabel("y")
plt.legend(loc="upper left")
plt.show()
Her ne kadar üstel model (mavi sürekli eğri) veriyi daha iyi temsil ediyor görünse de bir istatistik kullanmamızda fayda var. Örneğin lineer korelasyon katsayısını bu amaçla kullanabiliriz. Ders 4. Eğri Uyumlama dersinde yazdığımız korelasyon katsayısı fonksiyonunu bu amaçla kullanabiliriz. Yalnız fonksiyonu daha sonra yapacağımız F-testinde kullanmak üzere artıkları da döndürecek şekilde değiştirelim.
import numpy as np
def linkor_katsayi(x,y,ymod):
n = y.shape
St = np.sum((y - y.mean())**2)
Sr = np.sum((y - ymod)**2)
r2 = np.abs((St - Sr) / St)
return(r2,Sr)
# Gozlemsel zaman degerleri icin modeller ne ongoruyor?
ymod1 = A*np.exp(alpha*t)
ymod2 = B*t**(beta)
r2_ustel,Sr_ustel = linkor_katsayi(t,y,ymod1)
r2_kuvvet,Sr_kuvvet = linkor_katsayi(t,y,ymod2)
print("Ustel model icin lineer korelasyon katsayisi r^2 = {:.3f}".format(r2_ustel))
print("Kuvvet yasasi icin lineer korelasyon katsayisi r^2 = {:.3f}".format(r2_kuvvet))
Görüldüğü gibi üstel model gerçekten de $1$'e daha yakın bir lineer korelasyon katsayısı vermiştir. Yine de bu modelleri karşılaştırmak için bir de F-testine başvuralım. Çünkü lineer korelasyon katsayıları birbirine yakındır. Bunun için yine bir hipotez testi yapacağız. Bu amaçla sıfır ve alternatif hipotezlerimizi probleme uyun şekilde oluşturalım.
- Sıfır hipotezi ($H_0$): "1. model (üstel), istatistiksel olarak 2. modele (kuvvet yasası) göre daha iyi bir uyumlamadır".
- Alternatif hipotez ($H_1$): "1. model 2. modele göre istatistiksel olarak daha iyi bir uyumlama değildir”
Bu örnek için F-istatistiğini 1. modelden artıkları bölen, 2. modelden artıkları bölünen olarak almak suretiyle hesaplayalım. Test ettiğimiz sıfır hipotezi ($H_0$) 1. model daha iyi bir uyumlamadır şeklinde ifade edildiğine göre ondan artıkların daha küçük olmasını bekleriz. Dolayısı ile artıkları büyük olan 2. modeli bölümün üstüne almamız gerekir.
F_hesap = Sr_ustel/ Sr_kuvvet
print("Hesaplanan F-degeri = {:.3f}".format(F_hesap))
Yine bu örnek için her iki model de ikişer parametre ($(A,\alpha)$, $(B,\beta)$) ve aynı sayıda ölçüm ($n = 10$) içerdiğinden serbestlik dereceleri $df_1 = df_2 = n – 2 = 8$'dir.
Bir sonraki adım F-tablolarından verilen istatistik anlam seviyesi için, önreğin $\alpha = 0.01$ tek taraflı bir test için olanının seçilerek $df_1 = df_2 = 8$ serbestlik derecesi için verilen F-istatistiğine bakmaktır. Daha önce yaptığımız gibi bunun için scipy.stats.f.ppf() fonksiyonunu kullanacağız.
from scipy import stats as st
n = 10
df1 = df2 = n-2
alfa = 0.01
F_degeri = st.f.ppf(1-alfa, df1, df2)
print("Tablodan elde edilen F-degeri: {:.3f}".format(F_degeri))
Bu durumda $F_{hesap} = 0.089 < F_{tablo} = 6.029$ olduğundan $H_0$ hipotezi reddedilemez. %99 güven düzeyinde üstel modelin, kuvvet yasasıyla verilen modele üstün olduğu görülmektedir. Güven düzeyi %90’a ($\alpha = 0.10$) düşürülse dahi ($F = 2.5893$) $H_0$ reddedilemez.
Şimdi sıfır hipotezini ($H_0$): "2. model (kuvvet yasası), 1. modele (üstel) göre istatistiksel olarak daha iyi bir uyumlamadır", alternatif hipotezi ($H_1$) ise "2. model 1. modele göre istatistiksel olarak daha iyi bir uyumlama değildir" şekilnde değiştirelim ve sonuçta bir değişiklik olup olmayacağına bakalım.
Bunun için öncelikle F-istatistiğini fark kareler toplamının bölümde yerini değiştirerek hesaplayalım.
F_hesap = Sr_kuvvet / Sr_ustel
print("Hesaplanan F-degeri = {:.3f}".format(F_hesap))
Yine aynı şekilde her iki model de ikişer parametre ($(A,\alpha)$, $(B,\beta)$) ve aynı sayıda ölçüm ($n = 10$) içerdiğinden serbestlik dereceleri $df_1 = df_2 = n – 2 = 8$'dir. Bu serbestlik derecesi için verilen F-istatistiğini scipy.stats.f.ppf() fonksiyonunu kullanarak hesaplamıştık. $F_{hesap} = 11.212 > F_{tablo} = 6.029$ olduğundan $H_0$ bu kez reddedilir. Ancak bu kez reddedilen sıfır hipotezi kuvvet yasasının üstel modele göre veriye daha uygun olduğunu ifade etmektedir. Bu durum alternatif hipotez için, yani kuvvet yasasının üstel modele göre veriye daha uygun olmadığına yönelik görüşe kanıt teşkil eder.
Durbin Watson İstatistiği İle Model Başarımı Testi¶
Durbin-Watson istatistiği, bir regresyon analizinden artıklarda birinci dereceden gecikme grafiğinde otokorelasyonun varlığını tespit etmek için kullanılan bir test istatistiğidir. Durbin-Watson istatistiği aşağıdaki şekilde tanımlanır ve bu istatistik hesaplanıp, kullanılarak model uyumlamanın başarısı test edilebilir.
$$ \mathcal{D} = \frac{\sum_{i=2}^{N} [R_{i} - R_{i-1}]^2}{\sum_{i=1}^{N} [R_{i}]^2} $$
Burada $\mathcal{D}$, Durbin-Watson istatistiği, $R_i$, uyumlamadan (fitten) artıklar, $R_{i-1}$ ise artıkların $i-1$'den başlayarak dizilmesiyle oluşan artıklardır. $\mathcal{D}$’nin alacağı değerlere göre aşağıdaki çıkarımlar yapılabilir:
- $\mathcal{D} = 0$: artıklar sistematik olarak korelasyon göstermektedir.
- $\mathcal{D} = 2$: artıklar normal dağılım göstermektedir.
- $\mathcal{D} = 4$: artıklar sistematik olarak antikorelasyona sahiptir.
Örnek olarak Durbin-Watson istatistiğinin yanı sıra pek çok başka model istatistiği de veren bir paket olan statsmodel paketi fonksiyonlarını modelleme için kullanalım. Durbin-Watson istatistiği yukarıdaki formülden kolayca hesaplanabilen bir istatistiktir. Dolayısıyla bu amaçla bir fonksiyon yazabilir ve modelleme işlemlerinizde scipy ya da numpy fonksiyonlarının yanı sıra bu yazdığınız fonksiyonu da kullanabilirsiniz. Ancak, burada hem bu fonksiyonu yazmayı size alıştırma olarak bırakmak hem de gelişmiş bir istatistiksel model paketini tanıtarak örneklemek açısından statsmodel ile bir modelleme örneği yapılacaktır.
Öncelikle yine modellemek istediğimiz verimizi oluşturalım ve grafiğini çizdirelim.
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
np.random.seed(48)
x = np.linspace(0,10,100)
y = -3*x + 5 + np.random.rand(100)*10
plt.plot(x,y,'ro')
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Şimdi bu değişimi modellemek üzere statsmodel.api fonksiyonlarından en küçük kareler yöntemiyle model uyumlaması yapan OLS() (ing. Ordinary Least Squares) fonksiyonunu kullanalım. Öncelikle uyumlanacak doğrunun y-ekseninde bir kayma terimi (sabit) de içerdiğini belirtmek ve bu terimi de uyumlamak için statsmodel.api.add_constant fonksiyonu kullanılmaktadır. Bu terim bağımsız parametreye 1'lerden oluşan bir sabit terim ekler ve fit sırasında bu parametre de uyumlanır ($y = m*x + n*1$). Bu nedenle
X = sm.add_constant(x)
ifadesiyle bir $X$ iki boyutlu dizisi oluşturulmuştur. Bu ifade bağımsız değişken $x$'i $X = [x,1] \rightarrow y = m*x + n*1$ şeklinde ifade etmek ve uyumlama sırasında $m$ ve $n$ parametrelerinin bulunmasını sağlamaktadır.
OLS fonksiyonuna önce argüman olarak verinin bağımlı değişkenini içeren $y$ dizisi, sonra bağımsız değişken $x$ ve sabit parametreyi (1) içeren iki boyutlu $X$ değerleri numpy dizisi olarak sağlanmaktadır. Daha sonra uygulanan statsmodel.api.OLS() fonksiyonu bir model nesnesi döndürmekte bu model nesnesi aşağıdaki kod parçasında $model$ değişkenine atanmaktadır. Eğri uyumlamada gördüğünüz LMFIT fonksiyonunda olduğu gibi model nesnesinin uyumlama için kullanılan bir fit() fonksiyonu bulunmakta; model sonuçları ise bir sonuç nesnesi ile birlikte bu fonksiyonla döndürülmektedir. Bu nesnenin summary() metodu aralarında Durbin-Watson istatistiğinin de bulunduğu model istatistiklerini görüntülemek için kullanılmaktadır.
import statsmodels.api as sm
X = sm.add_constant(x)
model = sm.OLS(y, X)
sonuclar = model.fit()
print(sonuclar.summary())
İstendiği takdirde modele ilişkin istatistiklerin bazıları sonuç nesnesinin metodları ile ayrı ayrı da görüntülenbilir. Yapılması gereken yukarıdaki istenen istatistiği veren metodu sonuç nesnesi üzerine uygulamaktır. Örneğin model parametreleri params, bu parametrelerin hataları bse ve korelasyon katsayısı ($r^2$) rsquared metodları ile aşağıdaki şekilde görüntülenebilir.
print('Parametreler: ', sonuclar.params)
print('Standart hatalari: ', sonuclar.bse)
print('r^2: ', sonuclar.rsquared)
İstendiği takdirde uyumlanan fonksiyon yine bu nense üzerinden grafiğe aktarılabilir.
plt.plot(x,y,'ro', label='veri')
plt.plot(x,sonuclar.params[0]+sonuclar.params[1]*x, 'b-', label='model')
plt.xlabel("x")
plt.ylabel("y")
plt.legend(loc="best")
plt.show()
Oldukça gelişmiş model istatistikleri ve olanaklar sunan bu pakete daha ileri model örnekleri için tekrar döneceğiz. Ancak şu aşamada yaptığımız uyumlamanın Durbin-Watson İstatistiği'ne bakalım. $\mathcal{D} = 2.058~\sim~2$ bulunmuş olması artıkların normal dağılıma sahip olduğu iyi bir modele ulaşıldığını göstermektedir.
Akaike Bilgi Kriteri¶
Akaike Bilgi Kriteri (ing. Akaika Information Criterion,AIC), bir modelin başarısını veren ölçütlerden biridir. Doğrudan bir modelin başarısını verebildiği gibi, farklı modellerin karşılaştırılabilmesinde de olanak sağlar. Bir sıfır hipotezinin testine dayanan bir değerlendirme yapılmadığı için bir modelin doğruluğunun hangi istatistiksel anlamlılık seviyesinde reddedilip edilmeyeceğiyle ilgili bir bilgi vermez. Dolayısıyla ilgilenilen modellerin tamamının veriye uyum sağlayıp sağlamamalarının bir göstergesi olmayıp, sadece bir karşılaştırmalarıdır.
Akaike Bilgi Kriteri (AIC) aşağıdaki şekilde tanımlanır ve $AIC$ değeri küçük olan modelin başarısı daha yüksektir.
$$ AIC = 2~k - 2~ln(\hat{L}) $$
Burada $k$, modelde kullanılan parametrelerin sayısıdır. $\hat{L}$ ise modelin doğruluğu durumunda gözlem verisinin olabilirlik fonksiyonudur (ing. likelihood function).
Olabilirlik fonksiyonu, bir modelin önerdiği olasılık dağılımına göre tüm verilerin bu model ile oluşabilme olasılıklarının çarpımıdır. $AIC$ ya da bir çok hesapta olabilirlik fonksiyonunun kendisi kullanılabildiği gibi, kolay hesaplanabilirliği sebebiyle logaritması da kullanılmaktadır.
Görüldüğü gibi parametre sayısının artması durmunda $AIC$ değeri artar. Yani modelde kullanılan parametre sayısının fazlalığı, modelin gerçek değişimden çok gürültü modellemeye doğru eğilimde bulunacağı kabulu yapılır. Dolayısıyla daha az parametre sayısına sahip (daha basit) modellerin seçilmesi yönünde bir denge sağlamış olur.
Eğer uyumlama yapılmak istenen veri sayısı çok az ise $AIC$ daha fazla parametreye sahip olan modelin seçilmesine sebep olur. Bu sebep ile AIC düzeltmesi (AIC correction, $AICc$) tanımı yapılmıştır.
$$ AICc = AIC + \frac{2k~(k+1)}{n - k -1} $$
Düzeltme terimindeki $n$, gözlem verisi sayısı; $k$ ise parametre sayısıdır.
$AIC$ ya da $AICc$ kullanımı için sınır koşul olarak $n/k > 40$ kullanılır. Yani gözlem verisi sayısı, modeldeki parametre sayısından en az 40 kat büyük ise $AIC$ kullanılabilir. Bu durumda parametre sayısı fazlalığının etkisi ihmal edilebilecek kadar küçük kalır. Eğer bu oran 40’tan daha küçük ise mutlaka $AICc$ kullanılmalıdır.
$AIC$ ya da $AICc$ değeri pozitif ya da negatif olabilir. Başarılı model her zaman daha küçük değere sahip olandır. $AIC$, farklı sayıda veriye sahip uyumlama işlemlerinde kullanılamaz. $AIC$, bir sıfır hipotezi sınaması olmadığı için sonuçta anlamlılık düzeyi, hipotez reddi gibi ifadeler kullanılmamalıdır.
Model uyumlamada AIC hesabı aşağıdaki şekilde yapılabilir.
$$ AIC = n~ln~(\frac{S_r}{n}) + 2k $$
Daha önce de tanımlandığı gibi $S_r$ gözlemsel verinin modelden fark kareler toplamını ifade eder. $n$ gözlem verisi sayısı, $k$ modelin parametre sayısıdır.
Model seçimi durumunda:
- Modellerin tamamının $AIC$ değerleri yukarıdaki ifadeden hesaplanır.
- En düşük $AIC$ değerine ($AIC_{min}$) sahip olan model baz alınarak $AIC$ farkları ($Delta_i$) hesaplanır.
- Her bir modelin olabilirlik fonksiyonu değeri, o modelin $AIC$ değerinin minimum AIC değerine sahip modelden ($AIC_{min}$) farknın üstel bir ifadesi ile orantılıdır.
- Ardından modellerin Akaike ağırlıkları ($w_i$) hesaplanır. Bu değer, her bir model için bulunan olabilirlik fonksiyonunun karşılaştırma yapılan tüm modellerin toplam olabilirliklerine oranıdır.
$$ \Delta_i = AIC_i - AIC_{min} $$ olmak üzere,
$$ \mathcal{L}(g_i~|~y)~\alpha~e^{-\frac{1}{2} \Delta_i} $$
$$ w_i = \frac{e^{-\frac{1}{2} \Delta_i}}{\sum_{r=1}^r e^{-\frac{1}{2} \Delta_r}} $$
Örnek: Akaike Bilgi Kriteriyle Model Karşılaştırması¶
F-Testi ile model karşılaştırması için yaptığımız örneğe geri dönelim. Örneğimizin verisi aşağıdaki tabloda birinci sütunda verilen zaman değerlerine karşılık ikinci sütundaki bir niceliğin ölçüm değerlerinden oluşmaktaydı.
| $t_i$ | $y_i$ |
|---|---|
| 2.5 | 10.07 |
| 2.8 | 11.07 |
| 5.4 | 14.59 |
| 6.5 | 20.75 |
| 9.2 | 27.94 |
| 9.5 | 31.50 |
| 11.0 | 38.04 |
| 13.3 | 49.90 |
| 14.6 | 60.72 |
| 16.4 | 75.57 |
Bu veri setini biri üstel ($y = A~e^{\alpha~t}$), diğeri kuvvet yasasına uyan $y = B~t^{\beta}$ iki ayrı modelle modellemiş ve bu iki modeli F-testi ile karşılaştırmıştık. Bu kez bir de doğru modeli tanımlayalım ve bu modellerden hangisinin veriyi daha iyi temsil edebildiğini Akaike Bilgi Kriteri ($AIC$) kullanarak belirleyelim.
Uyumlamak istediğimiz her iki fonksiyon parametreleri bağlamında lineerleştirilebilirdi. Eklediğimiz son model de zaten lineer bir model olsun.
$$ y_1 = a_0 + a_1*t $$ $$ y_2 = A~e^{\alpha~t} \Rightarrow ln(y_2) = ln(A) + \alpha~t $$ $$ y_3 = B~t^{\beta} \Rightarrow ln(y_3) = ln(B) + \beta~ln(t) $$
Yine numpy.polyfit() gibi basit bir fonksiyon kullanarak bu modelleri verimize uyumlayabiliriz ve bu modellerden verimizin fark kareleri toplamı, veri sayısı ve parametre sayısını dikkate alarak Akaike Bilgi Kriteri'ni ($AIC$), kritierin veri sayımız az olduğu için düzeltilmiş değerini ($AICc$) kendi yazabileceğiniz bir Python fonksiyonuyla hesaplayıp, model karşılaştırması yapabiliriz. Yine önce verimizi tanımlayalım ve uyumlamalarımızı yapıp veri üzerinde görerek başlayalım.
import numpy as np
t = np.array([2.5, 2.8, 5.4, 6.5, 9.2, 9.5, 11.0, 13.3, 14.6, 16.4])
y = np.array([10.07,11.07,14.59,20.70,27.94,31.50,38.04,49.90,60.72,75.57])
lnt = np.log(t)
lny = np.log(y)
ks_y1 = np.polyfit(t,y,deg=1)
ks_y2 = np.polyfit(t,lny,deg=1)
ks_y3 = np.polyfit(lnt,lny,deg=1)
# Katsayilari donusturmeliyiz
a0 = ks_y1[1]
a1 = ks_y1[0]
A = np.exp(ks_y2[1])
alpha = ks_y2[0]
B = np.exp(ks_y3[1])
beta = ks_y3[0]
# Fit fonksiyonlarini ekrana yazdiralim
print("y1 = {:.2f} + {:.2f}t)".format(a0,a1))
print("y2 = {:.2f} e^({:.2f}t)".format(A,alpha))
print("y3 = {:.2f} t^{:.2f}".format(B,beta))
Fit fonksiyonlarını elde ettikten sonra bunları verimiz üzerine çizdirelim.
from matplotlib import pyplot as plt
%matplotlib inline
# Fit fonksiyonlarini biraz fazla nokta
# kullanarak cizdirelim
tyeni = np.linspace(t[0],t[-1])
y1 = a0 + a1*tyeni
y2 = A*np.exp(alpha*tyeni)
y3 = B*tyeni**(beta)
plt.plot(t,y,'ro',label="veri")
plt.plot(tyeni,y1,'m-',label="Dogrusal model")
plt.plot(tyeni,y2,'b-',label="Ustel model")
plt.plot(tyeni,y3,'g-',label="Kuvvet yasasi")
plt.xlabel("t")
plt.ylabel("y")
plt.legend(loc="upper left")
plt.show()
Şimdi modellerimiz için $AIC$ ve $AICc$ değerlerini hesaplamak üzere basit bir fonksiyon yazalım ve verilerimiz ile modelimizi bu fonksiyona gönderereek $AIC$ ve $AICc$ değerlerini hesaplayalım
import numpy as np
def aic(y,ymodel,katsayilar):
Sr = np.sum((y - ymodel)**2)
n = len(y)
k = len(katsayilar)
AIC = n*np.log(Sr / n) + 2*k
AICc = AIC + (2*k*(k+1)) / (n - k - 1)
return(AIC,AICc,Sr)
AIC_1,AICc_1,Sr1 = aic(y,a0+a1*t,(a0,a1))
AIC_2,AICc_2,Sr2 = aic(y,A*np.exp(alpha*t),(A,alpha))
AIC_3,AICc_3,Sr3 = aic(y,B*t**beta,(B,beta))
print("Dogrusal model icin AIC={:.2f}, AICc = {:.2f}".\
format(AIC_1, AICc_1))
print("Ustel model icin AIC={:.2f}, AICc = {:.2f}".\
format(AIC_2, AICc_2))
print("Kuvvet yasasi modeli icin AIC={:.2f}, AICc = {:.2f}".\
format(AIC_3, AICc_3))
Her üç modelin parametre sayısı ve veri noktası sayısı da aynı olduğu için karşılaştırma doğrudan düzeltilmiş Akaika Bilgi Kriteri ($AICc$) üzerinden de yapılabilir. Böyle bir karşılaştırma üstel modelin diğer iki modele göre üstün olduğunu göstermektedir. Ancak öğretici olması açısından yukarıdaki şemaya bağlı kalarak bir tablo hazırlayyalım. Bunun için yukarıda tanımlanan $\Delta_i$ ve $w_i$ değerlerini de hesaplamamız gerekir. Bu hesaplar için gerekli bütün parametrelere sahibiz.
AICc = np.array([AICc_1,AICc_2,AICc_3])
AICc_min = np.min(AICc)
Delta = np.array([aicc-AICc_min for aicc in AICc])
expDelta = np.exp(-0.5*Delta)
wi = expDelta / np.sum(expDelta)
print("Sr_i: ", (Sr1,Sr2,Sr3))
print("Delta_i: ", Delta)
print("w_i: ", wi)
Bu değerleri bir tablo şeklinde ifade edelim
| Model | k | Sr | AICc | $\Delta_i$ | $w_i$ |
|---|---|---|---|---|---|
| Dogrusal | 2 | 242.68 | 37.61 | 20.33 | 0.000039 |
| Ustel | 2 | 31.79 | 17.28 | 0.00 | 0.999558 |
| Kuvvet | 2 | 356.44 | 41.45 | 24.17 | 0.000006 |
Bu değerlerden hareketle artık modelleri birbirleriyle karşılaştırmak mümkündür. Bunun için model ağırlıklarını birbirlerine bölerek, bir modelin diğerine göre "olabilirliğinin" ($\mathcal{L}$) ne kadar daha yüksek (olası) olduğunu belirlemeliyiz.
L2_1 = wi[1] / wi[0]
L2_3 = wi[1] / wi[2]
L1_3 = wi[0] / wi[2]
print("Ustel model dogrusal modele gore {:.2f} kat daha yuksek olabilirlige sahiptir".\
format(L2_1))
print("Ustel model kuvvet yasasi modeline gore {:.2f} kat daha yuksek olabilirlige sahiptir".\
format(L2_3))
print("Dogrusal model kuvvet yasasi modeline gore {:.2f} kat daha yuksek olabilirlige sahiptir".\
format(L1_3))
Sonuç olarak bu üç modelden üstel modelin eldeki gözlemsel veriyi modellemek konusunda başarımı en yüksek model olduğu tüm ölçütlerde görülmektedir. Aslında bu anlamda fark-karelere bakmak bile yaterli olabilirdi. Ancak buradaki üçüncü karşılaştırma anlamlıdır. Veriyi uyumlamak anlamında başarımı yakın iki modelden doğrusal modelin kuvvet yasası modeline göre 6.83 kat daha olası olduğu Akaike Bilgi Kritieri kullanılarak ortaya konmuş durumdadır. $AIC$ statsmodels ve LMFIT gibi gelişmiş model paketleri tarafından da verilen bir çıktı parametresidir. Yukarıdaki örneği bu modüllerden en az birini kullanarak çözmeyi deneyiniz. Bu modüllerin asıl avantajının parametreleri açısından örneğimizdeki gibi lineer modeller değil lineer-olmayan (nonlinear) söz konusu olduğunda ortaya çıktığını belirtmekte yarar vardır.
Bayesian Bilgi Kriteri¶
Bayesian Bilgi Kriteri (ing. Bayesian Information Criterion, BIC) de bir modelin başarısını veren bir başka bir ölçüttür. $AIC$ ile yakından ilişkilidir. Ancak ilave parametrelere daha hassastır ve temel olarak aşağıdaki şekilde tanımlanır.
$$ BIC = ln~(n)~k - 2~ln~(\mathcal{L}) $$
Burada $n$, gözlem sayısı, $k$, modelde kullanılan parametrelerin sayısıdır. $L$ ise modelin doğruluğu durumunda gözlem verisinin olabilirlik fonksiyonudur (ing. likelihood function).
Tıpkı $AIC$ 'de olduğu gibi $BIC$ değerinin küçük olduğu modeller veriyi uyumlamada daha başarılıdır denilebilir. Eğri uyumlama işleminde kriter, aşağıdaki şekilde kullanılabilir. Burada $S_r$, yine modelden fark kareler toplamını ifade etmektedir.
$$ BIC = n~ln~(\frac{S_r}{n}) + k~ln~(n) $$
Farklı modellerin $BIC$ değerleri arasındaki fark $\Delta_{BIC}$ aşağıdaki şekilde yorumlanabilir.
- $0 < \Delta_{BIC} < 2$, yüksek $BIC$ değerli modele karşı güçlü bir kanıt yoktur.
- $2 < \Delta_{BIC} < 6$, yüksek $BIC$ değerli modele karşı "pozitif" bir kanıt vardır.
- $6 < \Delta_{BIC} < 10$, yüksek $BIC$ değerli modele karşı "güçlü" bir kanıt vardır.
- $\Delta_{BIC} > 10$, yüksek $BIC$ değerli modele karşı "çok güçlü" bir kanıt vardır.
Kendiniz bir $BIC$ hesaplama fonksiyonu yazabileceğiniz gibi, statsmodels ve LMFIT gibi güçlü modülleri kullanarak $BIC$ değerlerini elde edebilir ve yukarıdaki şemayı kullanarak karşılaştırma da yapabilirsiniz.
Kaynaklar¶
- Measurements and their Uncertainties, Ifan G. Hughes & Thomas P.A. Hase, Oxford University Press, 2010
- Data Reduction and Error Analysis for the Physical Sciences, Philip R. Bevington & D. Keith Robinson, MC Graw Hill, 2003
- scipy.stats Modülü Dokümantasyonu
- numpy.random Modülü Dokümantasyonu
- LMFIT Dokümantasyonu
- statsmodels Dokümantasyonu