800100715151 Astronomide Veritabanları

Ders - 07 İstatistiksel Yöntemler, Tahmin ve Çıkarım–II

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

Monte Carlo Yöntemleri

Giriş

Şu ana kadar değişkeni tek bir olasılık dağılmından değer alan problemler üzerinde çalıştık. Ancak fen bilimlerinde karşılaşılan problemler genellikle birden fazla değişken içerebildiği için, çözümleri de birden fazla olasılık dağılımının kombinasyonunu gerektirir. Söz konusu problemlerde değişkenlerin alabileceği değerlerin belirli aralıkları için olasılıkları, bu dağılımların bu aralıklardaki integrasyonu (ya da toplamı) ile elde edilir.

Monte Carlo Yöntemleri, değişkenlerin olasılık dağılımlarından rastgele örneklemler üretip, bunların kombine olasılıklarını belirlemek için çoklu integrasyonlar yapabilmeye olanak sağlayan tekniklere verilen isimdir. Tekrarlanan rastgele örneklemlerle belirlenen değişken olasılıklaırnın kombinasyonu, nümerik integrasyonla elde edilir.

Örneklemin rastgeleliğine dayandıkları için bu yöntemler Fransa’nın kumarhaneleriyle meşhur şehri Monte Carlo’nun adıyla bilinirler.

Markov Zinciri

Her ne kadar $\pi$ sayısının hesabı için verilen örnek oldukça basitse de Monte Carlo yöntemlerinin dayandığı temel mantığı anlamak açısından açıklayıcıdır. Temel mantık bir olasılık dağılımından (burada tekdüze dağılım) rastgele örneklemler üretip, bu örneklem için olasılık hesabı yapmaktan ibarettir.

Markov özelliği (ing. Markov property) herhangi bir rastgele durumun "hafızasızlık" (ing. memorylessness) özelliği taşıması, yani kendisinden önceki durumlardan bağımsız olması demektir.

Markov zinciri (ing. Markov chain) ise, Markov özelliği taşıyan durumlar arasındaki olasılık temelli süreksiz geçişler sürecine verilen isimdir.

Örneğin, zar oyunları birer Markov zinciridir. Ancak kart oyunları, oynadıkça destedeki kartların sayıları ya da türleri değiştiği için hafızaya sahiptirler ve Markov özelliğine sahip değildirler. Bu sebep ile bir Markov zinciri oluşturmazlar.

Rastgele Yürüyüş

Rastgele yürüyüş, bir matematiksel uzayda, rastgele adımlar atarak izlenilen yolu açıklayan rastgele bir süreçtir. Örnek olarak, bir akışkandaki parçacığın yolu, yemek arayan bir böceğin yolu, bir fotonun yılıdız içindeki yolculuğu rastgele yürüyüş yaklaşımı ile incelenebilir. Rastgele yürüyüş süreçleri Markov özelliğini taşıdığı sürece birer Markov zinciridir.

Rastgele Yuruyus

Monte Carlo Markov Zinciri MCMC

Bir olasılık dağılımı üzerinden örneklemler ile, Markov zincirinin denge dağılımına ulaşmayı hedefleyen algoritmalar grubuna verilen isimdir. Fizik, biyoloji, ekonomi, dilbilim gibi birçok farklı alanda kullanılmaktadır. Örneklenmek istenen soruna ilişkin boyut sayısının yüksek olması, diğer yöntemleri kullanışsız hale getirirken (bkz. ing. curse of dimensionality) MCMC yöntemleri çoğu zaman uygulanabilir tek yöntem olmaktadır. Genellikle rastgele yürüyüş mantığı ile çalışan algoritmalar kullanır. Bunlardan bazıları:

  • Metropolis-Hastings
  • Gibbs
  • Slice
  • Reversible-jump

Metropolis Hastings Algoritması

Bu algoritma, belirli bir olasılık dağılımından elde edilen rastgele örneklerin, doğrudan (ana) dağılımın kestirilmesinde ya da integral hesaplamasında kullanılmasına dayanmaktadır.

Ardarda üretilen rastgele değerler, bir olasılık karşılaştırmasına tabi tutularak seçilir. Bu şekilde Markov zincirinin denge dağılımına erişilene kadar yeni rastgele değerler üretilmeye devam edilir.

Tekrar sayısı ne kadar fazla olursa, denge dağılımını örnekleme kalitesi o oranda artar.

Algoritma aşağıdaki adımlar üzerine kuruludur:

1) Rastgele başlangıç parametre seti ($x_0$) oluşturulur. 2) Varsayılan bir olasılık dağılımından bir diğer parametre seti ($x_i$) seçilir. 3) İki setten hangisinin kabul edileceği aşağıdaki şekilde belirlenir: Eğer random(0,1) <= min(1, P(xi)/P(xi+1)) ise $x_i$ seti seçilir. Değilse; $x_{i+1}$ seti seçilir. 4) Hesaplama 2. adımdan istenilen sayıda tekrarlanır. 5) Hesap tekrarı (iterasyon) tamamlandıktan sonra seçilen parametrelerin histogramları oluşturulur. 6) Histogramlarda en çok tekrar eden değer aralıkları parametrelerin değer aralığı olarak kabul edilir.

Rastgele üretilen daha olası parametre setinin, $0-1$ arasında tekdüze bir dağılımdan rastgele seçilen bir sayı ile karşılaştırılması, yerel minimuma hapsedilmeyi önleyip global minimumun bulunabilme ihtimalini arttırır; ancak bu ihtimal her zaman vardır. Yerel minimuma hapsolmamak için farklı ilk parametre setlerine doğrudan "zıplamak" gerekebilir.

Bu algoritmada, kaç defa parametre seti oluşturulması gerektiği bilinmemektedir.

Üretilen ilk değerlerin, istenilen dağılımdan çok uzak değerler alması olasıdır. Bu sebep ile ilk tekrarların önemli bir kısmı hesaba katılmadan ayıklanır; bu aralığa "ısınma dönemi" (ing. Burn-in) adı verilir. Isınma dönemine dahil olan tekrar sayısı da önceden bilinmemektedir.

Metropolis-Hastings

Aşağıda bir çift yıldızın dönemi sürekli olarak değiştirilerek bir O-C analizinde en küçük farkı verecek yörünge dönemi aranmaktadır. x-ekseni adım sayısını; y-ekseni gün cisinden yörünge dönemini göstermektedir. Renkli her bir eğri başka bir rastgele yürüyüşü (ing. random walk) göstermekte, birden fazla yürüyücüyle en iyi çözüm aranmaktadır.

Metropolis-Hastings

Eğer üretilen parametre setleri birbirlerinden önemli miktarda farklı değerlere sahip olursa, setlerin karşılaştırılması sonucunda yeni setin reddedilme sıklığı artar. Bu durumda kullanılmayan setler üretilmiş olur.

Eğer parametre setleri birbirlerine çok yakın değerlere sahip olurlarsa, yeni setler büyük ihtimalle kabul edilmekle birlikte, istenilen dağılıma ulaşmak için aşırı set üretimi gerekecektir.

Yukarıdaki iki durumda da işlem süresi uzayacaktır. Buradaki en büyük sorun, yeni setlerin üretileceği dağılımların seçimi için genel bir kural bulunmamasıdır. Yani eldeki soruna özel çözümler üretmek gerekmektedir.

Rastgele üretilen ardışık parametre setleri birbirlerinden tam olarak bağımsız değildirler. Bu sebeple eğer birbirinden bağımsız setlerin incelenmesi ihtiyacı söz konusuysa, üretilen ve kabul edilen tüm setlerin sadece her $n$ tanesinden biri alınmalıdır. Bu durumda diğer setlerin kullanılmaması söz konusu olmaktadır.

Bayesian İstatistik Paradigması

Olasılık kavramına temel yaklaşımında fark bulunan iki önemli istatistik paradigması vardır: Bayesian yaklaşım ("Bayesci") ve Frekansçı (klasik) yaklaşım.

Frekansçı (klasik) yaklaşımda parametre değerleri için bir popülasyon varsayılır. Bu popülasyon uzun vadede, herhangi bir deney ya da olayın sonsuz kez gerçekleştirilmesi sonucu parametrelerin aldığı değerleri temsil eder. Gerçekte böyle bir popülasyon yoktur. Ancak örneğin bir yazı / tura atışında sonsuz kez bu deneyin yapılması sonucu atışların yarısının "Tura", diğer yarısının ise "Yazı" geleceğini varsayarak parametrenin her iki değerini $\%50$ olasılık atarız. Bir başka deyişle, klasik yaklaşımda yazı gelme olasılığı da tura gelme olasılığı da birbiriyle aynı koşullarda sonsuz kez gerçekleştirilmiş bir yazı-tura deneyinde yazı (ya da tura) gelme sıklıklarını değerlendirmektir. Bu durumda atılan her bir yazı-turayla değişen şey veridir, yazı ve tura gelme olasılıkları yani modelin parametreleri ise sabittir! Bu nedenle, frekansçı yaklaşım "objektif" bir yaklaşımdır. Zira deneylerin öncülü (ing. prior) olmadığı gibi, deneyin herhangi bir sonucuna yönelik bir "inanç" da söz konusu değildir. $n$ kez tekrarlanan bir deneyde bir $x$ olayının $n_x$ kez gözlenmesi durumunda olasılığı;

$$ P(x) = \frac{n_x}{n} $$

ile verilir ve gerçek olasılığı ancak deneyin sonsuz kez tekrarlanmasıyla ulaşılabilir ($n \rightarrow \infty$).

1700'lerin ortalarında Thomas Bayes (1701-1761) tarafından geliştirilen Bayesian yaklaşıma göre bir $x$ olayının gerçekleşmesi olasılığı ona duyulan inancın bir ifadesidir. Bu istatistik paradigmasında bir $x$ olayının olasılığı, deney ya da gözlem yapılmadan önce $x$ olayının gözlenmesi olasılığına olan inancımız (ing. prior belief) ve eldeki veri (ing. evidence) dikkate alınarak hesaplanır. Yani veri sabittir, ancak modelin parametreleri değişkendir. Daha çok veri aldıkça olayın gerçekleşme olasılığına ilişkin "inancımızı", yani o olayın olasılığını güncelleriz.

Bayes Formülü

Bayesian istatistiğin temeli aşağıdaki basit formüle dayanır:

$$ P(\theta~|~veri) = \frac{P(veri~|~\theta) \times P(\theta)}{P(veri)} $$

Burada;

$P(veri~|~\theta)$, Olabilirlik fonksiyonu (ing. Likelihood Function) : $\theta$ model olarak seçildiğinde söz konusu veri setini elde etme olasılığı,

$P(\theta)$, Öncül (ing. Prior): Herhangi bir ön bilgi (veri) olmaksızın $\theta$ gibi bir modelin (ya da teorinin) doğru olma olasılığı,

$P(veri)$, Veri (Kanıt, ing. Evidence): Tüm olası $\theta$ model seçimleri için söz konusu veri setini elde etme olasılığı (hesaplanmasının zor olacağı açıktır),

$P(\theta | veri)$, Ardıl (ing. Posterior): Eldeki veriye dayalı olarak $\theta$ modelinin geçerli olma olasılığı olarak tanımlanabilir. Aslında bu parametreler doğrudan olasılıklara karşılık gelmeyebilir, ancak başlangıç için bunları olasılık olarak düşünebilirsiniz.

Bayesian istatistik paradigmasının temel amacı $\theta$’nın bütün değerleri için birer olasılık değeri, bir başka deyişle, tüm $\theta$ modellerinin gerçekleşme olasılıklarını ayrı ayrı hesaplamaktır. Bu nedenle ardıl (posterior) bir olasılık dağılımıdır.

Bayes Formülünün Türetilmesi

Her ne kadar Bayesian istatistik frekansçı istatistikten farklı bir paradigma olsa da Bayes formülü klasik olasılık yaklaşımıyla türetilebilir.

Bir $A$ olayının $B$ olayına bağlı olarak gerçekleşme olasılığı olan koşullu olasılığın $B$ olayının gerçekleşmesi durumunda $A$ olayının gerçekleşmesinin olasılığı olduğunun görmüş ve aşağıdaki şekilde tanımlamıştık:

$$ P(A~|~B) = \frac{P(A,B)}{P(B)} $$

Bu durumda,

$$ P(A,B) = P(A~|~B) \times P(B) $$

Aynı şekilde,

$$ P(B~|~A) = \frac{P(B,A)}{P(A)} \rightarrow P(B,A) = P(B~|~A) \times P(A)$$

Diğer taraftan bileşke olasılıklar birbirine eşit ($P(A,B) = P(B,A)$) olacağından ;

$$ P(A,B) = P(B,A) \rightarrow P(A~|~B) \times P(B) = P(B~|~A) \times P(A) $$

$$ P(A~|~B) = \frac{P(B~|~A) \times P(A)}{P(B)} $$

elde edilir ki, bu Bayes formülüdür.

Kavramsal Örnek: Hileli Madeni Parayla Yazı / Tura Atışı

Madeni Para

Bayesian istatistik paradigmasıyla düşünme, yani deneyler ya da gözlemlerle elde edilen veriye ve arka plan bilgimize dayalı olarak bir olayın olası tüm sonuçlarının olasılıklarını hesaplamaya yönelik iyi bir örnek yazı / tura atış problemidir.

Bir kumarhanede bir seri yazı / tura atışı problemine tanıklık ettiğimizi düşünelim. Öğrenmeye çalıştığımız şey yazı / tura atışını yapan kişinin bir şekilde atışın sonuçlarını belirleyip belirlemediği, eğer belirliyorsa bunu hangi başarımda yaptığı olsun.

Problemi frekansçı (klasik) istatistik yaklaşımı ile ele alırsak madeni paranın hangi miktarda hileli olduğuna ilişkin bir sıfır hipotezimiz ($H_0$) olur (null hypothesis). Elimizdeki veriye (bir seri yazı / tura atışının sonucuna) ve yazı / tura atışları için olasılık dağılımı belirleyen bir dağılıma (örn. Bernoulli dağılımı) bakarak, belirlediğimiz anlamlılık düzeyinde (örn. 0.05) bu hipotezi reddedip reddedemeyeceğimizle ilgileniriz. Zarın olası tüm hilelilik durumlarının (hep yazı gelir, hep tura gelir, 3/4 oranında tura gelir vs.) olasılıklarını belirlemekle ilgilenmeyiz.

Şimdi aynı probleme Bayesian istatistik yaklaşımıyla yaklaşalım ve bu yolla Bayesian istatistiğin temellerini anlamaya çalışalım.

Problem: Bir kumarhanede yapılan 4096 seri yazı-tura atışının çoğunun (kaç tanesinin?) yazı gelmesiyle yapılan yazı / tura atışlarının hileli olmasından şüpheleniyoruz. Taraflı olabileceğinden (ing. bias-weighted) kuşkulandığımız bu atışlar konusundaki kuşkularımızı gidermek üzere, tüm yazı atışlarının $T = 0$ ve tüm tura atışlarının $T = 1$ ile temsil edildiği bir skala oluşturuyoruz. $T = 0.5$ bu durumda adil bir atışı temsil etmiş oluyor. Probleme Bayes teoremi ile yaklaşarak elimizdeki veri (4096 yazı / tura atışının sonucu), hiç veri görmeden her bir hileli durumun olasılığına dair inancımızı gösteren öncül (ing. prior) ve arkaplan (ing. background) bilgimiz (bir atışının sonucunun diğerini belirlemediği: $I$) ışığında atışların adil olup olmadığını belirleyen $T$ fonksiyonu için seçtiğimiz herhangi bir değerin ne kadar olası olduğudur (ardıl: posterior): $P(T | veri, I)$.

Öncül Fonksiyon (Prior): Elimizde hiçbir veri olmaksızın $T$’nin her bir değeri için sadece arkaplan bilgimize dayanarak atadığımız olasılık değeri bizim $T$’nin her bir değerinin olasılığına dair inancımız ya da önyargımızdır. Başlangıç olarak tüm hileli durumların eşit olasılıkta olduğunu varsayalım. Sonuçta bir kumarhanede neyin nasıl olacağını bilemeyiz. O nedenle $T$’nin her bir değerine aynı olasılığı veren tekdüze (ing. uniform) bir olasılık dağılım fonksiyonu kullanalım. Bu fonksiyon aşağıdaki şekilde ifade edilebilir.

$$ P(T | I) = 1, 0 \leq T \leq 1 $$, $$ P(T | I) = 0, T < 0 , T > 1 $$

Yazi Tura Deneyi - I

Yukarıda belirli bir sayıda deney yapılmadan önce paranın hangi olasılıkla tura (T) geleceği konusundaki inancımızı gösteren bir tekdüze dağılım görülmektedir (ing. prior distribution). $0$ tüm atışların yazı, $1$ tüm atışların tura gelme olasılığını göstermektedir. $n$ deney sonucunda $0$ ile $n$ arasında tura gelme olasılığı eşittir. Bu tür, herhangi bir teoriye ya da öngörüye dayanmayan öncüllere bilgi vermeyen (ing. uninformative) öncüller adı verilir. Sağ üst köşedeki 0, henüz atış yapmadığımızı göstermektedir.

Olabilirlik Fonksiyonu (Likelihood Function): Verilen bir hileli durum için (T’nin herhangi bir değeri için) elimizdeki veri setini (4096 yazı / tura atışının sonuçlarını) oluşturabilme olasılığımızı belirleyen fonksiyon olabilirlik fonksiyonudur. Arkaplan bilgimiz ışığında bunu hesaplayabiliriz. Arkaplan bilgimiz (I) bir atışın sonucunun bir sonrakini belirlemediğidir. Tıpkı klasik istatistikte olduğu gibi her bir hileli durum için elimizdeki veri setini üretebileceğimiz bir dağılım öngörüyoruz. Yazı / tura atışı gibi iki seçenekli durumlarda olasılıkların dağılımını Bernoulli dağılımı ile belirlediğimizi görmüştük ($P(x) = p^m (q)^{n-m}, q = 1 – p)$). Problemimizi kolaylıkla buna adapte edebiliriz. $p$ önermesini yazı / tura atışına atadığımız hilelilik durumu $T$, atışların tura gelmeye ayarlanma oranını gösterdiği için; $1-T$ de diğer seçenek yani yazı gelmeye ayarlanma oranını (yani $q$) göstermektedir ve bu ikisinin toplamı 1’dir. Atış ne şekilde hileli olursa olsun para ya yazı ya da tura gelir, üçüncü bir seçenek yoktur. Paranın $R$ denemede tura geldiğini ve $N$ kez atıldığını varsayarak $N-R$ denemede de yazı geldiği sonucu çıkar. Bu durumda olabilirlik fonksiyonu

$$ P(veri~|~T, I) = T^R~(1 – T)^{N-R} $$

olarak belirlenmiş olur.

Kanıt (Evidence): Her bir $T$ durumu için elimizdeki veri setini (4096 yazı / tura atışının sonuçlarını) üretebilme ihtimalidir. Gördüğümüz gibi kanıt olabilirlik fonksiyonuna çok benzemekte, ondan tüm durumlar için olasılıkların toplamını belirlemesi bağlamında ayrılmaktadır. Yani olabilirlik fonksiyonu herhangi bir hileli durum ($T$’nin herhangi bir değeri) için bu veri setinin türetilme ihtimalini belirlerken, kanıt (evidence) tüm hileli durumlar için bu veri setinin türetilme olasılıklarının toplamıdır. Yani aslında T’nin herhangi bir değeri için olabilirlik fonksiyonunu da içeren ve tüm olasılıkları kapsayan bir normalizasyon terimidir. Sonuç olarak belilremek istediğimiz elimizdeki veriyle her bir $T$ durumunun (tek tek) gerçekleşme olasılğını gösteren ardıl olasılık dağılımıdır (ing. posterior probability distribution). Bu nedenle seçilmiş bir $T$ durumu için öncülümüzün (prior) ne öngordüğü ile olabilirlik fonksiyonunun bu $T$ durumu için elimizdeki veriyi üretme konusunda belirlediği olasılığın çarpımını tüm olası $T$ değerleri için bu veriyi üretme olasılığının toplamını ifade eden kanıta bölerek, söz konusu $T$ durumunun elimizdeki veri ile üretlimiş olma olasılığını (ardıl olasılık) bulmuş oluruz. Bunu tüm $T$ durumları için yapar ve toplarsak, toplamda 1 elde etmeliyiz! Görüldüğü gibi her bir ardıl olasılığın hesabı için hep aynı kanıtı (evidence) kullanıyoruz. O nedenle her bir $T$ durumunun doğrudan olasılığı yerine bunların birbirlerine göre ne kadar olası olduklarını karşılaştırmak istersek kanıtı hesaplamaya ihtiyacımız olmaz. Bayes teoremini bir eşitlik olarak değil bir bağıntı olarak yazmamız yeter:

$$ P(T~|~veri, I)~\alpha~P(veri~|~T, I) \times P(T~|~I) $$

İlk yazı / tura atışı: Örnek olarak ilk yazı / tura atışının ardıl (posterior) dağılımını hesaplayalım ve ilk atışın tura geldiğini varsayalım. Bu durumda, $P(veri | T, I) = T^R~(1 – T)^{N-R}$ şeklinde verilen olabilirlik fonksiyonunda $R = 1$ ve $N = 1$’dir. Olabilirlik fonksiyonu bu durumda,

$$ P(veri~|~T, I) = T^1~(1 – T)^{1-1} = T $$

olarak belirlenmiş olur. Öncülümüz T’nin her değeri için eşit ve 1 idi.

$$ P(T~|~I) = 1 $$

Ardıl olasılık dağılımı (posterior probability function) bu durumda aşağıdaki gibi bulunur.

$$ P(T~|~veri, I)~\alpha~P(veri~|~T, I) \times P(T~|~I) = T \times 1 = T $$

Yazi Tura Deneyi - II

Şeklilde y-ekseni olasılık olarak belirtilmiş olmasına karşılık herhangi bir ölçek verilmemiştir. Bunun nedeni eksenin gerçekte bir Olasılık Yoğunluk Fonksiyonu (ing. Probability Density Function, PDF) biriminde olması ve olasılıkların eğrinin altında kalan alanlarla belirlenmesi ve kanıt hesaplanmadığı için bunların göreli olasılıklar olmasıdır. Eğer her bir $T$ durumu için veriyi elde etme olasılıkları (kanıt) hesaplanır ve hesaba dahil edilirse eğrinin altında kalan toplam alan 1’e eşit çıkar. Sol üst köşede en son yapılan atış ($T$: tura), sağ üst köşede ise toplam deney sayısı (yazı / tura atışı sayısı) verilmektedir. Görüldüğü gibi atışların tümünün yazıya ayarlanmış olma ihtimali kalmamıştır ($P (T = 0 | veri, I) = 0$). En büyük olasılık ise tüm atışların turaya ayarlanmış olması durumuna verilmiştir, zira şu ana kadar yapılmış olan tek atışın sonucu turadır.

İkinci yazı / tura atışı: Şimdi ikinci yazı / tura atışın sonrası ardıl (posterior) dağılımını hesaplayalım ve yine tura geldiğini varsayalım. Bu durumda $R = 2$ ve $N = 2$’dir. Olabilrlik fonksiyonu,

$$ P(veri~|~T, I) = T^2 (1 – T)^{2-2} = T^2 $$

olarak belirlenmiş olur. Öncülümüz $T$’nın her değeri için eşit ve 1 idi.

$$ P(T~|~I) = 1 $$

Ardıl olasılık dağılımı (posterior probability function) bu durumda aşağıdaki gibi bulunur.

$$ P(T~|~veri, I)~\alpha~P(veri~|~T, I) \times P(T~|~I) = T^2 \times 1 = T^2 $$

Yazi Tura Deneyi - II

Görüldüğü gibi $T > 0.5$ durumları (atışların çoğunun tura gelmeye ayarlanmış olma durumu) giderek daha olası hale gellirken $T < 0.5$ durumlarının olasılığı giderek azalmaktadır. Zira ilk 2 atışın 2’si de tura gelmiştir. Yine en büyük olasılık tüm atışların turaya ayarlanmış olması durumundadır ($T = 1$).

Üçüncü yazı / tura atışı: Şimdi üçüncü yazı / tura atışı sonrası ardıl (posterior) dağılımını hesaplayalım ve bu kez (sonunda!) yazı geldiğini varsayalım. Bu durumda $R = 2$ ve N = 3’tür. Olabilrlik fonksiyonu;

$$ P(veri~|~T, I) = T^2~ (1 – T)^{3-2} = T^2~(1 - T) $$

olarak belirlenmiş olur. Öncülümüz T’nin her değeri için eşit ve 1 idi.

$$ P(T~|~I) = 1 $$

Ardıl olasılık dağılımı bu durumda aşağıdaki gibi bulunur.

$$ P(T~|~veri, I)~\alpha~P(veri~|~T, I) \times P(T | I) = T^2~(1 - T) \times 1 = T^2~(1 - T) $$

Yazi Tura Deneyi - III

Artık tüm atışların tura gelmeye ayarlanmadığını da biliyoruz. Zira son atışımız yazı (Y) gelmiştir: $P(T = 1~|~veri , I) = 0$.

Ancak hala atışların çoğunun tura gelmeye ayarlandığı düşüncesi hakimdir ve ardıl dağılımda açıkça görülmektedir. $P(T > 0.5~|~veri, I) > P(T < 0.5~|~veri, I)$.

64. yazı / tura atışı sonrası: Deneyi yapmayı sürdüyor, her durumda ardıl dağılımı hesaplıyor ve grafiğe geçiriyoruz. Burada y-ekesninde olasılık yoğunluğu dağılımının göreli bir ifadesinin olduğunu hatırlatalım. Aksi durumda hesaba kanıtı da dahil etmemiz gerekirdi. Diyelim ki 64. atış sonunda aşağıdaki gibi bir ardıl olasılık dağılım fonksiyonu elde etmiş olalım.

Yazi Tura Deneyi - IV

Gördüklerimiz giderek bizi yazı / tura atışlarının adil olmadığına "ikna" ediyor. Başlangıçta tüm olasılıklara eşit şans tanırken 64 atış sonunda tüm atışların neredeyse 3/4’ünün yazı gelmeye ayarlanmış olduğu kanısına varıyoruz. Hala diğer durumlar da olası ama olasılıkları daha az!

4096. yazı / tura atışı sonrası: Tüm atışlar tamamlandığında artık atışların hileli olduğundan neredeyse eminiz ve bu hilelilik oranı da $T = 0.25$ (ki bu dağılımın standart sapması üzerinden belirsizliğini de belirleyebilirdik). Yani atış dört atışın birinde Tura, üçünde Yazı gelmeye ayarlanmış.

Yazi Tura Deneyi - IV

Başa Dön

Farklı Öncül Fonksiyonların Seçimi

Acaba tekdüze (ing. uniform, flat) bir öncül (prior) yerine farklı öncüller kullansak sonuç değişir miydi? Örneğin parayı adil varsaysak ve maksimum olasılığı $T = 0.5$ olan ve $T = 0.35$ ile $T = 0.65$ arasında olasılığın daha fazla olduğu bir dağılımla yola çıksak (aşağıdaki şekilde --- kesikli eğri), ya da parayı inanılmaz yanlı olarak hayal ettiğimiz ve maksimum ağırlığı paranın tamamen yazı ya da tamamen tura gelme olasılıklarına doğru dağıttığımız bir olasılık dağılımı (aşağıdaki şekilde .... noktalı eğri) ile yola çıksak ne olur? Karşılaştırma için her durumda elde edilen ardıl olasılık dağılım fonksiyonları (posterior pdf) aşağıdaki şekilde verilmektedir. Deney sayısı arttıkça aynı sonuçları elde etmemiz durumunda (yani veri setimizin değişmemesi durumunda) dağılımın aynı ardıla (posterior pdf) doğru evrildiği açıktır.

Yazi Tura Deneyi - IV

Şekle bakılarak her dağılımın aynı maksimum olasılık değerlerini verdiği düşünülmemelidir. Daha iyi bir karşılaştırma için ardıl (posterior pdf) değerleri tekdüze dağılımla elde edilene ölçeklenmiş durumdadır. Ancak sonuç olarak öncül bilgimiz ve yaklaşımımız ne olursa olsun veri arttıkça olabilirlik fonksiyonunun öncülü domine etmesi nedeniyle olayın doğasının gerektirdiği olasılık dağılımına yakınsanmaktadır. Başlangıçtaki öncüle bağımlılık giderek azalmakta başlangıç inancımız ne olursa olsun deneysel kanıtın artmasıyla aynı sonuçlara ulaşılmaktadır.

Şekilde görülen bir başka sonuç ancak bin küsür atış sonrası paranın doğasına yönelik olan inancımızın güçlendiğidir. Diğer bir husus da paranın tüm hileli durumlara eşit olasılıkla sahip olduğu durum (sürekli eğri) ile paranın neredeyse tamamen hileli varsayıldığı (noktalı eğri) durum hızla birbirine yaklaşır ve bir $T$ değerini yakınsarken, hemen hemen adil olduğunun varsayıldığı öncülle yola çıkıldığında (kesikli eğri) yakınsamanın ancak 2000 deneyden sonra gerçekleşmesidir. Birinci gözlemimizin nedeni paranını hilelilik derecesinin çok yüksek olmamasıdır, yoksa hemen farkedilirdi. İkinci gözlememiz ise kesikli eğriyle gösterilen öncülün tekdüze ve hileliği yüksek para yaklaşımlarına göre daha az belirsizlik içermesiyle açıklanabliir. Zira tekdüze öncül zaten tamamıyla düzdür. Hileliği yüksek zar yaklaşımı (noktalı eğri) da $T$'nin çok fazla değeri için neredeyse düzdür. Kesikli eğri ise Gaussyen'e benzeyen bu nedenle de adil bir zar yaklaşımına daha yakın bir yaklaşımdır. O nedenle de zarın hileli olduğuna çok daha zor "ikna" olmaktadır.

Öncüller zarın hileli olduğuna ilişkin başlangıç varsayımları (önyargıları) farklı 3 insan gibidir. Başlangıçta bütün veriyi kendi önyargılarının gözünden değerlendirmektedirler. Ancak bu 3 insan da ne kadar önyargılı başlasalar da "makul" insanlardır ve veriyi gördükçe yaklaşımlarını değiştirmektedirler ve ikna olmaya açıktırlar. En zor ikna olan doğal olarak başlangıçta zarların hileli olmadığı konusunda en önyargılı olan olacaktır. Dolayısı ile öncül bütünüyle önemsiz değildir. Veri sayısına ve bilgi düzeyimize göre yanlılıklarımız çıkarımlarımızda önemli rol oynayabilir. Doğal olan da budur ve gerçekten biz bu deneyimi yaşarız. Gerçekten en uzakta yanlılığa sahip olanlar inatçılıkları oranında gerçeğe en az uyananlardır. Ama gerçek orada durur ve daha çok veri gördükleri vakit o gerçeğe eninde sonunda "uyanırlar".

Başa Dön

Adım Adım Ya Da Tek Adımda Veri Analizi

Acaba $N$ tane örnekten oluşan $D_k = (D_1, D_2, D_3, .., D_N)$ veri setini (N atıştan oluşan tek bir yazı / tura deneyi) tek bir kere de alıp analiz etmekle $(P(T | D_k, I))$, biraz önce yazı / tura deneyinde yaptığımız gibi veri geldikçe analiz etmek ($P(T | D_1, I)$, $P(T | D_2, I)$, ... ) arasında bir fark var mıdır? Öncelikle durumu sadece 2 veri açısından ($D_1$ ve $D_2$) değerlendirelim.

İlk olarak bu iki veriyi aynı anda analiz ettiğimizi varsayalım. Bu durumda ardıl olasılık yoğunluk fonksiyonu:

$$ P(T | D_2, D_1~|~I)~\alpha~P(D_2, D_1~|~T,I) \times P(T~|~I) $$

Bu kez diyelim ki öncelikle $D_1$'i gözledik ve $D_1$ ile öncülümüze dayalı olarak $D_2$ verisi de alındıktan sonra ardıl olasılık yoğunluk fonksiyonunu hesaplıyoruz. Bu kez,

$$ P(T | D_2, D_1~|~I)~\alpha~P(D_2~|~T,D_1,I) \times P(T~|~D_1,I) $$

Bu ifadenin en sağındaki öncülün $D_1$ ile yapılan analizin sonucunda elde edilen ardıl olasılık yoğunluk fonksiyonu (posterior) olduğu kolayca görülebilir. Bu ifadedeki olabilirlik fonksiyonu ise klasik kullandıklarımızdan bir miktar farklı olmakla birlikte verilerin birbirinden bağımsız oldukları, yani bir veriyle ulaşacağımız sonuç ikincisiyle ulaşacağımızdan (ya da tersi) bağımsızdır. Örneğin beşinci kez para atışı soncunda yapacağımız analiz dördüncünün sonucundan bağımsız olmalıdır. Bu durum matematiksel olarak şu şekilde ifade edilebilir:

$$ P(D_2~|~T,D_1,I) = P(D_2~|~T,I) $$

Bu ifade yukarıdaki ikinci ifadede yerine konduğunda 1. ifade elde edilir ki bu adım adım ilerleyerek analizle, tek bir seferde tüm verinin analizi arasında bir fark olmadığını ortaya koyar.

Başa Dön

İteratif Bir Yöntemle Daha İyi Sonuç Elde Edilebilir Mi?

Burada iteratif yöntemden kasıt aynı verinin analizinde bir adımda elde edilen ardılın (posterior) bir sonraki adımda öncül (prior) olarak kullanılmak suretiyle izlenen bir yöntemdir. Böyle bir yöntem yanlış ve yanıltıcı sonuçlar verir! Bu tür bir "boot-strapping" yöntemi bir ardıl olasılık yoğunluk fonksiyonunu tekrar kendisine bağlar. Veri değişmediği için bu iki ardıl fonksiyon birbirine ancak eşit olabilir, bir olabilirlik fonksiyonu ile bağlanamaz. Eğer ısrar edilirse iterasyonun birinci adımında ulaşılan ardılın öngördüğü en muhtemel değer etrafında giderek keskinleşen dağılım fonksiyonları elde edilir ve veri analizcisi yaptığının doğru olduğuna yanlış bir şekilde giderek daha fazla inanır. Bir analizi iyileştirmenin tek yolu daha fazla veri almak ya da veri az ise gerçekçi öncül ve olabilirlik fonksiyonları kullanmaktır.

Parametreler İçin En İyi Tahmini Değerler, Hataları ve Güven Aralıkları

Ardıl olasılık fonksiyonunun (ing. posterior pdf) analiz sonuçlarını nasıl ortaya koyduğunu gördük: parametrenin her olası değer aralığı için bir olasılık değeri belirleyerek! Ancak bazen bir parametre için en olası değere ve bu değer için de bir güven düzeyine (ing. confidence level) ihtiyacımız olur. En iyi tahmin açık ki ardıl olasılık dağılımının aldığı maksimum değere denk gelen parametre değeridir. Matematiksel olarak bu değeri ($X_0$) ardılın ($P(X~|~veri, I)$) türevini alıp $0$'a eşitleyerek bulabiliriz.

$$ \frac{dP}{dX}|_{x_0} = 0 $$

Yine matematiksel olarak bu noktanın maksimum olduğunu güvence altına almak için ikinci türevin $0$'dan küçük olup olmadığı da kontrol edilmelidir. Süreksiz durumda iş daha kolaydır zira y-ekseni doğrudan olasılıktır (posterior probability mass function). En olası değer de o maksimum olasılığın olduğu $X$ değeri olur.

Bu değerin güven düzeyi için ise ardılın $x_0$ civarında nasıl dağıldığına bakılır. Herhangi bir fonksiyonun bir nokta etrafında nasıl dağıldığına bakmanın iyi bir yolu, fonksiyonun o nokta civarında Taylor açılımını elde etmektir. Aslında bu o nokta civarında fonksiyona düşük dereceden bir polinom fiti yapmakla özdeştir. Genellikle bu durumda pdf'in kendisinin yerine daha yavaş değiştiği için onun logaritması ile ilgilenilir ($L = ln(P(X~|~veri,I)$). Bu durumda ardılın $x_0$ civarındaki Taylor açılımı:

$$ L = L(x_0) + \frac{1}{2} \frac{d^2L}{dX^2}|_{x_0} (X - x_0)^2 + ... $$

Bu tanım dahilinde en iyi tahmin $dL / dX = 0$'ın karşılık geldiği $x_0$ değeri olur. Maksimum civarında bir seri açılımı olduğu için lineer terim bu seride bulunmaz. Birinci terim ($L(x_0)$) da bir sabit olup ardılın şekliyle ilişkili değildir. Dolayısı ile ardılın şeklini (en olası değer etrafındaki saçılmasını) kuadratik terim temsil eder. Diğer tüm üst dereceden terimleri çok küçük olacakları gerekçesiyle ihmal edecek olursak

$$ P(X~|~veri,I) = A~exp[\frac{d^2L}{dX^2}|_{x_0} (X - x_0)^2] $$

Tekrar ardıla geçmek için $L$'nin eksponansiyelini almalıyız. Bu durumda $x_0$ civarında $pdf$ yukarıdaki fonksiyonla ifade edilmiş olur. Aslında bu şekilde $x_0$ civarında seriyi Taylor serisinie açmakla $x_0$ civarında ardılı bir Gaussyen fonksiyonla ifade etmiş oluruz. Burada $A$ bir normalizasyon sabitidir. Ardıl böylece bir Gaussyen, bir başka deyişle bir normal dağılım olmuş olur.

Gaussyen PDF

Bu dağılımın standart sapması $L$'nin $X = x_0$'daki ikinci türevinin karekökü ile ters orantılıdır.

$$ \sigma = (-\frac{d^2L}{dX^2}|_{x_0})^{-1/2} $$

$X$ parametresi için en iyi tahmin (en olası değer; best estimate) ise klasik şekilde verilir $X = x_0 \pm \sigma$. Doğal olarak normal dağılımın gereği $X$'in gerçek değerinin $x_0 \pm 1~\sigma$ arasında olma ihtimali $\%67$'dir. Örneklemde yeterince veri olması durumunda medyan da en olası değer için kullanılabilir.

$$ P(x_0 - \sigma \leq X \leq x_0 + \sigma~|~veri,I) = \int_{x_0 - \sigma}^{x_0 + \sigma} P(X~|~veri,I) dX \approx 0.67 $$

Örnek: Yazı / Tura Probleminde En İyi Değer Hesabı

Tekdüze (ing. uniform, flat) öncül (prior) ile binomial olabilirlik fonksiyonlarının (likelihood function) kullanarak hileli madeni parayla yazı / tura atışı örneğinde ardıl olasılık yoğunluk dağılımı fonksiyonunu

$$ P(T~|~veri,I)~\alpha~T^R~(1 - T)^{N-R} $$

şeklinde elde etmiştik. Bu ifadenin doğal logaritmasını alarak $L$ fonksiyonunu hesaplayacak olursak

$$ L = sbt~+~R~ln(T) + (N-R)~ln(1 - T) $$

ifadesi elde edilir. Bu ifade değişkeni olan $T$'ye göre iki kez türevlenirse,

$$ \frac{dL}{dT} = \frac{R}{T} - \frac{N-R}{1-T} \rightarrow \frac{d^2L}{dT^2} = -\frac{R}{T^2} - \frac{N-R}{(1-T)^2} $$

elde edilmiş olur. En olası değeri bulmak için birinci türev $0$'a eşitlenir ve cebirsel düzenleme yapılacak olursa,

$$ T_0 = \frac{R}{N} $$

bulunmuş olur. Bu basit bir şekilde en olası hilelilik değerinin tura sayısının toplam yazı/tura atış sayısına oranı olduğunu gösterir. Dağılımın standart sapması ise $L$'nin ikinci türeviyle verilir.

$$ \frac{d^2 L}{dT^2}|_{T_0} = -\frac{N}{T~(T - 1)}$$

Standart sapma hatırlanacağı gibi bu değerin tersinin karaköküdür.

$$ \sigma = \sqrt{\frac{T~(T - 1)}{N}} $$

$T$ (atışın turaya yanlılık oranı) bir süre deney yaptıktan sonra $T_0$'a doğru yakınsar ve ardılın maksimumu $T_0$ civarına gelir. Ancak bu maksimum etrafındaki dağılımın standart sapmasının azalması (yani bu maksimumun daha güvenilir hale gelmesi) daha çok veri almakla ($N$'i büyütmekle) mümkündür. Ayrıca bu denklemden adil bir atışın ($T \sim 0.5$) yanlılık oranını bulmak çok daha zordur; zira $T_0 = 0.5$ için standart sapma ifadesinin bölüm çizgisinin üstünde kalan tarafı ($T_0 \times (1 - T_0)$) maksimum olur. Dağılımın genişliği de bu değer için doğal olarak maksimumdur!

Astronomide Modelleme ve Bayesian İstatistik

Astronomide (aslında tüm bilimlerde) istatistik yaklaşımlara başvurmamız gerektiğinde (örneğin geçerli bir teoremiz olmadığında) en temel problem verinin yetersizliğidir. Örneğin galaksimizdeki ötegezegenlerin oluşumuna ilişkin bir çalışma yapmak istediğimizde tüm ötegezegen sistemlerine ilişkin tüm parametrelere (popülasyon) sahip olmadığımız için bugüne kadar keşfedilenlerin bilinen parametreleri (örneklem) üzerinden çıkarım yapmak zorunda kalırız.

Örneğin Güneş-benzeri salınımlar gözlenen yıldızlarda ortalama büyük ayrışma (ing. large separation) gibi bir parametreyle ilgileniyorsanız, ancak bu parametrenin hesaplandığı (örneğin Kepler örneklemi) örneklemle sınırlı kalırsınız. Örneklemeniz gelişip (örneğin TESS örneklemi) parametrenin yeni bir ortalama değerini elde ettiğinizde bu değer de tüm Güneş-benzeri yıldızlarınkiyle (popülasyon) aynı olmaz. Buna örneklem hatası (ing. sampling error) diyoruz. Frekansçı yaklaşımda popülasyonun ortalama büyük ayrışması sabit bir değerken, değişen şey veri setidir. Bayesian istatistiksel yaklaşımda ise popülasyon gibi bir bilgi olmadığı için, elde ettiğimiz şey elimizdeki veriye dayanarak ortalama büyük ayrışmanın alabileceği tüm değerlerin olasılıkları, yani bir olasılık dağılımıdır. Dolayısı ile değişenimiz aslında ortalama büyük ayrışma olur, sabit olansa elimizdeki veri setidir.

Ardil PDF

Dolayısıı ile hesaplamaya çalıştığımız şey elimizdeki veriyle herhangi bir ortalama ayrışma değerinin elde edilme olasılığıdır ($P(\Delta \nu~|~veri$).

Bunun için öncelikle $\Delta \nu$'nün herhangi bir değeri için elimizdeki verinin türemiş olma ihtimaline bakarız (olabilirlik fonksiyonu, likelihood function). Daha sonra eğer hiç veri almamış olsaydık ortalama ayrışma için her bir değerin olasılığının ne olması gerektiğini düşünürüz (öncül, prior). Ortalama ayrışma için her bir değerin olasılığının ne olması gerektiği teoriden geliyor olabilir (ing. informative), tahmini bir değer etrafında rastgele dağıldığı öngörülebilir (ing. random) veya her değeri alma ihtimali aynı olabilir (ing. uninformative). Son olarak da her bir ortalama ayrışma için bu ayrışma değerini elde etme olasılığımıza bakarız (kanıt, evidence). Sonuçta elde etmek istediğimiz de herhangi bir ortalama ayrışma değerinin elde edilme olasılığı (ardıl, posterior) ve hangi ortalama ayrışma değerinin en olası olduğudur (ardılın maksimum değeri).

$$ P(\Delta \nu~|~veri) = \frac{P(veri~|~\Delta \nu) \times P(\Delta \nu)}{P(veri)} $$

Gerçekte bir modelimiz olmadan ne herhangi bir ortalama ayrışma değeri için elimizdeki veriyi elde etme olasılığını (olabilirlik fonksiyonunu), ne her bir ortalama ayrışma değerinin olasılığını (öncül), ne de tüm ortalama ayrışma değerleri için olasılık hesaplamamız (kanıt) mümkündür. Bu nedenle bu ifadeye bir modelin varlığında bu olasılıkların hesaplanacağını girmeliyiz.

$$ P(\Delta \nu~|~veri,model) = \frac{P(veri~|~\Delta \nu,model) \times P(\Delta \nu~|~model)}{P(veri)} $$

MCMC ve Bayesian Analiz Örnekleri

Örnek 1. Basit Veri Örneği

Önce basit bir örnekle Bayesian analiz yaklaşınının parametrelerini öğrenmeye çalışalım. Bunun için öncelikle $\mu = 10$, $\sigma = 3$ ile tanımlı bir normal dağılımdan $N = 30000$ örnek sayılı bir örneklem oluşturalım ve bu örneklemden de sadece rastgele seçilmiş $n = 1000$ örneğin gözlenebildiğini varsayalım. Bu örneğin amacı, gözlenen bu 1000 örnekten yola çıkarak dağılımın standart sapması ($\sigma$) için bir dağılım bulmak ve bunu Bayesian istatistik paradigmasından yararlanarak yapmaktır.

Doğal olarak Gaussyen bir dağılımın standart sapmasının bir ifadesi olduğu ve bu ifadeyle kolaylıkla hesaplanacağı düşünülebilir. Ancak amacımız Bayesian yaklaşımın temellerini öğrenmek üzere seçilmiş bu basit örnekte standart sapmayı ($\sigma$) doğrudan hesaplamak yerine onun için bir olasılık dağılımı elde etmek! Çünkü, bu örnekte olduğu gibi popülasyon bilgisine her zaman sahip olmadığımız için basit ifadelerle standart sapmasını hesaplamak da mümkün olamayacaktır. Bu durumda yapabileceğiniz "en iyi şey" eldeki örneklemden (veriden) hareketle popülasyonun standart sapması için bir dağılım elde etmek olabilir.

In [1]:
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
mu = 10.0
sigma = 3.0
N = 30000
n = 1000
np.random.seed(758)
populasyon = np.random.normal(loc=mu,scale=sigma,size=N)
# bu populasyondan rastgele 1000 ornek secelim
orneklem = populasyon[np.random.randint(0, N, n)]
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1,1,1)
ax.hist(orneklem,bins=35)
ax.set_xlabel("Deger")
ax.set_ylabel("Frekans")
ax.set_title("$\mu = 10, \sigma = 3$ bir normal dagilimdan (N = 30000) rastgele secilmis 1000 gozlemin histogrami")
xbar = orneklem.mean()
str2wr = r'$\overline{x}$' + " = {:.2f}".format(xbar)
plt.text(17.5,80.0,str2wr)
plt.show()

Kendi yarattığımız gözlemsel veri (kanıt, ing. evidence) ve hiçbir veri görmeden olasılık dağılımını aradığımız parametreye ilişkin inancımızın derecesini ifade eden öncül olasılık dağılımını (prior) kullanarak ilgilendiğimiz parametrenin ($\sigma$) ardıl olasılık dağılımını (posterior) bulmak istiyoruz. Bu amaçla parametre uzayından rastgele örnekler seçip, ilgilendiğimiz parametrenin dağıılımını incelememizi sağlayan Monte-Carlo Markov Chain (MCMC) yöntemini daha önce gördüğümüz Metropolis-Hastings algoritmasına dayalı olarak kullanacağız.

Bu yöntemde parametre uzayından (standart sapmalardan) örnekler almamıza yardımcı olacak bir öneri dağılımına (ya da geçiş modeline) (ing. proposal distribution, transition model) ($Q(\theta^\prime / \theta)$) ihtiyaç duyulur. Bu dağılım ya da model her adımda parametre (örneğimizde standart sapma) için yeni bir değer önererek bir parametre değerinden ($\theta$) diğerine ($\theta^\prime$) geçilmesini sağlar.

Metropolis-Hastings algoritması bu $Q$ fonksiyonunu kullanarak parametre uzayında rastgele yürüyüşünü (ing. random walk) yapar ve her bir adımı atıp atmayacağına rastgeleliği göz ardı etmeden ancak bir sonraki adımın daha olası olup olmadığına bakarak karar verir. Bu olasılık (her bir adımın ne kadar olası olduğu) olabilirlik fonksiyonu (ing. likelihood function) ile ($f = P( veri = D | \theta)$) belirlenir.

İlgilenilen parametre için herhangi bir adımdaki parametrenin değeri ($\theta$) alınır ve $Q(\theta^\prime / \theta)$ öneri dağılımı kulanılarak bir sonraki adım için parametrenin yeni değeri belirlenir ($\theta^\prime$). Bu model genellikle $Q(\theta^\prime / \theta) = N(\bar{x},s)$) şeklinde ortalaması $\bar{x}$, standart sapması ise $s$ olan bir normal dağılımdır. Her bir adımın atılıp atılmamasına bu iki parametrenin ne kadar olası olduklarının (ardıl olasılıkların) oranlarına bakılarak karar verilir.

$$ \frac{P(\theta^\prime~|~D)}{P(\theta~|~D)} $$

Bayes formülü bu oranı hesaplamak için olabilirlik fonksiyonu ($P( veri = D~|~\theta)$) ve öncül ($P(\theta)$) temelinde kolaylıkla kullanılabilir. Elde edilecek olan ardıl dağılımların bir oranı olacağı ve her iki ardıl da aynı veriye dayanarak belirleneceği için tüm olası modeller (ya da parametre değerleri) için verinin türetilme olasılığını veren kanıta (ing. evidence) ihtiyaç duyulmaz. Zira oran alındığında kanıt terimleri sadeleşir.

$$ \frac{P(D~|~\theta^\prime) \times P(\theta^\prime)}{P(D~|~\theta) \times P(\theta)} $$

Sonuç olarak Metropolis-Hastings algoritması Bayes formülünü kullanarak olabilirlik fonksiyonu ($P( veri = D~|~\theta)$) ve öncül ($P(\theta)$) temelinde hesaplanan bu oranı parametre uzayında adım atıp atmamak üzere karar vermek amacıyla kullanır ve kabul edilen parametrelerin bir dağılımının elde edilerek ardıl dağılımın bir örneklemine ulaşılmasını sağlar.

$$ \frac{P(\theta^\prime~|~D)}{P(\theta~|~D)} = \frac{P(D~|~\theta^\prime) \times P(\theta^\prime)}{P(D~|~\theta) \times P(\theta)} $$

Bu oran $1$’den büyükse $\theta^\prime$ modelinin ya da parametresinin daha olası olduğu değerlendirmesi yapılarak adım kabul edilir. Değilse, bu durumda $(0,1)$ aralığında rastgele bir sayı üretilir. Eğer bu oran, üretilen sayıdan büyükse bu adım kabul edilir, değilse bu adım reddedilir ve parametre uzayında yeni bir $\theta$ değeri aranır. Bu algoritma istenen tekrar sayısında (ing. number of iterations) uygulanır.

Problemin teorik çerçevesini bu şekilde tanımladıktan sonra, problem için bir veri seti oluşturalım ve çözüm yöntemini bir Python koduna aktaralım.

Geçiş Fonksiyonu Seçimi

Parametre uzayından örnekler almamıza yardımcı olacak bir öneri dağılımı (ya da geçiş modeli) (ing. proposal distribution, transition model) ($Q(\theta^\prime / \theta)$) basit bir biçimde bir Gaussyen dağılım (normal dağılım) olarak seçilebilir. Ardıl olasılık dağılımını aradğımız parametrenin ($\sigma$) herhangi bir adımdaki değerini merkez alan standart sapması $\sigma^\prime = 0.5$ (aradığmız standart sapmayla ilgisi olmayıp, öneri dağılımının standart sapmasıdır) olan öneri dağılımı (ya da geçiş fonksiyonu) bir sonraki adımı atmak için aday bir standart sapmanın belirlenmesi için kullanılır.

$$ Q(\sigma_{i} / \sigma_{i-1} = Normal(\mu = \sigma_{i-1}, \sigma^\prime = 0.5) $$

In [2]:
# Bir sigmadan digerine gecmek icin kullanilacak oneri dagilimini olustur.
# Bu, gecis modeli (transition_model) olarak da isimlendirilebilir.
# Ornegimizde bu, merkezi bir onceki sigma degerimiz, standart sapmasi ise
# 0.5 olan bir normal fonksiyon olsun.
# Burada theta model olmak uzere theta[0]: ortalama, theta[1]: sigma,
# scale: 0.5 (parametremiz olan st.sapma'nn dagiliminin standart sapmasi :)
# size: 1, parametre sayisi
gecis_fonksiyonu = lambda theta: [theta[0],np.random.normal(loc=theta[1],scale=0.5,size=(1,))]

Olabilirlik Fonksiyonu

Olabilirlik fonksiyonunu (ing. likelihood function) belirlemek amacıyla veri setindeki ($D$) her bir nokta ($d_i$) için $f$ fonksiyonu;

$$ f(d_i~|~\mu,\sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-\frac{(d_i - \mu)^2}{2 \sigma^2}} $$

şeklinde olabilirlik fonksiyonunu Gaussyen bir dağılım olarak seçmemizin önünde bir engel yoktur. Ancak işlem kolaylığı açısından da bu fonksiyonu logaritmik olarak yazmakta fayda var. Modelle ($\theta = \mu,\sigma_a$) veriyi, modelin her yeni değeri için karşılaştırıp ne kadar olası olduğunu belirleyen aşağıdaki gibi bir logaritmik normal dağılım bu ihtiyacı karşılayacaktır.

In [3]:
# Hem mevcut adim, hem de bir sonraki adimin ne kadar olasi oldugunu
# belirlemek icin kullanilacak olabilirlik fonksiyonu (logaritmasi)
# Gauus fonksiyonu fomuluyle hesaplamayi yapmaktadir.
import numpy as np
def log_olabilirlik_normal(theta,veri):
    # theta: model, theta[0]=mu, theta[1]=sigma 
    # veri = orneklem
    return np.sum(-np.log(theta[1] * np.sqrt(2* np.pi) )-((veri-theta[0])**2) / (2*theta[1]**2))

Kabul kriteri

Tüm veri setini gözönünde bulundurursak her bir noktanın olabilirlik fonksiyonu $f$ iken tüm gözlem verisi için olabilirlik fonksiyonu bu fonksiyonların tüm noktalar üzerinden bir çarpımı olmalıdır. Ayrıca bu dağılımın (Gaussyen) ortalamasını örneklemin ortalamasına eşit ($\bar{x}$) kabul ettiğimiz için aslında dağılımı aranan tek değişken standart sapması ($\sigma$) olduğundan $\theta = \sigma$’dır.

$$ \mathcal{L}(D~|~\mu,\sigma_a) = \prod_{i=1}^{n} f(d_i~|~\mu,\sigma_a) $$

Burada $a$ indisi, $\sigma$'nın ve $i-1$ ve $i.$ adımlarına karşılık gelir. $i-1.$ adımdan $i.$ adıma geçme koşulu ise bu durumda aşağıdaki şekilde tanımlanabilir.

Sonuç olarak, rastgele yürüyüşü yaparken $i$ adımındaki parametre (ya da model) değerinin mi yoksa bir önceki adımdakinin mi kabul edileceğine ilişkin bir kritere ihtiyacımız vardır ve bu kriter, doğal olarak, bu adımlardaki olabilirlik fonksiyonu değerlerinin oranının 1 ile karşılaştırılmasına bağlı olacaktır.

$$ \frac{\mathcal{L}(D~|~\mu,\sigma_i) \times P(\mu,\sigma_i)}{\mathcal{L}(D~|~\mu,\sigma_{i-1}) \times P(\mu,\sigma_{i-1})} > 1 $$

Eğer oran 1'den büyükse yeni adım daha olası olduğundan kabul edilir; değilse bu kez işe rastgelelik de getirmek için $(0,1)$ arasında rastgele bir sayı üretilir ve bu sayı $i-1.$ adımın olabilirlik fonksiyonu ile karşılaştırılarak bir sonraki adım için karar verilir.

In [4]:
def kabul_kriteri(L_i_1, L_i):
    # Eger adimi atmakla elde edilecek dagilim daha olasiysa kabul et
    if L_i > L_i_1:
        return True
    # Degilse, 0-1 arasinda bir sayi uret 
    else:
        rastgele = np.random.uniform(0,1)
        # Hesapladigimiz olasiliklar logaritmik oldugu icin
        # karsilatirmayi yapmak uzere eksponansiyel al ve
        # urettigin sayi yeni adimin olasiligi ile mevcut adimin
        # olasiligi arasindaki farktan dusukse onu,
        # degilse yeni adim ile mevcut adim olasiklari farkini dondur.
        # Boylece olasiligi dusuk yeni adimlar daha az atilan adimlar olsun
        return (rastgele < (np.exp(L_i-L_i_1)))

Öncül Seçimi

Aradığımız parametrenin (standart sapma, $\sigma$) tanımı gereği pozitif olması dışında bir kısıtı olmadığından öncül onu pozitif yapmayacak bir tekdüze dağılım olarak aşağıdaki şekilde seçilebilir.

In [5]:
def oncul(theta):
    # theta[0] = mu, theta[1]=sigma
    # sigmanin pozitif (gecerli) degerleri icin 1,
    # negatif dolayisiyla gecersiz degerleri icin 0, dondurur
    # 1 ardil'i belirlerken yapilan carpma'ya etki etmez
    # 0 ardil'i imkansiz hale getirir!
    if(theta[1] <=0):
        return 0
    return 1

Metropolis Hastings Algoritması

Her ne kadar rastgele örneklem oluşturmak için Monte Carlo Markov Zinciri hesaplamak üzere geliştirilmiş ve araştırmacılar tarafından sıkça kullanılan Python paketleri olan emcee ve pymc Metropolis-Hastings algoritması yerine daha efektif başka algoritmalar tercih ediyor olasalar ve Metropolis-Hastings algoritmasıyla MCMC hesaplayan modüller bulunsa da biz pedagojik olması açısından kendi algoritmamızı kodlayalım. Bu amaçla Ders-6 Metropolis-Hastings Algoritması başlığındaki anlatımı takip ederek bir Python fonksiyonu yazalım.

Fonksiyonumuz yazdığımız olabilirlik fonksiyonunu (log_olabilirlik_normal), öncül fonksiyonunu (oncul), geçiş fonksiyonunu (gecis_fonksiyonu) ve kabul kriterini belirleyen fonksiyonun (kabul_kriteri) yanı sıra modelin (bizim örneğimizde sadece standart sapma $\sigma$ ve sabit olduğu halde normal dağılım oluşturabilmek için gereken $\mu$ değerinden oluşan $\theta$) başlangıç parametrelerini, iterasyon sayısı ve veriyi ($n = 1000$ 'lik örneklemimiz) argüman olarak kabul etmeli ve kabul kriterine uyan model parametre değerlerini ($\theta = \mu, \sigma$) ve karşılaştırma için reddedilen parametre değerlerini döndürmelidir. Fonksiyon basit olarak her adımda geçiş fonksiyonundan yeni bir parametre değeri ikilisi ($\theta = \mu, \sigma$) alır; olabilirlik fonksiyonuna göndererek bunların ne kadar olası olduğunu veriyle karşılaştırarak hesaplar. Birbirini takip eden iki adımdaki olasılıkların oranlarını öncül fonksiyonu ile birlilkte kabul kriteri fonksiyonuna gönderir ve bu fonksiyon tarafından kabul edilen parametre değerlerini $kabul$, reddedilenleri $red$ listelerine ekleyip, bu listeler tamamlandıktan sonra birer numpy dizisine dönüştürür ve programa döndürür.

In [6]:
def metropolis_hastings(olabilirlik, oncul, gecismodeli, baslangic_param, iterasyon, veri, kabul):
    # olabilirlik(theta,veri): verilen theta modeli (mu,sigma) (mu sabit) 
    # ile eldeki veriyi turetme olasiligini hesaplayan fonksiyon
    # gecismodeli(x): oneri fonksiyonu. Yaptigi tek is bir sonraki
    # adimi belirlemek uzere simetrik (Gauss) bir dagilimdan
    # rastgele bir ornek secip dondurmek
    # baslangic_param: baslangic ornegi (sigma icin baslangic degeri)
    # iterasyon: iterasyon sayisi
    # veri: modellemeye calisilan veri
    # kabul(L_i,L_i_1): adimi kabul edip etmemeye karar veren kabul kriteri
    theta = baslangic_param
    kabul = []
    red = []   
    for i in range(iterasyon):
        # oneri dagilimiyla yeni dagilimi (theta_yeni) belirlemek uzere
        # mevcut dagilimi gecis model fonksiyonuna gonderdik.
        theta_yeni =  gecis_fonksiyonu(theta)
        # Mevcut ve yeni dagilimlar icin olabilirlik fonksiyonu hesabi
        # Bu fonksiyon istenirse elle yazilan formul (log?) istenirse de
        # scipy.stats.norm olabilir.
        theta_L = log_olabilirlik_normal(theta,veri)
        theta_L_yeni = log_olabilirlik_normal(theta_yeni,veri)
        # yeni dagilimi kabul edip etmeyecegimize yonelik karar
        if (kabul_kriteri(theta_L + np.log(oncul(theta)),theta_L_yeni+np.log(oncul(theta_yeni)))):            
            # kabul edildiyse kabul edilenler listesine ekle
            theta = theta_yeni
            kabul.append(theta_yeni)
        else:
            # degilse reddedilenler listesine ekle
            red.append(theta_yeni)                            
    return np.array(kabul), np.array(red)

Kodun Çalıştırılması ve Ardıl Olasılık Dağılımının Hesabı

Artık Metropolis-Hastings algoritmasını çalıştırıp, ardıl dağılımları elde edebilir ve sonuçları inceleyebiliriz. İterasyon sayısını 50000 seçtikten sonra öncelikle sadece kabul edilen ilk 50 parametre ($\sigma$) değerine bakalım.

In [7]:
# Algoritmayi calistir ve sonuclari incele
# Dagilim icin baslangic degerimiz 1000'lik ornegin ortalamasi
# standart sapma icin baslangic degeri 0.1,
# iterasyon sayisi 50000
from matplotlib import pyplot as plt
%matplotlib inline
kabul,red = metropolis_hastings(log_olabilirlik_normal,oncul,gecis_fonksiyonu,[mu,0.1], \
                                50000,orneklem,kabul_kriteri)
print('Ardil dagilimdaki kabul edilen parametre sayisi:{:d}'.\
      format(len(kabul)))
plt.title('Standart Sapma ($\sigma$) Icin MCMC Ile Ornekleme (Ilk 50 Ornek)')
plt.xlabel('Iterasyon Sayisi')
plt.ylabel('$\sigma$')
plt.plot(kabul[:50,1],'b+')
plt.plot(red[:50,1],'rx')
plt.show()
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in log
  
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:24: RuntimeWarning: divide by zero encountered in log
Ardil dagilimdaki kabul edilen parametre sayisi:8556

Görüldüğü gibi ilk 50 adımda dahi hemen $\sigma = 3$ değerine yakınsama gerçekleşmektedir. Kabul edilen (mavi +) değerler küçük bir standart sapmayla bu değer etrafında saçılırken, reddedilen değerlerin (kırmızı x) bu değer etrafındaki saçılması büyük olsa da onlar da giderek popülasyonun standart sapmasına yaklaşmaktadır. Burada bazı değerler negatif olabileceğinden logaritma fonksiyonu hata verebilir ancak bunu şimdilik gözardı ediyoruz. İstendiği takdirde bu hata mesajı uygun bir şekilde yönetilebilir.

Şimdi de sadece son %25 adımı (12500) inceleyelim

In [8]:
# Sadece son %25 adima bakalim
from matplotlib import pyplot as plt
%matplotlib inline
plt.title('Standart Sapma ($\sigma$) Icin MCMC Ile Orneklemeler')
plt.xlabel('Iterasyon Sayisi')
plt.ylabel('$\sigma$')
baslangic_indeksi = int(0.75*kabul.shape[0])
son_indeks = kabul.shape[0]
xekseni = range(baslangic_indeksi,son_indeks)
plt.plot(xekseni,kabul[baslangic_indeksi:son_indeks,1], 'b+')
plt.plot(xekseni,red[baslangic_indeksi:son_indeks,1], 'rx')
plt.show()

Son %25 adımda artık ısınma periyodu geçilmiş, MCMC denge dağılımına oturmuştur. Kabul edilen parametre değerleri popülasyonun standart sapması etrafında küçük bir standart sapma ile saçılırken, reddedilenler daha büyük bir standart sapma ile de olsa yine popülasyonun standart sapması etrafında saçılmaktadır. İstenirse tüm dağılım bir grafikle gösterilebilir ancak bu çok fazla noktayı grafik etmemiz gerekeceğinden size bırakılmıştır.

Son olarak ardıl dağılımın ilk %25'lik bölümünü ısınma periyodu (ing.burn-in) kabul ederek dağılımın son %75'lik bölümünden hesaplayacağımız $\theta = \mu, \sigma$ modelini ($\mu$ sabit) kullanarak popülasyonla bir karşılaştırma yapalım.

In [9]:
# Son olarak basarimi gostermek icin ardil dagiimin (son %75inin)
# parametreleriyle (mu,sigma) olusturulacak 30000 degeri iceren populasyonu
# en basta olusturdugumuz populasyonla karsilastirmasi yapildiginda
from matplotlib import pyplot as plt
%matplotlib inline
mu_075post = kabul[baslangic_indeksi:son_indeks,0].mean()
sigma_075post = kabul[baslangic_indeksi:son_indeks,1].mean()
print("Ardilin son %75'inden elde edilen mu={:.4f}, sigma={:.4f}".\
      format(mu_075post, sigma_075post[0]))
# Bu parametrelerle uretilen Gaussyen model
model_075post = lambda t, mu_075post, sigma_075post:\
                np.random.normal(mu_075post,sigma_075post,t)
# Bu modelden bastaki populasyonla ayni sayida (30000) icerecek
# sekilde turetilmis gozlem verisi
orneklem_075post = model_075post(populasyon.shape[0],mu_075post,sigma_075post)
fig2 = plt.figure(figsize=(10,10))
ax2 = fig2.add_subplot(1,1,1)
ax2.hist(orneklem_075post, bins=70, label='Ardil dagilim parametreleriyle olusturulan 30000 gozlem verisi')
ax2.hist(populasyon, bins=70, alpha=0.5, label='Orjinal populasyon')
ax2.set_xlabel("$\mu$")
ax2.set_ylabel("Frekans")
ax2.set_title("Ardil Dagilimla Orjinal Populasyonun Karsilastirmasi")
ax2.legend()
plt.show()
Ardilin son %75'inden elde edilen mu=10.0000, sigma=3.0743

Örnek Problem: Güneş Lekeleri Modeline Bayesian Yaklaşım

Bu örnek bağlamında üzerinde çalışılacak veri Ocak 1749 ile Kasım 2018 arasında kaydedilen aylık ortalama güneş lekeleri sayısıdır. Veri Sunspot Modeling Index and Long Term Solar Observations internet adresinden istenen formatta (txt, csv) kolayca çekilebilir. Bu veriyi açıp çizdirdiğimizde göreceğimiz şey aylık ortalama leke sayısının değişimi olacaktır.

In [10]:
# Oncelikle veriye bakalim
from astropy.io import ascii
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
veri = ascii.read("veri/SN_m_tot_V2.0.csv",delimiter=";",guess=False)
# Sutun 3 (MidMonthYrFrac)'teki zaman
# Ayin ortasinin karsilik geldigi yil kesri cinsinden ifade edilmistir
zaman = veri['MidMonthYrFrac']
# Sutun 4 (MonthlyMean)
# Aylik ortalama leke sayisi
# 0.1'i veride yanlislikla negatif girilen degerleri duzeltmek icin ekliyoruz
leke_sayisi = veri['MonthlyMean'] + 0.1
plt.plot(zaman,leke_sayisi,label = 'Sekil 1. Gunes Leke Sayilari')
plt.title("Ocak 1748 - Kasim 2018 Aylik Ortalama Gunes Leke Sayilari")
plt.xlabel("Zaman (Yil Kesri)")
plt.ylabel("Aylik Ortalama Leke Sayisi")
plt.legend(loc="best")
plt.show()

Leke sayılarına daha kısa bir zaman aralığını dikkate alarak daha yakından bakalım. Bu şekilde bir çevrimdeki leke sayılarını hangi fonksiyonla modelleyebileceğimize ilişkin bir fikrimiz oluşabilir.

In [11]:
# Daha yakindan kucuk bir kismina bakalim
from matplotlib import pyplot as plt
%matplotlib inline
plt.plot(zaman[2650:3200], leke_sayisi[2650:3200], label = 'Sekil 2. Gunes Lekeleri')
plt.xlabel("Zaman (Yil Kesri)")
plt.ylabel("Aylik Ortalama Leke Sayisi")
plt.title("Ocak 1970 - Haziran 2010 Aylik Ortalama Gunes Leke Sayilari")
plt.ylim((-10, 350))
plt.legend(loc="best")
plt.show()

Bu değişime yakında bakılacak ve Güneş Fiziği dersinden Güneş’in leke çevriminin ortalama periyodunun $\sim11$ yıl olduğu hatırlanacak olursa, değişimin bir $\Gamma$ dağılımı ile temsil edilebileceği görülebilir. Bu durumda değişimin $X$: aylık ortalama leke sayısını göstermek üzere $X \sim \Gamma(k ; \theta)$ ile temsil edilebilmesi için bulunması gereken dağılımın $k$ ve $\theta$ parametreleridir.

Ancak $\Gamma(k ; \theta)$ dağılımı ile her bir çevrimi ayrı ayrı temsil edebilirken genel çevrimsel değişimi temsil edemeyiz. Bunun için $\Gamma(a) = (a-1)!$ fonksiyonunun ($\Gamma$ dağılımı değil; $\Gamma$ fonksiyonu) bu çevrimsel değişime imkan verdiği bir olasılık dağılımını kullanmalıyız. Bu fonksiyon aşağıdaki şekilde tanımlanabilir.

$$ f(x; a; b) = \frac{b^a~x^{(a-1)}~e^{-bx}}{\Gamma(a)} $$

Şimdi leke sayıları histogramına bakarak bu fonksiyonun ne kadar gerçekçi bir yaklaşım olduğunu değerlendirelim.

In [12]:
# Histograma bakalim
from matplotlib import pyplot as plt
%matplotlib inline
plt.hist(leke_sayisi, bins=40, density=True)
plt.xlabel("Leke Sayilari")
plt.ylabel("Normalize Frekans")
plt.title("1749-2018 Yillari Arasindaki Gunes Leke Sayilari Histogrami")
plt.show()

Görüldüğü gibi bu dağılım, verilen $\Gamma$ fonksiyonu ile temsil edilebilir.

Problemin Çözümü ve Ardıl Olasılık Dağılımının Elde Edilmesi

Öncelikle MCMC’nin Metropolis-Hastings algoritmasıyla parametre uzayında hangi dağılımdan yeni parametreler seçeceğini belirleyen geçiş fonksiyonuna (ing. proposal distribution, transition function) ihtiyaç duyulur. Bu öneri dağılımı basit bir Guass olarak seçilebilir. Bu kez 2 ayrı parametre ($\Gamma$ fonksiyonunun $a$ ve $b$ parametreleri) için normal dağılım üretilmelidir.

In [13]:
import numpy as np
gecis_fonksiyonu = lambda theta: np.random.normal(theta,[0.05,5],(2,))

Öncül (prior) yine bilgi vermeyen (ing. uninformative) bir fonksiyon seçilebilir. Buradaki tek kısıt yine $\Gamma$ fonksiyonunun yapısı gereği $a$ ve $b$’nin negatif olamamasıdır.

In [14]:
def oncul(theta):
    if(theta[0]<=0 or theta[1] <=0):
        return 0
    else:
        return 1

Olabilirlik fonksiyonu (likelihood function) ise $\Gamma$ dağılımı kullanılmalıdır. Bu fonksiyon için yine gamma dağılımı formülü kullanılabileceği gibi, scipy.stats fonksiyonlarından gamma da kullanılabilir. Biz bu fonksiyonu daha önceden görmüş olduğumuz için onu kullanalım. Bu durumda iş daha da kolaydır. Zira olasılık yoğunluk fonksiyonu gamma.pdf() fonksiyonu kullanılarak kolaylıkla elde edilebilir. Fonksiyonun parametreleri $a$ ve $scale$ ($b$) parametreleri ile verinin kendisidir.

In [15]:
import numpy as np
from scipy import stats as st
def olabilirlik_gamma(theta,veri):
    return np.sum(np.log(st.gamma(a=theta[0],scale=theta[1],loc=0).pdf(veri)))

Herhangi bir adımdaki değeri ve olasılığını kabul etmemiz için gerekli kabul kriteri daha önce tanımlandığı şekilde, ardışık iki adımdaki olabilirlik fonksiyonu değerlerinin oranına bağlı olarak aşağıdaki şekilde tanımlanır.

In [16]:
def kabul_kriteri(theta, theta_yeni):
    # Eger adimi atmakla elde edilecek dagilim daha olasiysa kabul et
    if theta_yeni > theta:
        return True
    # Degilse, 0-1 arasinda bir sayi uret 
    else:
        rastgele_sayi=np.random.uniform(0,1)
        # Hesapladigimiz olasiliklar logaritmik oldugu icin
        # karsilatirmayi yapmak uzere eksponansiyel al ve
        # urettigin sayi yeni adimin olasiligi ile mevcut adimin
        # olasiligi arasindaki farktan dusukse onu,
        # degilse yeni adim ile mevcut adim olasiklari farkini dondur.
        # Boylece olasiligi dusuk yeni adimlar daha az atilan adimlar olsun
        return (rastgele_sayi < (np.exp(theta_yeni - theta)))

Bundan sonra tıpkı Örnek 1. Basit Veri örneğinde olduğu gibi MCMC Metropolis-Hastings algoritması takip edilerek ve yukarıdaki fonksiyonlar temelinde hesaplanan olasılıklar karşılaştırılarak ardıl dağılıma ulaşılabilir.

In [17]:
# Metropolis-Hastings algoritmamiz
def metropolis_hastings(olabilirlik_fonk, oncul, gecis_fonksiyonu, baslangic_param, iterasyon_sayisi, \
                        veri, kabul_kriteri):
    # oablirlik_fonksiyonu(x,veri): verilen parametreler (a,b) dahilinde 
    # gamma dagiliminin veriyi turetme olasiligini hesaplayan fonksiyon
    # gecis_modeli(x): oneri fonksiyonu. Yaptigi tek is bir sonraki
    # adimi belirlemek uzere simetrik (Gauss) bir dagilimdan
    # her iki parametre icin (a,b) rastgele birer ornek secip dondurmek
    # baslangic_param: parametreler (a,b) icin baslangic degerleri
    # iterasyon_sayisi: iterasyon sayisi
    # gozlem_veri: Gunes leke sayilari
    # kabul_kriteri(x,x_yeni): adimi kabul edip etmemeye karar veren fonk.
    x = baslangic_param
    kabul = []
    red = []   
    for i in range(iterasyon_sayisi):
        # oneri dagilimiyla yeni dagilimi (x_yeni) belirlemek uzere
        # mevcut dagilimi gecis model fonksiyonuna gonderdik.
        x_yeni =  gecis_fonksiyonu(x)
        # Mevcut ve yeni dagilimlar icin olabilirlik fonksiyonu hesabi
        # Bu fonksiyon istenirse elle yazilan formul istenirse de
        # scipy.stats.norm olabilir.
        x_olabilirlik_degeri = olabilirlik_fonk(x,veri)
        x_yeni_olabilirlik_degeri = olabilirlik_fonk(x_yeni,veri)
        # yeni dagilimi kabul edip etmeyecegimize yonelik karar
        if (kabul_kriteri(x_olabilirlik_degeri + np.log(oncul(x)),x_yeni_olabilirlik_degeri+np.log(oncul(x_yeni)))):            
            # kabul edildiyse kabul edilenler listesine ekle
            x = x_yeni
            kabul.append(x_yeni)
        else:
            # degilse reddedilenler listesine ekle
            red.append(x_yeni)
    return np.array(kabul), np.array(red)

Örneğin sürekliği ve bağımsızlığını korumak açısından Metropolis-Hastings algoritması yukarıda verilmiştir ama daha önceki basit veri uygulamasındaki Metropolis-Hastings fonksiyonu da birebir kullanılabiir. Bundan sonra yapılması gereken aynı şekilde iterasyon sayısı (50000) ve $(a,b)$ parametreleri için birer başlangıç değeri $(4,10)$ belirledikten sonra algoritmayı çalıştırmaktır. Yine daha sonra, ısınma periyodu (burn-in) olarak belirleyip atacağımız ilk 50 adımda $a$ ve $b$ parametrelerinin nasıl değiştiğine bakalım.

In [18]:
# Algoritmayi calistir ve sonuclari incele
# Dagilim icin baslangic degerimiz 1000lik ornegin ortalamasi
# standart sapma icin baslangic degeri 0.1,
# iterasyon sayisi 50000
from matplotlib import pyplot as plt
%matplotlib inline
kabul, red = metropolis_hastings(olabilirlik_gamma,oncul,gecis_fonksiyonu,\
                                 [4,10], 50000, leke_sayisi, kabul_kriteri)
# Sonuclari inceleme:

plt.plot(kabul[:50,0], kabul[:50,1], label="Yuruyus")
plt.plot(kabul[:50,0], kabul[:50,1], 'b+', label='Kabul')
plt.plot(red[:50,0], red[:50,1], 'rx', label='Red')
plt.xlabel("a")
plt.ylabel("b")
plt.legend(loc="best")
plt.title("$\Gamma$ Dagiliminin Parametreleri (a,b) icin MCMC Ile Yapilan Orneklemenin Ilk 50 Adimi")
plt.show()

İlk 50 adım parametrelerin nereye doğru yakınsadığını göstermedi. Şimdi tüm adımlara bakalım.

In [19]:
from matplotlib import pyplot as plt
%matplotlib inline
plt.plot(kabul[:,0], kabul[:,1], label="Yuruyus")
plt.plot(kabul[:,0], kabul[:,1], 'b.', label='Kabul',alpha=0.3)
plt.plot(red[:,0], red[:,1], 'rx', label='Red',alpha=0.3)
plt.xlabel("a")
plt.ylabel("b")
plt.legend(loc = "best")
plt.title("$\Gamma$ Dagiliminin Parametreleri (a,b) icin MCMC Ile Yapilan Ornekleme (Tum Adimlar)")
plt.show()

Reddedilen parametreler de $a$ ve $b$ parametrelerinin sırasıyla $(0.8,1.2)$ ve $(60,100)$ değerleri arasında yoğunlaştığı için kabul edilen değerlerin nereye yakınsadığı tam görülememektedir. Reddedilen değerleri çizdirmeyedebiliriz ama pedagojik olması açısından tercihimiz bu yöntedir. Parametrelerin yakınsan değerlerini görmek açısından grafiğin ilgili bölümünü (son 50 değer) tekrar çizdirelim.

In [20]:
from matplotlib import pyplot as plt
%matplotlib inline
son_n = 50
plt.plot(kabul[-son_n:,0], kabul[-son_n:,1], label="Yuruyus")
plt.plot(kabul[-son_n:,0], kabul[-son_n:,1], 'b+', label='Kabul',alpha=0.5)
plt.plot(red[-son_n:,0], red[-son_n:,1], 'rx', label='Red',alpha=0.5)
plt.xlabel("a")
plt.ylabel("b")
plt.legend(loc="best")
plt.title("$\Gamma$ Dagiliminin Parametreleri (a,b) icin MCMC Ile Yapilan Orneklemenin Son 50 Adimi")
plt.show()

Bu da parametrelerin yakınsanan değerileri konusunda iyi bir fikir vermedi. Leke sayılarını temsil etmek üzere kullanılan $\Gamma$ dağılımı modelinin $a$ parametresini Metropolis-Hastings yöntemiyle atılan son %50 adımındaki değerleri ve bu değerlerin histogramına bakarak $a$ parametresinin değerini marjinalize ederek onun hakkında daha iyi bir fikir edinebiliriz.

In [21]:
# Yapilan isin basarimini test etmek iz grafigi (trace plot) ve
# ve parametrelerin histogramlari
# Isinma (burn-in) adimlarini ilk %25 olarak kabul edip atalim
from matplotlib import pyplot as plt
%matplotlib inline
son_075 = int(0.25*kabul.shape[0])
# 0.25 ile *0.75'ten sonrasi ayni olmali
hist_altlimit=int(-0.75*kabul.shape[0])
plt.plot(kabul[son_075:,0])
plt.title("a Parametresi icin Iz Grafigi (Trace Plot)")
plt.ylabel("a")
plt.xlabel("Iterasyon Sayisi")
plt.show()
plt.hist(kabul[hist_altlimit:,0], bins=20,density=True)
plt.ylabel("Normalize Frekans")
plt.xlabel("a")
plt.title("a Parametresinin Histogrami")
plt.tight_layout()
plt.show()

Şimdi artık $a$ parametresi konusunda daha iyi bir fikre sahibiz. Aynısını $b$ parametresi için de yapabiliriz.

In [22]:
from matplotlib import pyplot as plt
%matplotlib inline
plt.plot(kabul[son_075:,1])
plt.title("b Parametresi icin Iz Grafigi (Trace Plot)")
plt.ylabel("b")
plt.xlabel("Iterasyon Sayisi")
plt.show()
plt.hist(kabul[hist_altlimit:,1], bins=20,density=True)
plt.ylabel("Normalize Frekans")
plt.xlabel("b")
plt.title("b Parametresinin Histogrami")
plt.tight_layout()
plt.show()

$\Gamma$ dağılımının $a$ ve $b$ parametrelerini marjinalize ederek olasılık dağılımlarını gördük. Şimdi her iki parametre için ortak bir çözüm (olasılık dağılımı) ihtiyacımız olacağından bileşke olasılık dağılımına (ing. joint distribution) bakalım.

In [23]:
# a ve b'nin bileske olasilik dagilimlari (joint distribution)
from matplotlib import pyplot as plt
%matplotlib inline
# Hem a, hem de b parametresi icin x ve y eksenlerinde binler olustur
# x ekseni: a (0.8 ile 1.2 arasinda 30 bin yeterli
# a 1.0 civarinda oldugu ve cozunurluk acisindan da ilk %25'i attigimiz icin
# 50000*0.75 = 37500 adimda elde edilen degerleri 30'ar 30'ar birlestir.
# b 85 civarinda ve yine 30ar 30 ar birlestir.
xbin, ybin = np.linspace(0.8, 1.2, 30), np.linspace(75, 90, 30)
sayimlar, xekseni, yekseni, gorsel = plt.hist2d(
    kabul[hist_altlimit:, 0], kabul[hist_altlimit:, 1], normed=True, bins=[xbin, ybin])
plt.xlabel("a")
plt.ylabel("b")
plt.colorbar(gorsel)
plt.title("$a$ ve $b$ Parametreleri icin 2 boyutlu Histogram")
plt.show()

Bileşke dağılım $a$ ve $b$ için oluşturulan parametre uzayında en olası değerlerin sarı renkle en az olası değerlerin lacivertle renklendirildiği bir renk skalasıyla verilmiştir. En iyi çözümün $a$ için 0.97, $b$ iin 84 civarında olduğu görülmektedir. Klasik frekansçı yaklaşımdaki gibi parametreler ve onların standart hatalarını vermek yerine parametreler için bu şekilde bir olasılık dağılımı vermek Bayesian istatistiğin temel yaklaşımına daha uygundur. Buna karşın, her iki parametre için ortalama değer ve standart sapmalar da ilgili dizilerden alınıp ayrıca verilebilir.

Elde ettiğimiz ardıl olasılık dağılımını (posterior) orjinal dağılımla karşılaştırmak üzere, biraz da muhafazakar davranıp MCMC adımlarının bu kez %50'sini atarak kalan %50'sinde kabul edilen $a$ ve $b$ parametrelerinin ortalama ve standart sapmalarını alalım. Bu ortalama ve standart sapmaları kullanarak orjinal dağılımımızla aynı büyüklükte bir dağılım elde edelim ve orjinal dağılımla karşılaştıralım. Burada elde edilen ortalama ve standart sapmalar, $a$ ve $b$ parametreleri için en olası değerler ve standart sapmaları olarak kullanılabilir.

In [24]:
# Ardil dagilimin basarimini orjinal dagilimla karsilastirarark
# incelemek uzere histogram
# Bu kez daha muhafazakar davranip sadece son %50 adimi alalim
# ve bu %50 adimda bulunan parametrelerin ortalamasini ve
# standart sapmasini verelim.
ardil_altlimit = int(0.5*kabul.shape[0])
a = kabul[ardil_altlimit:, 0].mean()
b = kabul[ardil_altlimit:, 1].mean()
a_sigma = kabul[ardil_altlimit:, 0].std()
b_sigma = kabul[ardil_altlimit:, 1].std()
print("a = {:.2f} +/- {:.2f}, b = {:.2f} +/- {:.2f}".
      format(a, a_sigma, b, b_sigma))
a = 0.98 +/- 0.02, b = 84.00 +/- 2.58

Bu parametrelerle ($a$ ve $b$) tanımlı bir $\Gamma$ dağılımı oluşturmak üzere numpy.random.gamma() fonksiyonunu kullanalım. Sonra da gözlemsel veriden elde ettiğimiz orjinal leke sayıları (leke_sayisi) bu dağılımı birer histogramlarını çizerek karşılaştıralım.

In [25]:
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
def gamma_dagilimi(x, a, b): return np.random.gamma(a, b, x)
x = np.arange(leke_sayisi.shape[0])
orneklem_ardil = gamma_dagilimi(x.shape[0], a, b)
# Leke sayilarinin ardildan elde edilen tahmini dagilimi ile
# orjinal dagilimin karsilastirmasi
plt.hist(orneklem_ardil, bins=np.linspace(
    0, 500, 50), density=True, label="Tahmini dagilim")
plt.hist(leke_sayisi, bins=np.linspace(0, 500, 50),
         alpha=0.5, density=True, label="Orjnal dagilim")
plt.xlabel("Leke Sayisi")
plt.ylabel("Frekans")
plt.title("Ardil Dagilimdan Elde Edilen Tahmini ve Orjinal Gozlemler")
plt.legend(loc="best")
plt.show()

Görüldüğü gib bir $\Gamma$ fonksiyonundan hareketle MCMC kullanılarak $a$ ve $b$ parametreleri için elde edilen değerleri Bayesci olasılıklar bağlımında değerlendirilerek elde edilen dağılımların ortalama ve standart sapmalarından hareketle elde edilen leke sayıları orjinal dağılımla büyük benzerlik göstermektedir. Veriden haraeketle yapılan bu analiz, gelecekteki leke sayılarının tahmininde kullanılabilir. Doğal olarak bu tahminler yapılırken varsayım Güneş'in manyetik etkinliğinin çevrimsel değişimler dışında büyük bir değişim göstermediğidir. Ancak bu tür değişimler de yeni verilere yansıyacağından veri-yönelimli bir istatistik paradigması olan Bayesian istatistikle yağılan analizlerin doğasına uygun olarak $a$ ve $b$ parametrelerinin güncellenmesine ve değerlerine ilişkin "inancımızı" ifade eden olasılık dağılımlarının değişimine neden olacaktır.

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
  • Devinderjit Sivia, John Skilling, “Data Analysis: A Bayesian Tutorial” Oxford University Press, USA (2006)
  • Bayesian Statistics: A Comprehensive Course, Ox Education