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?¶
MCMC ve Bayesian Analiz Örnekleri¶
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.
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) $$# 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.
# 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.
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.
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.
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.
# 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()
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
# 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.
# 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()
sigmastd_075post = kabul[baslangic_indeksi:son_indeks,1].std()
print("Ardilin son %75'inden elde edilen mu={:.4f}, sigma={:.4f} +/- {:.4f}".\
format(mu_075post, sigma_075post[0],sigmastd_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()
Ö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.
# 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.
# 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.
# 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.
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.
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.
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.
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.
# 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.
# 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.
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.
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.
# 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.
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.
# 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.
# 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))
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.
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.
pymc3 ile Modelleme Örneği¶
Full-Bayesian bir analizi adım adım gördükten sonra, bu adımları "perde arkasında", bir miktar otomatize şekilde gerçekleştiren bir python modülü olan pymc3 ile bir modelleme alıştırması yapmak üzere bir ötegezegenin geçiş zamanları değişimi (Transit Timing Variations, TTV) grafiğini modellemeye çalışalım. pymc3 NUTS (No-U-Turn Sampling) gibi gradyent tabanlı MCMC algoritmaları ile hızlı ve yaklaşık Bayesian çıkarım için ADVI (Automatic Differentiation Variational Inference) kullanarak büyük veri setleri de dahil olmak üzere veri üzerinde Bayesian yaklaşımla efektif, parametrik olmayan modeller oluşturmak için yazılmış bir pakettir.
Paketin en önemli nesnelerinden biri olan $Model$ nesnesi büyük esneklik sağlar. Öncelikle Bayesian yaklaşımla fit etmeyi planladığımız TTV verisini açalım.
import pandas as pd
df = pd.read_csv("veri/ttv_data.csv")
df2 = pd.DataFrame({"epoch":df['epoch'],'TTV':df['TTV'],'Tc_err':df['sigma_Tc']})
df2.head()
Bir TTV verisi bir ötegezegenin referans bir geçiş zamanı ($T_0$) ve yörünge dönemine ($P_{orb}$) göre bu referans zaman sonrası kaçıncı geçişinin gözlendiği ($epoch$), gözlenen zamanın, referans $T_0$ ve $P_{orb}$ ikilisi kullanılarak hesaplanan zamana göre ne kadar farkla gerçekleştiği ($TTV$) ve geçiş zamanının ölçüm hatası ($Tc_{err}$) parametrelerinden oluşur ve $epoch$ parametresi çevrimi sayısını gösterdiğinden birimsiz iken diğerleri "genellikle" gün ya da dakika birimlerindedir.
Bu veriye Bayesian istatistik paradigması bağlamında bir doğru uyumlaması yapmak üzere doğru modelimizi oluşturalım. MCMC analizlerinde kullanılan geçiş algoritmaları (pymc3 'te NUTS) bir başlangıç noktasına ihtiyaç duyarlar. Bu başlangıç noktasını belirlemek üzere lineer ya da lineer olmayan en küçük kareler yöntemine dayalı uyumlamaya olanak sağlayan herhangi bir paket kullanılabilir. Bu örnekte scipy.optimization.curve_fit fonksiyonu bu amaçla kullanılmıştır.
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
# Dogru modeli
def dogru(x, m, n):
return m*x + n
nlls_model = curve_fit(dogru,df2['epoch'], df2['TTV'])
# Parametre degerleri
a = nlls_model[0][0]
b = nlls_model[0][1]
# Kovaryans matrisi ve parametre hatalari
cov=nlls_model[1]
a_err=np.sqrt(cov[0][0])
b_err=np.sqrt(cov[1][1])
print('a = {:g} +/- {:g}'.format(a,a_err))
print('b = {:g} +/- {:g}'.format(b,b_err))
plt.errorbar(df2['epoch'], df2['TTV'],df2['Tc_err'],fmt='.k')
E=np.linspace(min(df2['epoch']),max(df2['epoch']),1000)
plt.plot(E,dogru(E,a,b),"b-")
plt.xlabel('Epoch')
plt.ylabel('TTV (day)')
Şimdi bu parametreleri başlangıç parametresi olarak kullanmak suretiyle pymc3 modelini oluşturalım. Model parametelerinin baslangıç değerlerini lineer olmayan en küçük kareler uyumlamasından alalım ve bir bilgimiz olduğu icin bunlari normal bir dağılımın en olası değerleri, hatalarını bu dağılımların standart sapması olarak alallım.
import pymc3 as pm
import arviz as az
dogru_modeli = pm.Model()
with dogru_modeli:
# Model parametelerinin baslangic degerleri
# sigmalari biraz genis tutabiliriz.
a = pm.Normal('a', mu= 9.1e-08, sigma= 1.0e-07)
b = pm.Normal('b', mu= -3.6e-05, sigma= 1.0e-04)
# uyumlanan model
model = a*df2['epoch'] + b
# Bu model baglaminda bu verinin olusmus olma olasiligini
# bir normal dagilimdan belirleyelim
likelihood = pm.Normal('y', mu= model, sigma=df2['Tc_err'],
observed=df2['TTV'])
# Orneklendirme (sampling)
trace = pm.sample(chains=10, draws=5000,init='adapt_diag',cores=2,tune=1000,\
return_inferencedata=False,\
discard_tuned_samples=True,start=None,progressbar=True)
az.plot_trace(trace);
Görüldüğü gibi her iki parametrenin de ardıl olasılık dağılımları (solda) normal dağılıma yakındır ve her bir yürüyücünün attığı adımlar boyunca (analize anılmayan ilk 1000 adım hariç bu normal dağılımların ortalama değeri etrafında küçük sayılabilecek birer saçılma bulunmaktadır (sağda)
Sonuçları bir tablo halinde görüntülemek üzere arviz.trace fonksiyonunun summary metodunu kullanalım.
az.summary(trace, round_to=8)
Şimdi sonuçlarımızı grafik üzerinde görelim.
### Sonuclarimizi bir grafige aktaralim
plt.errorbar(df2['epoch'], df2['TTV'], df2['Tc_err'],fmt='o', label='data')
for a, b in zip(trace['a'][-1000:], trace['b'][-1000:]):
plt.plot(df2['epoch'], a*df2['epoch']+b, c='red', alpha=0.1)
Bulgularımızı raporlayalım:
print("Ilk 1000 adim atilmistir (burn-in)")
print("--"*50)
median_a = np.median(trace['a'])
perc16_a = np.percentile(trace['a'],16)
perc84_a = np.percentile(trace['a'],84)
median_b = np.median(trace['b'])
perc16_b = np.percentile(trace['b'],16)
perc84_b = np.percentile(trace['b'],84)
print("a icin medyan degeri {:g}".format(median_a))
print("a icin 16. persantil {:g}".format(perc16_a))
print("a icin 84. persantil {:g}".format(perc84_a))
print("a = {:g}^(+{:g})_(-{:g})".format(median_a,np.abs(median_a-perc16_a),np.abs(perc84_a-median_a)))
print("b icin medyan degeri {:g}".format(median_b))
print("b icin 16. persantil {:g}".format(perc16_b))
print("b icin 84. persantil {:g}".format(perc84_b))
print("b = {:g}^(+{:g})_(-{:g})".format(median_b,np.abs(median_b-perc16_b),np.abs(perc84_b-median_b)))
print("a icin ortalama deger {:g}".format(np.average(trace['a'])))
print("b icin ortalama deger {:g}".format(np.average(trace['b'])))
print("a icin ortalama degerin standart sapmasi {:g}".format(np.std(trace['a'])))
print("b icin ortalama degerin standart sapmasi".format(np.std(trace['b'])))
Köşe diyagramını (ing. corner plot) çizdirelim.
labels = ["dP","dT"]
import corner
samples = np.vstack([trace[k] for k in ["a", "b"]]).T
corner.corner(samples, bins=30, max_n_ticks=1, \
smooth=True, smooth1d=True, show_titles=False, \
scale_hist=True,title_fmt=".8f",range=None, weights=None, \
color="k",truths=None,labels=labels,label_kwargs={"fontsize": 20}, \
unit_kwargs={"fontsize":1})
Model başarımına ilişkin istatistikleri hesaplayalım.
bayes_lin_model = a*df2['epoch']+b
m=2 # lineer model
toplam = 0
n = len(df2['TTV'])
for i in range (0,n):
fark = (bayes_lin_model[i] - df2['TTV'][i])/df2['Tc_err'][i]
toplam += fark**2
chi2_nu = toplam/(n-m) #dividing summation by N-m
print ("chi^2_nu: {:g}".format(chi2_nu))
BIC=chi2_nu+(m * np.log(n))
AIC = chi2_nu + (2 * m)
print('BIC: {:g}'.format(BIC))
print('AIC: {:g}'.format(AIC))
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
- “Data Analysis: A Bayesian Tutorial”, Devinderjit Sivia, John Skilling, Oxford University Press, USA, 2006
- Bayesian Statistics: A Comprehensive Course, Ox Education
- pymc3 Dokümantasyonu
- "Bayesian Analysis in Python", Osvaldo Martin, Packt Publishing, 2nd ed., 2018