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?¶
Dağılımlar¶
- Olasılık Dağılım Fonksiyonlarının Özellikleri
- Python'da Dağılımlar
- Tekdüze Dağılımlar
- Normal Dağılım ve Gauss Fonksiyonu
- Log Normal Dağılım
- Binom Dağılımı
- Bernoulli Dağılımı
- Poisson Dağılımı
- Cauchy Dağılımı ve Voigt Profili
- Gamma Dağılımı
- Üstel Dağılım
- Weibull Dağılımı
- Kuvvet Yasası Dağılımları
- Pareto Dağılımı
- Ödev 1
- Kaynaklar
Olasılık Dağılım Fonksiyonlarının Özellikleri¶
Rastgele Değişken (ing. random variable): Olası değerleri rastgele bir olayın sayısal sonuçları olan bir değişkendir. İki tür rastgele değişken vardır, süreksiz (kesikli) (ing. discrete) ve sürekli (ing. continuous) değişkenler.
Süreksiz bir rastgele değişken, sadece sayılabilir değerler alabilen ve böylece nicelleştirilebilen bir değişkendir. Örneğin, adil bir zarın her atılıştaki değeri rastgele bir $x$ değişkeni olarak tanımlanabilir. Bu durumda $x$, $[1,2,3,4,5,6]$ değerlerini alabilir. Bu nedenle de $x$, süreksiz rastgele bir değişkendir.
Süreksiz rastgele bir değişkenin olasılık dağılımı, olası değerlerinin her biri ile ilişkili olasılıkların bir listesidir. Buna olasılık fonksiyonu ya da daha doğru bir ifadeyle Olasılık Kütle Fonksiyonu (ing. Probability Mass Function, PMF) denir.
Rastgele bir $X$ değişkeninin k farklı değer alabileceğini varsayalım. $X = x_i$ 'ye karşılık gelen $P (X = x_i) = p_i$ olasılıkları aşağıdaki koşulları sağlar:
- $0 \lt p_i \lt 1$, $i = 1, 2, 3, ...$
- $p_1 + p_2 + ... + p_k = 1$
Sürekli rastgele değişken, sonsuz sayıda olası değeri alabilen değişkendir. Örneğin, bir sınıftaki öğrencilerin boyları rastgele bir $x$ değişkeni olarak tanımlandığında sürekli bir rastgele değişkendir Sürekli rastgele değişken bir değerler aralığı içinde tanımlandığından, o aralıkta olma olasılığı, olasılık yoğunlu fonksiyonunu ifade eden eğrinin (veya integralin) altındaki alanla temsil edilir.
Olasılık Yoğunluk Fonksiyonu (ing. Probability Density Function, PDF olarak bilinen sürekli rasgele değişkenin olasılık dağılımı, sürekli değerleri alan fonksiyonlardır. Herhangi bir tek değeri gözlemleme olasılığı, 0'a eşittir, çünkü rastgele değişken tarafından alınabilecek değerlerin sayısı sonsuzdur ($p(x_i) = 0$). Olasılık yoğunluk fonksiyonu için $P (a \lt x \lt b) = p$ olasılıkları aşağıdaki koşulları sağlar:
Olasılık her zaman pozitiftir ($x$'in bütün değerleri için $p \gt 0$).
$pdf(x)$'in altında kalan toplam 1'dir.
Ayrıca, tüm rastgele değişkenlerin (sürekli ya da değil) birikimli (kümülatif) dağılım karşılğı da vardır. Birikimli dağılım fonksiyonu, her $x$ değeri için rastgele değişken $X$'in $x$'den küçük veya ona eşit olma olasılığını veren bir fonksiyondur ($P(X \le x) = p(x))$. Kesikli rastgele bir değişken için olasılık dağılımları toplanarak birikimli (kümülatif) dağılım fonksiyonu bulunur.
Python'da Dağılımlar¶
Python'da tüm olasılık dağılım fonksiyonlarına ilişkin fonksiyonların yer aldığı modüller bulunmaktadır. Bu modüllerden en sık kullanılanları scipy.stats ve np.random paketleridir. Bu modüllerdeki fonksiyonlarla ilgilenilen dağılımın yapısına uygun örnekler ya da dağılım fonksiyonları oluşturularak incelenebilir. Çizim için matplotlib.pyplot fonksiyonları hist, bar ve plot en sık kullanılanlardır. Ayrıca seaborn modülü fonksiyonları da kolay kullanılabilmeleri nedeniyle önemli bir alternatif oluşturmaktadır.
scipy.stats paketindeki neredeyse tüm dağılım fonksiyonlarında bulunan bazı argüman ve parametreler, Normal Dağılım örneğinde (stats.norm) aşağıda listelenmeye çalışılmıştır.
rvs(loc=0, scale=1, size=1, random_state=None): Söz konusu dağılımdan bir örnek oluşturmak için kullanılır.
pdf(x, loc=0, scale=1): Olasılık yoğunluk fonksiyonu oluşturmak için kullanılır.
logpdf(x, loc=0, scale=1): Olasılk yoğunluk fonksiyonun logaritması (doğal logaritma, ln)'nı verir.
cdf(x, loc=0, scale=1): İlgilenilen dağılım türü temelinde bir birikimli (kümülatif) yoğunluk fonksiyonu oluşturmak için kullanılır.
logcdf(x, loc=0, scale=1): Birikimli yoğunluk fonksiyonunun logaritmasını (ln) verir.
sf(x, loc=0, scale=1): Kümülatif fonksiyon $P(X \le x)$ olasılığını verirken $sf$ (ing. survival function) $P(X \gt x) = 1 - (P(X \le x)$ olasılığını verir.
logsf(x, loc=0, scale=1): sf fonksiyonunun logaritmasıdır.
ppf(q, loc=0, scale=1): (ing. percent point function), birikimli yoğunluk fonksiyonunda elde edilen olasılığın karşılık geldiği $X = x$ değeridir. Bu anlamda birikimli yoğunluk fonksiyonunun tersidir. Birikimli yoğunluk fonksiyonu $X = x$ değeri argüman olarak sağlandığında $P(X \ge x)$ olasılığını verirken, $ppf$ verilen $P(X \ge x)$ olasılık değeri için $X = x$ değerini verir.
isf(q, loc=0, scale=1): (ing. inverse survival function), ppf'e benzer şekilde $1 - cdf$ olasılığının karşılığı olan $X = x$ değerini verir.
moment(n, loc=0, scale=1): Dağılımın n. momentini verir. Bir olasılık dağılımının n = 1. momenti beklenen değer (ortalama), n = 2. momenti varyans, n = 3. momenti çarpıklık ve n = 4. momenti ise basıklığını verir. Daha fazla bilgi için
stats(loc=0, scale=1, moments=’mv’): Dağılıma ilişkin momentleri verir. İstenen momentler
momentsparametresine baş harfiyle sağlanır. (Ortalama: Mean(‘m’), Varyans: variance(‘v’), Çarpıklık: skew(‘s’), ve/veya kurtosis(‘k’)).entropy(loc=0, scale=1): Dağılımın entropisini verir.
fit(data, loc=0, scale=1):
dataparametresi üzerinden verilen bir veriye ilgili dağılımı uyumlar (fit eder) ve bu uyumalamanın dağılımın türüne göre ilgili parametrelerini döndürür.expect(func, args=(), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds): Tek değişkenli bir fonksiyonun bir dağılıma göre beklenen değerini döndüren fonksiyondur.
median(loc=0, scale=1): Dağılımın ortancasını (medyan değerini) döndürür.
mean(loc=0, scale=1): Dağılımın ortalama değerini döndürür.
var(loc=0, scale=1): Dağılımın varyansını ($\sigma^2$) döndürür.
std(loc=0, scale=1): Dağılımın standart sapmasını ($\sigma$) döndürür.
interval(alpha, loc=0, scale=1): Tüm dağılımın $\alpha$ yüzdelik bölümüne denk gelen aralığının uç noktalarını verir.
Tekdüze Dağılımlar¶
Bir değişkenin tekdüze dağılımı onun bir aralık içerisindeki değerlerin tamamının aynı olasılıkla alabileceği sürekli bir dağılımdır.
numpy.random fonksiyonlarından uniform ile tekdüze bir dağılımdan giderek daha fazla sayıda örnek seçilerek tekdüze bir dağılımın nasıl oluştuğunu örnekleyen bir dizi kod aşağıda verilmiştir. Ayrıca tekdüze dağılımın her aşamada bir normal fonksiyona göre ne kadar daha çarpık (ing. skewness) ve basık olduğuna ilişkin değerler (dağılımın $skew$ ve $kurtosis$ parametreleriyle) de verilmiştir. Tekdüze bir dağılımın simetrik olması beklendiği için örnek sayısı arttırıldıkça çarpıklığının 0'a, ortancaya (medyan), basıklığı gösteren Kurtosis değerinin ise -1.2'ye yaklaşması beklenmelidir.
import numpy as np
from matplotlib import pyplot as plt
from scipy import stats as st
%matplotlib inline
x = np.random.uniform(-5,5,20)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=10)
plt.show()
x = np.random.uniform(-5,5,100)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=20)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
print('Normal bir dagilima gore basiklik: ', st.kurtosis(x))
print('Normal bir dagilima gore carpiklik: ',st.skew(x))
x = np.random.uniform(-5,5,1000)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=20)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
print('Normal bir dagilima gore basiklik: ', st.kurtosis(x))
print('Normal bir dagilima gore carpiklik: ',st.skew(x))
x = np.random.uniform(-5,5,100000)
print(x.mean(), x.std(), np.median(x))
plt.hist(x, bins=20)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
print('Normal bir dagilima gore basiklik: ', st.kurtosis(x))
print('Normal bir dagilima gore carpiklik: ',st.skew(x))
Tekdüze bir dağılımdan bir örneklemi scipy.stats.uniform fonksiyonlarından rvs ile de oluşturabilirsiniz. Tüm dağılımlar için bu fonksiyon bir örnek oluşturmak için kullanılır. Tekdüze dağılım için loc parametresi dağılımın başlangıç, scale parametresi son değerini, size parametresi için örneklem büyüklüğünü vermek için kullanılan argümanlarıdır.
# scipy.stats.uniform
x = st.uniform.rvs(loc=0,scale=2.5,size=10000)
plt.hist(x, bins=20)
plt.show()
Dağılımı bir sıklık dağılımı olarak göstermek yerine olasılık kütle fonksiyonu (PMF) cinsinden de gösterebilirsiniz. Bu durumda plt.hist fonkiyonunun density parametresini $True$ değerine ayarlamanız yeter.
N = 100000
x = np.random.uniform(-5,5,N)
print(x.mean(), x.std(), np.median(x))
plt.xlabel('x')
plt.ylabel('PMF')
plt.hist(x, bins=20, density=True)
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.show()
Normal Dağılım ve Gauss Fonksiyonu¶
Gauss dağılımı olarak da bilinen Normal Dağılım, deneysel bütün bilimlerde yeri bulunan sürekli bir olasılık dağılımıdır. Normal dağılım, ortalama, $\mu$ ve standart sapma, $\sigma$) ile tanımlanan çan şeklindeki yoğunluk eğrisine sahiptir. Yoğunluk eğrisi (ing. Probability Density Function, PDF), merkezi ortalamaya yakın, yayılımı standart sapması ile belirlenen; ortalamaya yakın verilerin, ortalamadan çok uzaktaki verilere göre daha sık olduğunu gösteren simetrik bir fonksiyondur. Belirli bir $x$ noktasında ortalama $\mu$ ve standart sapma $\sigma$ ile normal yoğunluk eğrisinin olasılık dağılım fonksiyonu şu şekilde verilir:
$$ p(x)~dx = \frac{1}{\sigma\sqrt{2\pi}} e^{\frac{-1}{2} (\frac{x - \mu}{\sigma})^2}~dx $$def gauss(x,mu,sigma):
a = 1/(np.sqrt(2*np.pi)*sigma)
b = np.exp(-0.5*((x - mu)/sigma)**2)
return a*b
x = st.norm.rvs(loc=0,scale=1,size=100000)
# Ayni seyi numpy.random ile de yapabilirsiniz.
#x = np.random.randn(1000)
print(x.mean(), x.std(), np.median(x))
plt.show()
x2 = np.linspace(-5,5,100000)
fx2 = gauss(x2, x.mean(), x.std())
plt.hist(x, bins=50, density=True)
plt.plot(x2,fx2,'r-')
plt.axvline(x = x.mean(), ls='--', c='r')
plt.axvline(x = np.median(x), ls='--', c='g')
plt.ylabel("PDF")
plt.xlabel("x")
plt.show()
print('Normal bir dagilima gore basiklik: ', st.kurtosis(x))
print('Normal bir dagilima gore carpiklik: ',st.skew(x))
Standart Normal Dağılım¶
Standart normal dağılım ($\mu = 0$, $\sigma = 1$) için verilen normal dağılımdır. Tüm normal dağılımlar standart normal dağılıma dönüştürülebilir. Bunun için o dağılım için verilen ortalama değeri ($\mu$) ve standart sapmayı ($\sigma$) kulanarak
$$ z = \frac{x - \mu}{\sigma} $$değişken değiştirmesini yapmak gerekir.
Bu durumda olasılık yoğunluk fonksiyonu
$$ p(z)~dz = \frac{1}{\sigma\sqrt{2\pi}} e^{\frac{z^2}{2}}~dz$$şeklinde ifade edilir.
Aslında normal dağılımı tanımlarken oluşturduğumuz örneklem ortalama değerin 0, standart sapmasının 1 olduğu bir dağılım kullanarak standart normal dağılım şeklinde bir tanımlama yapmış olduk. Şimdi bu dağılımın bazı özelliklerine bakalım.
Herhangi bir $x$ değişkeninin ortalamadan z = $\pm \Delta x$ kadar sapması olasılığı scipy.stats.norm paketi kümülatif dağılım fonksiyonu (ing. cumulative distribution function) cdf kullanılarak elde edilebilir.
import scipy.stats as st
z = 1.645
pz = st.norm.cdf(z)
print("Standart normal dagilimda z = +{:.2f} sigma sapma olasiligi : {:.2f}". \
format(z,pz))
Fonksiyon $x$'in $-\infty$'dan $\Delta x = +1.65\sigma$ 'ya kadarki aralıkta olma olasılığını verir (standart normal dağılımın bu değere kadar sol tarafının altında kalan alan). Eğer x'in $\Delta x = \pm1.65\sigma$ aralığında olma olasılığı isteniyorsa bu kez bu değerin
z1 = -1.645
z2 = 1.645
pz_z1_z2 = st.norm.cdf(z2) - st.norm.cdf(z1)
print("Standart normal dagilimda z'nin = +/-{:.2f} sigma dahilinde olma olasiligi : {:.2f}". \
format(z2,pz_z1_z2))
Dolayısıyla normal dağılımı temsil eden eğrinin en yüksek eğime sahip tanjantlarının eğriyi ve x eksenini kestiği noktalar; sırasıyla 1 ve 2 $\sigma$ verilerek öğrenilebilir.
z1 = -1.0
z2 = 1.0
pz_z1_z2 = st.norm.cdf(z2) - st.norm.cdf(z1)
print("Standart normal dagilimda z'nin = +/-{:.2f} sigma dahilinde olma olasiligi : {:.2f}". \
format(z2,pz_z1_z2))
z1 = -2.0
z2 = 2.0
pz_z1_z2 = st.norm.cdf(z2) - st.norm.cdf(z1)
print("Standart normal dagilimda x'in = +/-{:.2f} sigma dahilinde olma olasiligi : {:.2f}". \
format(z2,pz_z1_z2))
Tersten giderek verilen bir olasılığın merkezden kaç standart sapma uzaklığa denk geldiği (yani $z$ değeri) de bulunabilir. Bunun için ppf (percent point function) fonkiyonu kullanılır. Fonksiyonun yine $+\Delta x$'e kadar olasılık verdiği (eğrinin sol tarafının altında kalan) akılda tutulmalıdır.
pz = 0.95
z = st.norm.ppf(pz)
print("Standart normal dagilimda {:.2f} olasiligin karsilik geldigi sapma degeri : {:.2f}". \
format(pz,z))
İki taraflı olasılık isteniyorsa, $pz$ argümanı buna uygun olarak verilmelidir.
pz = 0.05
z = st.norm.ppf(pz)
print("Standart normal dagilimda +\-{:.2f} olasiligin karsilik geldigi sapma degeri : +/-{:.2f}". \
format(pz,z))
Aynı amaçla daha kullanışlı olan interval fonksiyonu da kullanılabilir. norm.interval fonksiyonu argüman verilen bir olasılığın $x$'in hangi z = $\pm \Delta x$ değerleri arasında olma olasılığı olduğunu verir. Bir başka deyişle, verilen olasılığa karşılık gelen alanın uç noktaları olan $z$ değerlerini verir.
pz = 0.90
z = st.norm.interval(pz)
print(z)
print("Standart normal dagilimda +\-{:.2f} olasiligin karsilik geldigi sapma degeri : +/-{:.2f}". \
format(pz,z[1]))
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Polygon
fig, ax = plt.subplots()
x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 1000)
ax.plot(x, st.norm.pdf(x), 'r-', lw=2, alpha=0.6, label='norm pdf')
xi = st.norm.ppf(0.001)
xf = st.norm.ppf(0.05)
ix = np.linspace(xi,xf)
iy = st.norm.pdf(ix)
verts = [(xi, 0), *zip(ix, iy), (xf, 0)]
poly = Polygon(verts, facecolor='b')
ax.add_patch(poly)
ax.set_xlabel("z")
ax.set_ylabel("pdf")
plt.arrow(-2.75, 0.10, 0.5, -0.05, capstyle="round",
head_width=0.02, head_length=0.05)
plt.text(-2.6,0.115,"p(z)")
plt.show()
Bu alanın değerini birikmli (kümülatif) dağılım fonksiyonu cdf üzerinden görebiliriz. Bir başka deyişle, kümülatif dağılım fonksiyonu, normal dağılım fonksiyonunun altında kalan alanın z değeri ile değişimini verir. $z \rightarrow \infty$ iken $cdf \rightarrow 1.0$ 'dır.
from matplotlib.patches import Polygon
fig, ax = plt.subplots()
x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 1000)
plt.plot(x, st.norm.cdf(x), 'r-', lw=2, alpha=0.6, label='norm pdf')
plt.axvline(x=st.norm.ppf(0.05))
plt.axhline(y=st.norm.cdf(st.norm.ppf(0.05)))
x_degeri = "z = {:.3f}".format(st.norm.ppf(0.05))
y_degeri = "p(z) = {:.3f}".format(st.norm.cdf(st.norm.ppf(0.05)))
plt.text(1, 0.1, y_degeri)
plt.text(-2.8,0.15, x_degeri)
plt.xlabel('z')
plt.ylabel('cdf')
plt.show()
st.norm fonksiyonları sadece standart normal dağılım için değil, herhangi bir normal dağılım için de kullanılabilir. Örneğin RR Lyrae yıldızlarının zonklama dönemleri dağılımından hareketle bir RR Lyrae yıldızının zonklama döneminin 17 saatten fazla olma olasılığını bulmak isteyebiliriz. RR Lyrae yıldızları için ortalama dönemin $\mu = 14$ saat, standart sapmasının ise $\sigma = 5$ saat olduğunu belirlemiş ve zonklama döneminin normal dağıldığını varsaymıştık.
Port = 14
Pstd = 5
x = 17
print("Bir RR Lyrae yildizi icin zonklama doneminin {:d} saatten uzun olma olasiligi: {:.2f}".\
format(x, 1 - st.norm.cdf(x, loc=Port, scale=Pstd)))
x = 3
print("Bir RR Lyrae yildizi icin zonklama doneminin {:d} saatten kisa olma olasiligi: {:.2f}".\
format(x, st.norm.cdf(x, loc=Port, scale=Pstd)))
Yukarıdaki ilk örnekte ($P_{zonk} > 17$) görüldüğü gibi cdf normal dağılım eğrisinin sağ tarafında kalan alanı verdiği için belirli bir değerden büyük alanı (olasılığı) hesaplamak için $1 - cdf$ kullanmak gerekebilir. Bunun yerine sf fonksiyonu (survival function) da kullanılabilir.
x = 17
print("Bir RR Lyrae yildizi icin zonklama doneminin {:d} saatten uzun olma olasiligi: {:.2f}".\
format(x, st.norm.sf(x, loc=Port, scale=Pstd)))
Aynı şekilde olasılıklardan hareketle zonkalama dönemleri de bulunabilir. Bu fonksiyona daha çok güven aralığı (ing. confidence level) oluşturmak için başvurulur.
px = 0.95
aralik = st.norm.interval(px,loc=Port, scale=Pstd)
print("Secilen bir RR Lyr yildizi %{:.1f} ihtimalle {:.2f} ile {:.2f} saatleri arasinda bir donemle zonklar."
.format(px*100, aralik[0], aralik[1]))
Log-Normal Dağılım¶
Bir değişimin logaritmasının normal dağılım göstermesi durumunda olasılık yoğunluk fonksiyonu log-normal dağılım fonksiyonu ile temsil edilir. Bir başka deyişle, bir $x$ değişkeni normal dağılıma uygun değerler alıyorsa, $e^x$ log-normal dağılıma uygun değerler alır. log-normal dağılım aşağıdaki şekilde fiade edilir.
$$ p(x)dx = \frac{1}{x~\sigma \sqrt{2\pi}}~exp(-\frac{(ln(x) - \mu)^2}{2\sigma^2})~dx$$Birden fazla pozitif değişkene, o değişkenlerin çarpımıyla bağlı bir değişkenin log-normal dağılım göstermesi beklenir. Eğer $y = exp(x)$ log-normal dağılıyorsa, bu dağılımın ortalama değeri $\bar{y} = exp(\mu + \sigma^2 / 2)$, medyanı $exp(\mu)$, modu ise $exp(\mu - \sigma^2)$ ile verilir.
Örnek olarak $\mu = 2$, $\sigma = 1$ parametreleriyle bir log-normal dağılım varsayalım.
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Port = 2.0
Pstd = 1.0
# Ortalamasi Port, standart sapmasi Pstd
# olan log-normal bir dagilimdan alinmis
# 1000 rastgele ornekle bir histogram olusturalim
# Kodunuzu her calistirdiginizda ayni dagilimi elde etmek icin np.random.seed
# ya da np.random.RandomState fonksiyonlarini kullanabilirsiniz
np.random.RandomState(123)
s = np.random.lognormal(Port, Pstd, 1000)
count, bins, ignored = plt.hist(s, 50, density=True, align='mid')
# histogramin sinirlari dahilindeki 1000 x
# degerini log-normal dagilimin olasilik yogunluk
# fonksiyonunda (pdf) yerine koyalim
# ve histogramin uzerine cizdirelim
x = np.linspace(min(bins), max(bins), 1000)
pdf = 1. / (x * Pstd * np.sqrt(2 * np.pi)) *\
(np.exp(-(np.log(x) - Port)**2 /(2 * Pstd**2)))
plt.plot(x,pdf,'r-')
plt.ylabel('PDF')
plt.axis('tight')
plt.show()
Olasılık yoğunluk fonksiyonunu, dağılımın tanımında (ve kodun içinde $pdf$ değişkeni ile) verilen ifadeyle sağlayabileceğimiz gibi $scipy.stats$ fonksiyonlarından $lognorm$ 'un $pdf$ fonksiyonunu da kullanabiliriz.
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats as st
%matplotlib inline
Port = 2.0
Pstd = 1.0
# Ortalamasi Port, standart sapmasi Pstd
# olan log-normal bir dagilimdan alinmis
# 1000 rastgele ornekle bir histogram olusturalim
#s = np.random.lognormal(Port, Pstd, 1000)
# Isterseniz ayni orneklemi sekil parametresini x'in standart sapmasi
# scale parametresini x'in ortalama degerinin eksponansiyeli olarak
# vermek kosuluyla st.lognorm.rvs fonksiyonu ile de olusturabilirsiniz.
s = st.lognorm.rvs(s=Pstd, scale=np.exp(Port), size=1000)
count, bins, ignored = plt.hist(s, 100, density=True, align='mid')
# histogramin sinirlari dahilindeki 1000 x
# degerini log-normal dagilimin olasilik yogunluk
# fonksiyonunda (pdf) yerine koyalim
# ve histogramin uzerine cizdirelim
# Bu kez scipy.stats.lognorm kullanalim
x = np.linspace(min(bins), max(bins), 1000)
pdf = st.lognorm.pdf(x,s=Pstd,scale=np.exp(Port))
plt.plot(x,pdf,'r-')
plt.ylabel('PDF')
plt.axis('tight')
Görüldüğü üzere st.lognorm.pdf diğer olasılık yoğunluk fonksiyonlarına göre bir miktar farklı çalışmaktadır. Bir log-normal dağılımın iki parametresi ($\mu$ ve $\sigma$) olmakla birlikte scipy.stats paketindeki diğer dağılım fonksiyonlarında loc dağılımın x-ekseninde kayma, scale parametresi ise ölçeğini verirken; scipy.stats.lognormal fonksiyonlarında $\sigma$, şekil parametresi; $e^{\mu}$, ölçek parametresi olarak verilir. Zira scipy.stats.lognormal fonksiyonları için olasılık yoğunluk fonksiyonu,
olarak tanımlanır. x'in logaritması değil kendisi $\mu$ ve $\sigma$ parametreleri ile normal dağılsaydı;
$$ y = exp(\mu + \sigma ~x) $$$$ y = exp(\mu) ~ exp(x)^\sigma $$$$ y = exp(\mu) ~ exp(x-x_0)^{\sigma} $$olurdu. Bu ifadeyi
$$ Y = m ~ (X - X_0)^n $$şeklinde yazarsak, $m = exp(\mu) = e^{\mu}$ ölçek parametresi, $n = \sigma$ şekil parametresi $(x - x_0) = e^{X - X_0}$ da $loc = X_0 = ln(x_0)$'dır.
Şimdi ortalamayı sabit tutup standart sapmayı değiştirdiğimiz vakit lognormal dağılımın normal dağılıma nasıl yaklaştığını inceleyelim.
mu = 1.0
sigma = [0.25,0.5,0.75,1.0]
x = np.linspace(0.01, 5.0, 1000)
for sig in sigma:
pdf = (np.exp(-(np.log(x) - mu)**2 / (2 * sig**2))/\
(x * sig * np.sqrt(2 * np.pi)))
leg = "$\sigma$= {:.2f}".format(sig)
plt.plot(x,pdf, label=leg)
plt.xlabel('x')
plt.ylabel('PDF')
plt.grid(True)
plt.legend(loc="upper left")
plt.show()
Aynı grafiği $scipy.stats.lognorm$ fonksiyonlarından $pdf$ ile elde etmeye çalışalım.
mu = 1.0
sigma = [0.25,0.5,0.75,1.0]
x = np.linspace(0.01, 5.0, 1000)
for sig in sigma:
pdf = st.lognorm.pdf(x,s=sig,scale=np.exp(mu))
leg = "$\sigma$= {:.2f}".format(sig)
plt.plot(x,pdf, label=leg)
plt.xlabel('x')
plt.ylabel('PDF')
plt.grid(True)
plt.legend(loc="upper left")
plt.show()
Binom Dağılımı¶
Başarı ya da başarısızlık, kazanç ya da kayıp, "yazı" ya da "tura" gibi sadece iki sonucun mümkün olduğu ve tüm denemeler için başarı (p) ve başarısızlık olasılığının (1-p) toplamının 1 olduğu süreksiz dağılıma Binom Dağılımı denir. Deneyin tekrar sayısının 1 olması durumunda elde edilen özel dağılıma ise Bernoulli Dağılımı adı verilir. Bununla birlikte, sonuçların eşit derecede olası olması gerekmez ve her deneme birbirinden bağımsızdır. Binom dağılımının parametreleri n ve p'dir; burada n, toplam deneme sayısıdır; p, her denemede başarı olasılığıdır. Süreksiz bir dağılım olduğu için olasılık kütle fonksiyonu (PMF) şu şekilde verilir:
$$ f(k, n, p) = P (k, n, p) = {n \choose k} p^k (1 - p)^{n-k} $$Burada $n \choose k$, n'in k'lı kombinasyonlarının sayısıdır ve
$$ {n \choose k} = \frac{n!}{k!~(n-k)!} $$ile verilir.
Python'da şekil parametreleri olarak n (deneme sayısı) ve p (başarı olasılığı) argümanları olan scipy.stats modülünün binom.rvs() yöntemini kullanarak bir binom dağılımı oluşturulabilir. Dağılımın ortalama değerini değiştirmek loc parametresi kullanılır. Kodun aynı sonuçları verecek şekilde tekrarlanabilmesi için, fonksiyona random_state argümanı eklenebilir.
Örnek: Basketbol oyuncularının faul atış başarıları¶
Diyelim ki profesyonel basketbolcular için ortalama faul başarı yüzdesi %80 olsun. Bu bilgi ışığında hepsinde 10 faul atışının yapıldığı 10000 maçlık rastgele bir Binom dağılımı oluşturalım. Bu dağılımın histogramını sıklık dağılımı yapısında çizdirerek bu 10000 maçtaki 10 faul atışından toplamda n = 0, 1, 2, ..., 10 başarılı atışın gerçekleştiği maç sayılarına bakalım.
from scipy.stats import binom
from matplotlib import pyplot as plt
# scipy stats rastgelelik icin numpy.random
# kullandigindan numpy.random.seed ile seed edebilirsiniz
np.random.seed(seed=233423)
# 10 bagimsiz deney sonucunun p = 0.8 olasilik icin
# 10000 bagimsiz denemede alinma olasiliklari
ornek_binom = binom.rvs(n=10,p=0.8,size=10000)
fig,ax = plt.subplots()
n,bins,patches = ax.hist(ornek_binom, align="right",
range=(0,10), bins=11)
ax.set(xlabel='Binom Dagilimi', ylabel='Frekans')
plt.show()
bins
Görüldüğü üzere de 10 denemenin (örneğin bir basketbol maçında denenen faul atışlarının) hepsinin birden başarıyla sonuçlanma (basketle sonuçlanması) frekansı 10000 ayrı denemede (maçta) ~1000 kez gerçekleşirken, 10 deneminin 8 isabetle sonuçlanması en sık rastlanılan durum olup ~3000 denemede (maçta) bu durumla karşılaşılmıştır. Örnekte basketbolcuların ortalama faul yüzdelerinin %80 olarak belirtilmiş olması nedeniyle bu durum doğaldır. Bir faul atışı ya başarılı ya da başarısız olacağı için bu bir Bernoulli deneyidir. Başarı olasılığı $p = 0.8 > 0.5$ olduğu için dağılım sola doğru çarpıklığa (pozitif) sahiptir (ing. left skewed). $n \rightarrow \infty$, her bir deneme için başarı olasılığı $p \rightarrow 0$ olduğu halde $np = \lambda$ sınırlı bir değerdir. $p$ ve $q = 1-p$'nin 0.5'e yakın olduğu bu durumda dağılım normal dağılıma yakınsar.
Herhangi bir faul atış başarı sayısının kaç maçta gerçekleştiğini belirlemek için dağılımın alındığı numpy diizisindeki bu sayının kaç kez geçtiği belirlenebilir. Bunun için farklı fonksiyonlarla farklı çözümler yapılabilir. İlgili değerler histogramdan doğrudan da alınabilir ($n[7]$ 8.bindeki sayıyı verecektir).
# Herhangi bir basarili atis sayisina kac macta ulasilmis?
np.count_nonzero(ornek_binom == 8)
n[8]
cnts = np.bincount(ornek_binom)
cnts
Bernoulli Dağılımı¶
Bernoulli dağılımı, binom dağılımının tek bir denemenin yapıldığı (n = 1) özel bir hali olup, olasılık kütle fonksiyonu (ing. Probability Mass Function, PMF) şu şekilde verilir:
$$ f(k, p) = p^k q^{(1-k)} = p^k (1-p)^{(1-k)} $$scipy.stats modülünün bernoulli.rvs() fonksiyonu şekil parametresi olarak p (başarı olasılığı) verilmek suretiyle bir Bernoulli Dağılımı oluşturmak için kullanılabilir.
Örneğin hileli (diyelim ki "Yazı" sonucunun her bir atışta 0.6 olasılıkla gerçekleştiği) bir yazı-tura atışının 10000 kez yapılması sonucu oluşan bir rastgele dağılımı inceleyelim.
from scipy.stats import bernoulli
# Tek bir denemenin 10000 kez yapilmasiyla
# p = 0.6 olasilikli sonucun gerceklesme sikligi dagilimi
ornek_bernoulli = bernoulli.rvs(size=10000,p=0.6)
fig,ax = plt.subplots()
n, bins, patches = ax.hist(ornek_bernoulli, bins=2)
ax.set(xlabel='Bernoulli Dagilimi', ylabel='Frekans')
plt.show()
n
Doğal olarak başarı = $1.0$ sonucu $p = 0.6$ olduğu için 10000 denemede 6000'e yakın sayıda gerçekleşmiştir.
Poisson Dağılımı¶
Poisson rasgele değişkeni nadir gerçekleşen bir olayın bir zaman aralığında kaç kez meydana geldiğini modellemek için kullanılır. Örneğin, gökbilimdeki CCD gözlemlerinde kısa poz süreleri dahilinde dedektör üzerinde toplanan foton sayıları ya da bir web sitesinde belirli aralıklarla ziyaret edilen kullanıcı sayısı bir Poisson süreci olarak düşünülebilir.
Bir olayın bir Poisson süreci olarak tanımlanması ve bu olayın verilen bir aralıkta n = 0, 1, 2, 3, ... kez gerçekleşme olasılıklarını belirlemek için Poisson dağılımının kullanılması için 3 gerekli koşul bulunmaktadır.
- Ölçümler nadir görülen olayların sayısıdır.
- Tüm ölçümler birbirinden bağımsızdır.
- Ölçümlerin görülme sıklığı ilgili zaman aralığında değişim göstermemektedir.
Poisson dağılımı olayların meydana gelme oranı ($\lambda$) cinsinden tanımlanır. Bir olay bir aralıkta $n = 0, 1, 2, … $ kez oluşabilir. Verilen bir aralıktaki ortalama olay sayısı $\lambda$ olarak belirtilir. $\lambda$, oran parametresi (ing. rate parameter) olarak da adlandırılan olayın gerçekleşme hızı; bir başka deyişle belirli bir aralıktaki gerçekleşme sayısıdır. Bir olayı verilen belirli bir aralıkta k kez gözleme olasılığı;
$$ P(k) = e^{-\lambda}\frac{\lambda^k}{k!} $$ile verilir.
Aslında Poisson dağılımı, Binom dağılımının, nadir gözlenen bir olayın (örneğin foton sayımı) olasılığının çok küçük olması ($p \lt\lt 1$) nedeniyle ortalama gerçekleşme sayısının ($\lambda$) olası tüm olayların sayısından ($n$) çok daha küçük olduğu ($k = \lambda \lt\lt n$) için özel bir durumudur. Poisson dağılımı ($p \lt\lt 1$) için Binom dağılımının bir yaklaşımı olduğundan ortalama etrafında asimetrik dağılan bir dağılımdır ve şekli 10 kez zar atma deneyinde x kez 1 gelme olasılığının dağılımına benzer.
Poisson dağılımının ortalama değeri aslında tanımında vardır. Ancak ortalama değer, $x$ için beklenen değer ($\bar{x}$) tanımından da çıkarılabilir. Standart sapması da ona eşit olduğundan Poisson dağılımı tek bir parametreyle karakterize edilir ($\mu$).
$\lambda \rightarrow \infty$ limit durumunda Poisson dağılımı Normal dağılıma yakınsar.
scipy.stats paketi fonksiyonlarından poisson.rvs() kullanılarak Poisson dağılımına uyan rastgele bir dağılım oluşturulabilir. Olayın gerçekleştiği aralık dahilinde ortalama olarak olayın kaç kez gözlendiğini veren $\lambda$ parametresi, bir şekil parametresidir ve scipy.stats.poisson foksiyonlarına $mu$ argümanı ile verilr. Bu nedenle Poisson dağılımı tek parametreleri bir dağılımdır.
Örnek: Diyelim ki bir yıldız gözleminde CCD'nin herhangi bir pikseline belirli bir dalgaboyu aralığında yıldızdan saniyede ortalama 5 foton geliyor olsun. Dedektörün bu dalgaboyu aralığındaki kuantum etkinliği de %80 olsun. Yani gelen 5 fotonun 4'ü dedektör tarafından algılansın ($\lambda = 4$). Gece boyunca 1'er saniyelik poz süreleriyle alınan 1000 görüntüde algılanan foton sayılarının değişimi bir Poisson dağılımı sergiler.
from scipy.stats import poisson
ornek_poisson = poisson.rvs(mu=4, size=1000)
fig,ax = plt.subplots()
ax.hist(ornek_poisson)
ax.set(xlabel='Poisson Dagilimi', ylabel='Frekans')
plt.show()
Foton sayıları için gözlenme sayıları (frekans) yerine, normalize yoğunluk (o aralıktaki ölçüm sayısının toplama bölümü) verilmek istenirse pyplot.hist fonksiyonunda density parametresi $True$ değerine ayarlanır. Poisson dağılımının süreksiz bir fonksiyon olduğunu düşünerek olasılık yoğunluk fonksiyonu (PDF) yerine, istenen her değer için Poisson olasılık kütle fonksiyonunu (PMF) örnek dağılımın üzerine çizdirmek bu şekilde mümkündür. Burada PMF, fonksiyon süreksiz olduğunu vurgulamak üzere, özellikle az sayıda (10) nokta için, 'o' sembolüyle verilmiştir.
from scipy.stats import poisson
lamda = 4. # oran parametresi
ornek_poisson = poisson.rvs(mu=lamda, size=1000)
fig,ax = plt.subplots()
ax.hist(ornek_poisson, density = True)
# %1 ile %99 olasiliklar arasina karsilik gelen x degerleri
# arasindaki 10 nokta icin olasilik kutle fonksiyonu (pmf)
x = np.linspace(poisson.ppf(0.01, mu=lamda), poisson.ppf(0.99, mu=lamda),10)
pmf = poisson.pmf(x, mu=lamda)
ax.plot(x, pmf, 'ro')
ax.set(xlabel='Poisson Dagilimi', ylabel='PMF')
plt.show()
Diyelim ki poz süresini uzatıp 5 saniyeye çıkarıyoruz. Ancak birim zamanda algılanan foton sayısı da değişmiyor. Buna karşın artık $\lambda = 5~x~4 = 20$ foton algılarız. Doğal olarak gece boyunca alacağımız görüntü sayısı düşecektir. Poz süresini 5 kat arttırdığımız için toplam görüntü sayısı da 5 kat azalmış ve 200'e düşmüş olsun.
from scipy.stats import poisson
ornek_poisson = poisson.rvs(mu=20, size=200)
fig,ax = plt.subplots()
ax.hist(ornek_poisson)
ax.set(xlabel='Poisson Dagilimi', ylabel='Frekans')
plt.show()
Görüldüğü gibi dağılım normal dağılıma daha fazla yaklaşmıştır. Astronomide gözlemlerinde zaman çözünürlüğünden ve toplam görüntü sayısından ödün verilebilecek durumlarda poz süresinin arttırılmasıyla Poisson dağılımı gösteren foton istatistiğinin Gaussyen dağılımı yakınsaması amaçlanır. Zira görüntü işleme rutinleri, sayımların Normal dağıldığı varsayımı üzerine kuruludur.
Parlak yıldızların gözlemlerinde çok kısa poz süreleri verilebildiğinden Poisson gürültüsü (bu olaya özgü olarak foton gürültüsü adı da verilir) önemli bir problem teşkil eder. Bu nedenle poz sürelerini uzatmak için yıldızın dedektör üzerindeki profilini geniş bir alana yaymak üzere teleskobun odak dışına alındığı odak dışı gözlem tekniğine başvurulabilir. Tekniğin diğer getirileri ile birlikte ötegezegen geçiş gözlemlerine getirdiği hassasiyet için bkz. Baştürk vd. 2015, Baştürk vd. 2020.
Cauchy Dağılımı ve Voigt Profili¶
Astronomide yıldızların atmosfer parametrelerinden elementlerin kimyasal bolluklarına, yüzey parlaklık anomalilerinden dikine hızlarına kadar pek çok özelliklerini belirleyebilmek üzere sıklıkla çalıştığımız tayfsal çizgi profilleri de özünde birer dağılımdır. Dolayısı ile gözlemsel tayflar modellenirken bu dağılımların yapısına uygun dağılım fonksiyonları kullanılır. Parametreleri çizginin dalgaboyu, ve yarı yükseklikteki tam genişliği olan Lorentz, Gauss ve Voigt olasılık dağılım fonksiyonları uygunlukları nedeniyle bu amaçla en sık başvurulan fonksiyonlardır.
Cauchy ya da Lorentz Dağılımı¶
Heisenberg belirsizlik ilkesi gereğince tayfsal çizgi mutlaka bir miktar genişler çünkü enerji düzeylerinin ($\Delta$E) ve parçacıkların bu enerji düzeylerinde bulunma sürelerinin ($\Delta$t) üzerinde her zaman bir miktar belirsizlik bulunur.
$$ \Delta E \times \Delta t \ge \hbar $$Bu genişleme Cauchy ya da fizikte ve astrofizikte daha sık kullanılan adıyla Lorentz profili ile daha iyi temsil edilir. Bu nedenle de "Lorentzian" bir genişleme olarak bilinir.
Cauchy ya da Lorentz dağılımı:
$$ p(x; x_0, \gamma) = \frac{1}{\pi \gamma} (1 + \frac{x - x_0}{\gamma})^2 = \frac{1}{\pi \gamma} \frac{\gamma^2}{(x - x_0)^2 + \gamma^2} $$ile verilir. Burada $x_0$ profil (çizgi) merkezini, $\gamma$ ise yarı-yükseklikteki yarı-genişliği gösterir. Daha sık kullanılan yarı-yükseklikteki tam genişlik (FWHM) bu durumda 2$\gamma$'ya eşittir. Bu parametre doğal olarak çizginin genişliği ile ilişkilendirilir.
Dev yıldızlarda çizgilerin daha dar, cüce yıldızlardaki tayfsal çizgilerin daha geniş olmasına neden olan parçacıkların plazma içinde birbirine yakınlıkları ile artan basınç genişlemesi de Lorentzian bir profil üretir.
Buna mukabil parçacıkların plazma içinde Maxwellian dağılan hareketlerinin neden olduğu genişleme ise Gaussyen (ya da normal dağılıma uyan) bir genişleme yaratır.
Voigt Profili¶
Çoğunlukla bu genişleme mekanizmalarının bileşkesi gözlenir. Lorentz (ya da Cauchy) dağılımı ile Gauss dağılımının konvolüsyonları bir Voigt profili verir ve hassas atmosfer parametrelerine ulaşabilmek için tayfsal modellemede genellikle Voigt profilleri tercih edilir.
$$ V(x; \sigma, \gamma) = G(x; x_0, \sigma) \times L(x; x_0, \gamma) = \int_{-\infty}^{+\infty} {G(x, \sigma) \times L(x - x_0, \gamma) dx} $$Lorentz dağılımının tanımında görüldüğü gibi $\gamma$ yarı-yükseklikteki yarı genişlik, Gauss dağılımında $\sigma$ ise standart sapmadır. Bir Gauss profili için yarı yükseklikteki tam genişlik (FWHM, $\alpha$) ile standart sapma arasındaki ilişki;
$$ \alpha = \sigma \sqrt{2 ln2} $$ile verilir.
Voigt profilinin elde edildiği konvolüsyon işlemi Lorentz ve Gauss fonksiyonlarının yapısı gereği kompleks bir ifade üretir. Bu nedenle Voigt profilinin kapalı bir ifadesi yoktur ancak Faddeva fonksiyonunun Reel kısmı üzerinden ($Re[\omega (z)]$) tanımlanabilir.
$$ V(x; \sigma, \gamma) = \frac{Re[\omega (z)]}{\sigma \sqrt{2 \pi}}, z = \frac{x + i\gamma}{\sigma \sqrt{2}} $$Örnek Voigt Profili¶
Aşağıdaki kod aynı merkezli ($x_0, \mu = x_0 = 5$) yarı yükseklikteki tam genişliği $\gamma = 0.5$ bir Lorentzian profil ile aynı genişlikteki bir Gaussyen profili ($\alpha = 0.5$) karşılaştırmak ve bu iki profilin konvolüsyonundan oluşan Voigt profilini örneklemek için yazılmıştır. Faddeva fonksiyonunun Reel kısmını türetmek için scipy.special.wofz fonksiyonu kullanılmıştır.
import numpy as np
from scipy.special import wofz
def Gauss(x, mu, alpha):
sigma = alpha / np.sqrt(2*np.log(2))
return 1 / (np.sqrt(2*np.pi)*sigma)* np.exp(-0.5*((x-mu) / sigma)**2)
def Lorentz(x, x0, gamma):
return gamma / np.pi / ((x-x0)**2 + gamma**2)
def Voigt(x, x0, alpha, gamma):
sigma = alpha / np.sqrt(2 * np.log(2))
return np.real(wofz(((x-x0) + 1j*gamma)/sigma/np.sqrt(2))) / sigma\
/np.sqrt(2*np.pi)
alpha, gamma = 0.5, 0.5
x0 = 5.
x = np.linspace(-10,10,1000)
plt.plot(x, Gauss(x,x0,alpha), c="b", ls=':', label='Gauss')
plt.plot(x, Lorentz(x,x0,gamma), c="r", ls='--', label='Lorentz')
plt.plot(x, Voigt(x,x0,alpha, gamma), c="m", label='Voigt')
plt.xlim(0,10)
plt.xlabel("x")
plt.ylabel("PDF")
plt.legend(loc="best")
plt.show()
Gamma Dağılımı¶
Gamma dağılımı, iki parametreli, sürekli olasılık dağılımlarının ailesidir. Nadiren ham formunda kullanılırken, üstel (ing. exponential), ki-kare ($\chi^2$) gibi diğer yaygın olaraak kullanılan dağılımlar ve Erlang dağılımları Gamma dağılımının özel durumlarıdır. Gamma dağılımının parametreleri, şekil parametresi, $\alpha = k$ ve oran parametresi olarak adlandırılan, ters ölçek parametresi $\beta = 1 / \theta$ 'dır. Gamma dağılımı için $\Gamma $ sembolü kullanılır ve $\Gamma (n) = (n-1)!$ olarak tanımlanır.
$$ f(x, \alpha, \beta) = \frac{\beta^{\alpha} x^{\alpha - 1} e^{-\beta x}}{\Gamma(\alpha)} $$scipy.stats paketi fonksiyonlarından gamma.rvs() kullanılarak gamma dağılımına uyan rastgele bir dağılım oluşturulabilir. $\alpha$ bir tamsayı olduğunda, gamma dağılımı Erlang dağılımına, $\alpha$ = 1 ise üstel dağılıma karşılık gelir.
numpy.random.gamma fonksiyonu bir gamma dağılımı oluşturmak üzere kullanılabilir. Ölçek parametresinin yukarıdaki tanımlar itibarıyla ters ölçek ($\beta$) değil, ($\beta = 1 / \theta$) tanımındaki $\theta$ parametresi olduğu dikkate alınmalıdır.
Gamma-dağılımının galaksi ve kuazarların ışınım gücünü fonksiyolarını temsil etmek üzere kullanıldığı bir çalışma için bkz. Zaninetti (2019).
np.random.seed(128)
alpha = 5. # sekil parametresi
theta = 0.5 # olcek parametresi
ornek_gamma = np.random.gamma(shape=alpha, scale=theta, size=1000)
fig,ax = plt.subplots()
ax.hist(ornek_gamma, bins=50)
ax.set(xlabel='Gamma Dagilimi', ylabel='Frekans')
plt.show()
Aynı şekilde olasılık yoğunluk fonksiyonu (PDF), birikimli olasılık fonksiyonu (CDF), herhangi bir olasılığın karşılık geldiği aralık gibi gamma sürekli dağılımına ilişkin parametreler scipy.stats.gamma fonksiyonlarıyla elde edilebilir.
alpha = 3.10 # sekil parametresi
# bir gamma dagilimi olusturalim
gamma_ornek = st.gamma.rvs(alpha, size=1000)
# histogramina bakalim
plt.hist(gamma_ornek, density=True, alpha=0.8)
# bir gamma dagiliminin istatistiksel parametrelerine bakalim
ortalama, varyans, carpiklik, kurtosis = st.gamma.stats(alpha, moments='mvsk')
print("Dagilimin ortalama degeri: {:.2f}, standart sapmasi: {:.2f}, carpikligi: {:.2f} ve basikligi: {:.2f}".\
format(ortalama, varyans, carpiklik, kurtosis))
# Olasilik yogunluk fonksiyonu
# %1 ile %99 oalsilik arasini alalim. Tum dagilimlarda
# ppf bu olasiliklara karsilik gelen x degerlerini verir
x = np.linspace(st.gamma.ppf(0.01, alpha), st.gamma.ppf(0.99, alpha), 500)
pdf = st.gamma.pdf(x, alpha)
plt.xlabel("x")
plt.ylabel("PDF")
plt.plot(x, pdf, 'r-', lw=4)
# Bu dagiilma uyan bir ornek olusturalim
plt.show()
Üstel Dağılım¶
Üstel dağılım (ing. exponential distribution), bir Poisson sürecindeki (sürekli ve bağımsız ("hafızasız") bir şekilde sabit bir ortalama oranında meydana gelen) olaylar arasındaki sürenin olasılık dağılım fonksiyonu olarak tanımlanır. Gamma dağılımının özel bir durumu, süreksiz bir olasılık dağılım fonksiyonu olan geometrik dağılımın sürekli durumudur. Normal, binom, Poisson, gamma gibi dağılımların üyesi olduğu üstel dağılımlar ailesinin bir üyesidir. Oran (ing. rate) parametresi olarak adlandırılan bir $\lambda$ parametresine sahiptir ve bu parametre dağılımın azalma (artma) hızını belirler. Olasılık yoğunluk fonksiyonu,
$$ p(x)~dx = \lambda~e^{-\lambda x}~dx $$ile verilir.
Burada $\lambda$ paramatresi Poisson dağılımındaki oran parametresi ($\lambda$)'ya benzer. Örneğin, CCD dedektörün bir pikseline saniyede gelen foton sayısı, Kepler Uzay Teleskobu'nun baktığı alandaki Güneş-benzeri yıldızların sayısı, bir saat içinde bir mağazaya gelen müşteri sayısı, yıllık deprem sayısı gibi nicelikler, Poisson dağılımının parametresi olan oran parametresine karşılık gelir.
Ancak, günlük hayatta olaylar arasındaki geçen süre söz konusu olduğunda, oran yerine zaman birimiyle konuşmyı yeğleriz. Örneğin, 10000 ADU sayım üretecek kadar fotonun algılanması için 20 saniye poz süresinin gerekmesi (20 / 10000 = 0.002 bir oran yerine), bir makinenin hata vermeden çalışacağı yıl sayısının 10 yıl olması (0.1 başarısızlık / yıl demek yerine) gibi. Örnekler çoğaltılabilir, bu nedenle üstel dağılımın “ortalama” terminolojisi Poisson süreçlerdeki oran parametresi $\lambda$ yerine, bunun tersi olan zaman parametresi $1 / \lambda$ üzerine kuruludur. Örneğimizden hareketle 10000 ADU 20 saniyede algılanıyorsa, "Poisson oranı" 0.002'dir, 10000 ADU algılamak için 20 saniye poz süresi vermek gerekir. Dolayısı ile $X \sim exp(0.002)$ verildiğinde kastedilen bu Poisson oranıdır.
Üstel dağılım, bir Poisson olayının gerçekleştiği iki durum arasındaki süre üzerinden tanımlandığından bu iki olay arasındaki bekleme süresinde herhangi bir olay gerçekleşmez. Bu durum aşağıdaki şekilde ifade edilebilir.
$$ P(x = k) = \frac{\lambda^k~e^{-\lambda}}{k!} $$$$ P(x = 0) = \frac{\lambda^0~e^{-\lambda}}{0!} = e^{-\lambda} $$Tek bir zaman biriminde değil de bir $t$ zaman aralığında hiçbir olayın gerçekleşmeme olasılığı bilinmek isteniyorsa P(x=0, 1/t aralıkta) P(x=0, 1/t aralıkta) ... * P(x=0, 1/t aralıkta) çarpımı t kez yapılmalıdır. Bu durumda
$$ P(x = 0, t) = e^{-\lambda~t} $$elde edilir. Bu şekilde elde ettiğimiz olasılık "bir olay gözlendikten sonra bir sonraki olaya kadar olan bekleme süresinin t zaman biriminden büyük olma olasılığı" olup $P(t > T)$ birikimli olasılığıdır. Birikimli olasılık $P(t < T$) ile ifade olunduğundan, bu nicelik $1 - CDF$'e eşittir.
$$ P(t \le T) = 1 - P(t \gt T) = CDF = 1 - e^{-\lambda~t} $$Birikimli yoğunluk fonksiyonunun türevi olasılık yoğunluk fonksiyonunu verecektir.
$$ PDF = \frac{d~CDF}{dt} = \frac{d(1 - e^{-\lambda~t})}{dt} \rightarrow p(t)~dt = \lambda~e^{-\lambda~t}~dt $$olasılık yğunluk fonksiyonu elde edilmiş olur.
Sonuç olarak, birim zaman başına bir Poisson sürecine uyan olay sayısı bir Poisson dağılımını izliyorsa, olaylar arasındaki süre üstel dağılımı takip eder.
lamda = 1.5 # oran (olcek) parametresi
ornek_ustel = np.random.exponential(scale=lamda, size=1000)
fig,ax = plt.subplots()
ax.hist(ornek_ustel, bins=50)
ax.set(xlabel='Ustel Dagilim', ylabel='Frekans')
plt.show()
Aynı şekilde üstel bir dağılım için olasılık yoğunluk fonksiyonu (PDF), birikimli yoğunluk fonksiyonu (CDF), herhangi bir olasılığın karşılık geldiği aralık gibi sürekli dağılıma ilişkin parametreler scipy.stats.expon fonksiyonlarıyla elde edilebilir.
lamda = -1.50 # sekil parametresi
# bir ustel dagilimi olusturalim
ustel_ornek = st.expon.rvs(lamda, size=1000)
# histogramina bakalim
plt.hist(ustel_ornek, density=True, alpha=0.5)
# bir ustel dagilimin istatistiksel parametrelerine bakalim
ortalama, varyans, carpiklik, kurtosis = st.expon.stats(lamda, moments='mvsk')
print("Dagilimin ortalama degeri: {:.2f}, standart sapmasi: {:.2f}, carpikligi: {:.2f} ve basikligi: {:.2f}".\
format(ortalama, varyans, carpiklik, kurtosis))
# Olasilik yogunluk fonksiyonu
# %1 ile %99 oalsilik arasini alalim. Tum dagilimlarda
# ppf bu olasiliklara karsilik gelen x degerlerini verir
x = np.linspace(st.expon.ppf(0.01, lamda), st.expon.ppf(0.99, lamda), 500)
pdf = st.expon.pdf(x, lamda)
plt.xlabel("x")
plt.ylabel("PDF")
plt.plot(x, pdf, 'r-', lw=4)
plt.show()
Weibull Dağılımı¶
log-normal dağılımlar gibi asimetrik dağılımlar ailesinin bir üyesi de Weibull dağılımıdır. Başka parametrizasyonları da bulunmakla birlikte olasılık yoğunluk fonksiyonunun (PDF) standart parametrizasyonu aşağıdaki gibidir.
$$ f(x, \lambda, k) = \frac{k}{\lambda} (\frac{x}{\lambda})^{k-1}~e^{-(x/\lambda)^{k}} $$$k \gt 0$ olmak üzere şekil parametresi, $\lambda \gt 0$ ölçek parametresidir. X bağımsız parametresi herhangi bir niceliğin azaldığı zamanı (başarısızlığa uğrayana kadarki zamanı) gösteriyorsa, Weibull dağılımı onun zamana kuvvet yasasıyla bağlı olan azalma hızını (ing. failure rate: başarısızlık hızı) veren bir dağılımdır. Bu durumda $k-1$ bu kuvvete karşılık gelir.
Örnek olarak $k \lt 1$ durumu, "başarısızlık hızının" azaldığını gösterir. Örneğin bebeklerdeki "ani bebek ölümü" olarak adlandırılan ölüm olasılığı önemli ölçüde yüksekse bu bebeklerin popülasyondan ayrılmasıyla da hızla azalır. Aynı şekilde herhangi bir şekilde defolu ürünlerin erken bozularak popülasyondan ayrılmaları ürünün başarısız bulunma hızını giderek düşürür. $k = 1$ durumu başarısızlık hızının sabit olduğu duruma karşılık gelir (üstel dağılım). Bu durumda örnekteki bebek ölüm oranlarının ve ürünlerin bozulmalarının dışarıdan bir faktörden etkilendiği düşünülür. $k > 1$ ise başarısızlık oranı zamanla artıyor demektir. Örnekte ürünlerin zamanla yaşlandığı ve çalışmaz hale gelebildiği bir senaryo bu duruma uygundur.
Weibull dağılımı, $k = 1$ durumu üstel dağılıma (ing. exponential distribution), $k = 2$, $\lambda = \sqrt{2} \sigma$ ise Rayleigh dağılımına karşılık gelir.
numpy.random.weibull fonksiyonu bir Weibull dağılımına uyan rastgele bir örnek dizisi oluşturmak için kullanılır.
k = 8. # basarisizlik hizinin zamanla arttigi durum
lamda = 1. # olcek parametresi
wb_ornek = np.random.weibull(k, 1000)
def weib(x,l,k):
return (k / l) * (x / l)**(k - 1) * np.exp(-(x / l)**k)
sayi, bins, _ = plt.hist(np.random.weibull(k,1000))
x = np.arange(1,50.)/25.0
olcek = sayi.max()/weib(x, lamda, k).max()
plt.plot(x, weib(x,lamda,k)*olcek)
plt.show()
Olasılık yooğunluk fonksiyonunu örneklem üzerine çizdirmenin iki yolundan biri yoğunluk fonksiyonunu örnekleme ölçeklemek (yukarıdaki örnek), diğer yolu da histogramı density parametresinin değerini $True$ yaparak, olasılık yoğunluğu cinsinden çizdirmektir (aşağıdaki çözüm).
k = 8. # basarisizlik hizinin zamanla arttigi durum
lamda = 1. # olcek parametresi
wb_ornek = np.random.weibull(k, 1000)
def weib(x,l,k):
return (k / l) * (x / l)**(k - 1) * np.exp(-(x / l)**k)
sayi, bins, _ = plt.hist(np.random.weibull(k,1000), density=True)
x = np.arange(1,50.)/25.0
plt.plot(x, weib(x,lamda,k))
plt.show()
Tıpkı diğer dağılım fonksiyonlarında olduğu gibi numpy.random fonksiyonlarıyla dağılıma uyan bir örnek oluşturabilir, olasılık yoğunluk fonksiyonu, herhangi bir değerin herhangi bir aralıka olma olasılığı, birikimli dağılım fonksiyonu gibi özellikler ise scipy.stats.weibull_min fonksiyonu ile elde edilebilir ve incelenebilir; dağılıma uygun bir yapı gösterdiği düşünülen örneklere söz konusu dağılım fonksiyonu uyumlanabilir (fit edilebilir). Şimdi $k < 1$ durumu için bir örneğe bakarak "başarısızlık hızının" azalışını görelim.
k = 0.85 # basarisizlik hizinin zamanla azaldigi durum
l = 1. # olcek parametresi
c = k/l # weibull_min fonksiyonlari icin c = k / lamda parametresi tanimlanir.
ortalama, varyans, carpiklik, kurtosis = st.weibull_min.stats(c, moments='mvsk')
print("Dagilimin ortalama degeri: {:.2f}, standart sapmasi: {:.2f}, carpikligi: {:.2f} ve basikligi: {:.2f}".\
format(ortalama, varyans, carpiklik, kurtosis))
# Olasilik yogunluk fonksiyonu
# %1 ile %99 oalsilik arasini alalim. Tum dagilimlarda
# ppf bu olasiliklara karsilik gelen x degerlerini verir
x = np.linspace(st.weibull_min.ppf(0.01, c), st.weibull_min.ppf(0.99, c), 500)
pdf = st.weibull_min.pdf(x, c)
plt.plot(x, pdf, 'r-', lw=4)
# Bu dagiilma uyan bir ornek olusturalim
wb_ornek = st.weibull_min.rvs(c, size=1000)
plt.hist(wb_ornek, density=True, histtype='stepfilled', alpha=0.6)
plt.show()
# Hangi olasiliklarin hangi x degerlerine
# karsilik geldiklerine bakalim
k = 0.85
xler = st.weibull_min.ppf([0.001, 0.5, 0.999], k)
print("x: ", xler)
Başarısızlık hızının ne kadar hızlı düştüğü hem graikten hem de %1 ve %99 olasılığa karşılık gelen x değerleri arasındaki büyük farkten görülmektedir.
Kuvvet Yasası Dağılımları¶
Kuvvet yasası (zaman zaman "ölçekleme yasası" olarak da adlandırılır), bir nicelikteki göreli bir değişikliğin, ona kuvvet yasasıyla bağlı diğer bir nicelikte orantılı bir değişme yol açtığını belirtir. Yasanın en basit örneği örneklerinden biri ters-kare yasasıdır. İki cisim arasındaki mesafe arttıkça aralarındaki çekim kuvveti mesafenin karesiyle azalır ya da iki özdeş yıldızdan yakın olanın parlaklığı, uzak olanın parlaklığından aralarındaki mesafenin karesi kadar fazladır. Yasaya uyan bir dağılımın olasılık yoğunluk fonksiyonu,
$$ p(x)~dx = a~x^{a - 1} $$ile verilir.
Örneğin $1 M_{\odot}$ ile $10 M_{\odot}$ arasındaki anakol yıldızlarının kütlelerinin olasılık dağılımı $f(M)~\alpha~M^{-3.5}$ ile bağıntılıdır, zira bir yıldızın anakol yaşı kütlesinin -3.5. kuvvetiyle orantılıdır. Yaşamakta olduğumuz CoVid19 hastalığına sebep olan SARS-CoV-2 virüsünün başlangıçtaki yayılma hızının üstel olduğu, sonrasında kuvvet yasasına uyduğu ve yine üstel olarak azaldığı tespit edilmiştir.
numpy.random.power fonksiyonu $[0, 1]$ aralığında verilen $a$ parametresinin 1 eksiğiyle (($a - 1$). kuvvetiyle) kuvvet yasasına uyan bir dağılım oluşturmak için kullanılır. $a \gt 0$ olması gerektiğine, bu anlamda fonksiyonun sadece pozitif kare yasaları (artış) için çalıştığına dikkat edilmelidir. Dağılımı ölçeklendirmek için istenen sayıyla çarpma/bölme yapılabilir.
a = 4 # dorduncu kuvvete
ornek_4kuvvet = np.random.power(a = a, size=1000)*2.5
fig,ax = plt.subplots()
ax.hist(ornek_4kuvvet, bins=50)
ax.set(xlabel='x', ylabel='Frekans')
plt.show()
Diğer dağılımlardakine benzer şekilde olasılık yoğunluk fonksiyonu (PDF), birikimli olasılık fonksiyonu (CDF), herhangi bir olasılığın karşılık geldiği aralık gibi parametreler scipy.stats.powerlaw fonksiyonlarıyla kuvvet yasasına uyan dağılımlar için de elde edilebilir.
a = 3.00 # kuvvet
# bir kuvvet yasasi dagilimi olusturalim
ornek_kuvvet = st.powerlaw.rvs(a, size=1000)
# histogramina bakalim
plt.hist(ornek_kuvvet, density=True, alpha=0.5)
# bir kuvvet yasasi dagiliminin istatistiksel parametrelerine bakalim
ortalama, varyans, carpiklik, kurtosis = st.powerlaw.stats(a = a, moments='mvsk')
print("Dagilimin ortalama degeri: {:.2f}, standart sapmasi: {:.2f}, carpikligi: {:.2f} ve basikligi: {:.2f}".\
format(ortalama, varyans, carpiklik, kurtosis))
# Olasilik yogunluk fonksiyonu
# %0.1 ile %99.9 oalsilik arasini alalim. Tum dagilimlarda
# ppf bu olasiliklara karsilik gelen x degerlerini verir
x = np.linspace(st.powerlaw.ppf(0.001, a), st.powerlaw.ppf(0.999, a), 500)
pdf = st.powerlaw.pdf(x, a)
plt.xlabel("x")
plt.ylabel("PDF")
plt.plot(x, pdf, 'r-', lw=4)
plt.show()
Eksenler logaritmik ölçekte verildiğinde logaritmanın özelliği gereği dağılımın logaritmik ölçekte eğimi kuvvete eşit şekilde doğrusal olduğu görülebilir.
plt.hist(ornek_kuvvet, density=True, alpha=0.5)
plt.plot(x, pdf, 'r-', lw=4)
plt.xscale('log')
plt.yscale('log')
plt.xlabel("x")
plt.ylabel("log(PDF)")
plt.show()
Pareto Dağılımı¶
Pareto dağılımı, kuvvet yasası dağılımları olarak bilinen ve bir $x$ değişkeninin herhangi bir değeri alma olasılığının kuvvet yasasıyla (Pareto için $a = log_45 \sim 1.16$) değiştiği dağılımlar ailesinin bir üyesidir. Adını, İtalyan mühendis ve ekonomist Vilfredo Pareto tarafından öne sürülen ve herhangi bir niceliğe etkilerin %80'inin nedenlerin sadece %20'sinden geldiğini belirten, bilimsel olmayan bir “yasa” olan Pareto İlkesi'nden (ya da Matthew İlkesi) alır. Bu ilkeye göre örneğin bir toplumunun refahının %80'i nüfusun %20'sinin elindedir.
Olasılık Yoğunluk Fonskiyonu'nun (PDF) matematiksel ifadesi,
$$ f_x (x) = \frac{\alpha x^{\alpha}_m}{x^{\alpha} + 1}, x \ge x_m $$$$ f_x (x) = 0, x \lt x_m $$ile verilir.
numpy.random.pareto fonksiyonuya bir örnek dağılım oluşturulup, Pareto fonksiyonu ifadesiyle uyumlanabileceği gibi scipy.stats.pareto fonksiyonları da dağılımın yapısına ve olasılığa ilişkin parametreleri elde etmek üzere kullanılabilir.
import numpy as np
alfa, m = 3., 2. # sekil ve mod parametreleri
s = (np.random.pareto(alfa, 1000) + 1) * m
import matplotlib.pyplot as plt
n, x, _ = plt.hist(s, 100, density=True)
pdf = alfa*m**alfa / x**(alfa+1)
plt.plot(x, max(n)*pdf/max(pdf), linewidth=2, color='r')
plt.xlabel('x')
plt.ylabel('Normalize PDF')
plt.show()
b = 5.80 # sekil parametresi
# bir Pareto dagilimi olusturalim
ornek_pareto = st.pareto.rvs(b, size=1000)
# histogramina bakalim
plt.hist(ornek_pareto, density=True, alpha=0.5)
# bir Pareto dagiliminin istatistiksel parametrelerine bakalim
ortalama, varyans, carpiklik, kurtosis = st.pareto.stats(b = b, moments='mvsk')
print("Dagilimin ortalama degeri: {:.2f}, standart sapmasi: {:.2f}, carpikligi: {:.2f} ve basikligi: {:.2f}".\
format(ortalama, varyans, carpiklik, kurtosis))
# Olasilik yogunluk fonksiyonu
# %1 ile %99 oalsilik arasini alalim. Tum dagilimlarda
# ppf bu olasiliklara karsilik gelen x degerlerini verir
x = np.linspace(st.pareto.ppf(0.01, b), st.pareto.ppf(0.99, a), 500)
pdf = st.pareto.pdf(x, a)
plt.xlabel("x")
plt.ylabel("PDF")
plt.plot(x, pdf, 'r-', lw=4)
plt.show()
# Hangi olasiligin %80'ee karsilik geldigine bakalim
# Pareto ilkesi geregi bu deger x'in
# ilk %20'lik dilimi icerisinde gerceklesir.
b = 5.80
x = st.pareto.ppf(0.20, b)
print("x: ", x)
Pareto ilkesine benzer bir ilke Zipf yasası olarak adlandırılır ve genellikle bir olayın sıralamasına göre sıklığını ifade eder. Zipf yasası, keyfi bir kitapta en sık kullanılan kelimelerin bir listesi verildiğinde, en sık kullanılan kelimenin en sık kullanılan ikinci kelimeden iki kat daha sık göründüğünü belirtir. En basit haliyle, Zipf yasası güç yasasına eşittir. Zipf yasasına uyan dağılımlar kuvvet dağılımları ailesinden Zipf Dağılımı [Zeta Dağılımı](https://en.wikipedia.org/wiki/Zeta_distribution) olarak da adlandırılır. Python'da numpy.random.zipf ve scipy.stats.zipf bu dağılımı çalışmak için gerekli fonskiyonları sunarlar.
Ödev 1¶
Teslim Tarihi: 15 Mart Salı 13:30¶
Çift yıldız ya da ötegezegen sistemlerinde yörüngenin bakış doğrultusuna dik düzlemle yaptığı açı, yörünge eğim açısı (i) olarak adlandırılır ve $0^{\circ}$ ile $90^{\circ}$ arasında değerler alır. Söz konusu parametrenin hangi dağılımla en uygun temsil edilebileceğini düşünerek ve uygun
numpyve / veyascipyfonksiyonlarını kullanarak 10000 sistemden oluşan bir örnek dağılım oluşturunuz. Bu dağılımın sıklık histogramını (y-ekseninde frekans olmalı) uygun sayıda aralık (bin) kullanarak çizdiriniz ve temel istatistiki parametrelerini (ortalama, ortanca, çarpıklık ve basıklık) ekrana yazdırınız.Bu ödev ekindeki odev02_veri_yildizlar_parlakliklar.csv dosyasının birinci sütununda parlaklık değerleri ikinci sütununda ise o parlaklık değerine kadarlık bir görünen parlaklığa sahip yıldızların sayısı virgüllerle ayrılmış şekilde verilmektedir. Dosyayı açarak inceledikten sonra bu parlaklık dağılımının ne tür bir dağılım fonksiyonuyla temsil edilebileceğini belirleyiniz. Uygun
maptlotlib.pyplot,numpyvescipyfonksiyonlarını kullanarak bu dağılımın bir histogramını y-ekseninde Olasılık Yoğunluk Fonksiyonu (PDF) olacak şekilde çizdiriniz. Bu histogramı en iyi temsil edecek dağılımıscipyfonksiyonlarını kullanarak histogramın üzerine çizdiriniz.52 adet draje çikolata içerdiği bilinen tanınmış bir markanın bir paketindeki çikolataların %14'ünün kahverengi renkte olduğu belirlenmiştir (Maddison 2013). Bu bilgiden hareketle 10000 paket çikolatadaki kahverengi çikolataların sayılaırnın dağılımını uygun bir dağılım fonksiyonu kullanarak rastgele oluşturunuz. Oluşturduğunuz bu örnek dağılımı kullanarak a) 6 tane kahverengi draje çikolata içeren paket sayısını, b) En az 10 kahverengi çikolata içeren paket sayısını c) Hiç kahverengi çikolata içermeyen paket sayısını bulunuz d) Örnek dağılımınızda tamamı kahverengi çikolata içeren paket sayısını bularak böyle hiçbir paket olmaması duruumunda aslında bunun mümkün olup olmadığını yorumlayınız.
Bir HII bölgesinin tayfında 656.28 nanometre dalgaboyunda salma (emisyon) çizgisi olarak gözlenen H$_{\alpha}$ çizgisini a) yarı-yükseklikteki yarı-genişliği 10 nm olan bir Lorentz dağılım fonksiyonuyla, b) standart sapması 10 nm olan bir Gauss fonksiyonuyla, c) aynı genişliğe (10 nm) sahip bir Voigt profiliyle ifade ederek oluşturduğunuz dağılım fonksiyonlarını tek bir grafik üzerine çizdiriniz.
Kaynaklar¶
- Python'da Dağılımlar
- Weibull Dağılımı
- Üstel Dağılım
- Üstel Dağılım - 2
- Kuvvet Yasası Dağılımları
- Binom Dağılımı
- scipy.stats Modülü Dokümantasyonu
- numpy.random Modülü Dokümantasyonu
- Ivezic, Z., Connolly, A.J., VanderPlas, J.T., Gray A., 2014, "Statistics, Data Mining, and Machine Learning in Astronomy: A Practical Guide for the Analysis of Survey Data", Princeton University Press