AST416 Astronomide Sayısal Çözümleme - II

Ders - 02b Python'da Belirsizliğin İfadesi ve Yayılımı

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

Belirsizliğin İfadesi

Pek çok kaynağı olabilecek bir ölçüm hatası ya da bir model üzerinden belirlenmiş bir parametrenin tahmini değerinin üzerindeki belirsizlik genel olarak

$$ x \pm \Delta x $$

gösterimiyle ifade edilir.

Buradaki $\Delta x$ değerinin ne olduğu (ölçüm aletinin limiti, standart hata, standart sapma gibi) ayrıca açıkça belirtilmelidir. Hatta bilimsel çalışmalarda parametre belirsizliklerinin nasıl türetildiği için özel bölümler ayrıldığı ve konunun ayrıntısıyla anlatıldığı da olur.

Parametrenin olası değerlerinin olasılık dağılımının elde edildiği ve dağılımın asimetrik olduğu durumlarda negatif ve pozitif yönde belirli bir persantil değeri (genellikle 16 ve 84. persantil tercih edilir) kullanılarak parametre belirsizlikleri verilir.

$$ x^{+\Delta_1}_{-\Delta_2} $$

Belirsizliğin Yayılımı

Bir parametrenin üzerindeki belirsizlik ona bağlı olarak türetilen bir diğer parametreye belirsizliğin yayılmasına ilişkin kurallar takip edilerek yansıtılmalıdır. Bu noktada parametreler arası ilişkiler (korelasyonlar) de dikkate alınır.

Python'da da bu kurallara uyularak belirsizlik hesabı parametre hesabıyla birlikte yapılabileceği, bunun için özel fonksiyonlar yazılabileceği gibi bu amaçla özel olarak geliştirilen uncertainties modülü fonksiyonlarına da başvurulabilir. uncertainties gelişmiş bir paket olup, dokümantasyonu da detaylıdır. Aynı amaçla yazılmış ve veri görselleştirmeleri için de kullanılabilen uncertainty-toolbox gibi başka paketler de bulunmaktadır.

uncertainties Modülü

uncertainties paketi, üzerinde belirsizlik bulunan parametre değerleriyle ($3.14 \pm 0.01$ gibi) hesap yapılmasına ve belirsizliğin bu hesapların yayılmasına olanak sağlayan ücretsiz, platform-bağımsız bir Python modülüdür. Ayrıca herhangi bir ifadenin türevlerini de verebilir.

Bir hesaplamanın karmaşıklığı ne olursa olsun, bu paket, parametre belirsizliklerini hata yayılım teorisinin öngördüğü şekilde uygulayarak hesap sonucunu bir belirsizlik değeriyle döndürür. Türevleri otomatik olarak hesaplar ve belirsizlikleri hesaplamak için kullanır. Hemen hemen tüm belirsizlik hesaplamaları analitik olarak yapılır. Değişkenlerin tanımlanma şekline bağlı olarak aralarındaki korelasyonu da otomatik olarak yönetir.

Birkaç basit örnekle modülün nasıl çalıştığını görelim. Her şeyden önce üzerinde belirsizlik bulunan parametre değerleri ufloat fonksiyonuyla birer uncertainties.core.Variable nesnesi olarak tanımlanmalıdır.

In [1]:
from uncertainties import ufloat
x = ufloat(3,1)
print(x**2 - 2*x + 4)
print(type(x)) 
7+/-4
<class 'uncertainties.core.Variable'>

Parametre hassasiyetini belirtmek için "u" notasyonu (precision modifier) yer tutucu olarak aşağıdaki şekilde kullanılabilir.

In [2]:
x = ufloat(0.2,0.01)
print('1 anlamli rakam : {:.1u}'.format(x))
print('3 anlamli rakam : {:.3u}'.format(x))
print('1 anlamli rakam ussel notasyon : {:.1ue}'.format(x))
print('1 anlamli rakam yuzde notasyonu : {:.1u%}'.format(x))
1 anlamli rakam : 0.20+/-0.01
3 anlamli rakam : 0.2000+/-0.0100
1 anlamli rakam ussel notasyon : (2.0+/-0.1)e-01
1 anlamli rakam yuzde notasyonu : (20+/-1)%

math paketindeki fonksiyonların belirsizliklerle işlem yapabilen ve bu belirsizlikleri hatanın yayılımı prensipleri çerçevesinde sonuçlara yansıtan benzerlerine uncertainties.umath fonksiyonları ile erişilir. Ancak matematiksel sabitlerin math kütüphanesinden indirilmesi gerekir.

In [3]:
from math import pi,e
from uncertainties.umath import *
r = ufloat(2.552,0.005)
print("Cevre = ",(2*pi*r))
Cevre =  16.035+/-0.031
In [4]:
theta = ufloat(pi/4,pi/32)
print(sin(2*theta))
1.000000000000000000+/-0.000000000000000012

Parametreler arası korelasyonlar da otomatik olarak yönetilier. Örneğin $\theta$ parametresi ölçüm belirsizliği ile tanımlanmış olmakla birlikte $\theta - \theta$ işlemine giren parametreler korele olduğu için sonucun üzerindeki belirsizlik 0'dır!

In [5]:
theta - theta
Out[5]:
0.0+/-0

Parametrenin belirsizlik belirleme için tabi tutulduğu türev alma işleminin sonucuna da erişilebilir.

In [6]:
(sin(theta-pi/4)).derivatives[theta]
Out[6]:
1.0

Kovaryans ve Korelasyon Matrisleri

İlgilenilen parametreler arasında bir korelasyonun olması durumunda bu korelasyon; kovaryans ve korelasyon matrisleri incelenerek araştırılabilir. Örneğin bbirbirlerinden bağımsız olarak belirlenen iki sayının toplamı bu sayılara bağlı olacak ve onlarla korelasyon gösterecektir.

In [7]:
u = ufloat(1, 0.1, "u degiskeni")  # Tag (ya da label)
v = ufloat(10, 0.1, "v degiskeni")
toplam = u + 2*v
print(toplam)
21.00+/-0.22
In [8]:
import numpy as np
import uncertainties
kovaryans_matrisi = uncertainties.covariance_matrix([u, v, toplam])
np.set_printoptions(precision=2)
print(np.array(kovaryans_matrisi))
[[0.01 0.   0.01]
 [0.   0.01 0.02]
 [0.01 0.02 0.05]]

Kovaryans matrisinin ilk iki sütunundaki 0 değerleri u ve v değişkenlerinin bağımısız olduğunu , diğer sütunun 0 olmaması ise toplam değişkeninin bu iki değişkene bağlı olduğunu, dolayısıyla aralarında bir korelasyon bulunduğunu göstermektedir.

In [9]:
toplam - (u + 2*v)
Out[9]:
0.0+/-0

İstendiği takdirde kovaryans matrisine erişilebildiği gibi korelasyon matrisine de erişilebilir.

In [10]:
korelasyon_matrisi = uncertainties.correlation_matrix([u, v, toplam])
korelasyon_matrisi
Out[10]:
array([[1.  , 0.  , 0.45],
       [0.  , 1.  , 0.89],
       [0.45, 0.89, 1.  ]])

Aralarında kovaryans matrisinin işaret ettiği ilişkiler bulunan ve değerleri verilen değişkenlerin üzerindeki belirsizlikler correlated_values fonksiyonuyla hesaplanır.

In [11]:
(u2, v2, toplam2) = uncertainties.correlated_values([1, 10, 21], kovaryans_matrisi)
print(u2)
print(v2)
print(toplam2)
1.00+/-0.10
10.00+/-0.10
21.00+/-0.22
In [12]:
toplam2 - (u2 + 2*v2)
Out[12]:
0.0+/-4.712160915387243e-09

Aslında buradaki belirsizliğin de 0 olması gerekirdi ama yuvarlama hataları 0'a yakın ama 0 olmayan bir belirsizliğin elde edilmesine neden olmuştur.

numpy Dizileriyle İşlemlerde Belirsizliğin Yayılımı

uncertainties fonksiyoları numpy dizileriyle birlikte tanımlanarak kullanılabilir.

In [13]:
import numpy as np
x = np.array([ufloat(2.00, 0.01), ufloat(4.00, 0.1)])
x**2
Out[13]:
array([4.0+/-0.04, 16.0+/-0.8], dtype=object)
In [14]:
3*x[0] + x[1]**2 - x[0]*x[1]
Out[14]:
14.0+/-0.6000833275471

Nesnenin "__print()__" işlem sonundaki anlamlı rakam sayısını dikkate alarak sonuç üzerindeki belirsizliğin ifadesine yardımcı olur.

In [15]:
print(3*x[0] + x[1]**2 - x[0]*x[1])
14.0+/-0.6

Ancak istenirse formatlama seçenekleri kullanılarak da sonucun hassasiyetinin ifadesi yönteilebilir. Yuvarlama işlemlerinde Parçacık Fiziği Grubu'nun takip ettiği kurallar takip edilmektedir.

In [16]:
print("{:.4f}".format(3*x[0] + x[1]**2 - x[0]*x[1]))
14.0000+/-0.6001

unumpy Alt Modülü

Her ne kadar ufloat üzerinde uygulanan fonksiyonlar bir numpy dizisinin elemanlarına uygulanabilir olsa da unumpy alt modülü daha ileri düzey ihtiyaçlar için kolaylıklar getirir.

Bir dizi tanımlarken önce parametre değerleri (nominal values) sonra da belirsizlikleri (uncertainties) liste, tuple ya da numpy dizisi gibi bir "iterable" ile sağlanır.

In [17]:
from uncertainties import unumpy
dizi = unumpy.uarray([0, 1], [0.1, 0.005])
dizi
Out[17]:
array([0.0+/-0.1, 1.0+/-0.005], dtype=object)

Bir dizinin içindeki parametre değerleri nominal_values, belirsizlikleri std_devs fonksiyonlarıyla alınabilir.

In [18]:
print(unumpy.nominal_values(dizi))
print(unumpy.std_devs(dizi))
[0. 1.]
[0.1  0.01]

Tıpkı ufloat için olduğu gibi unumpy için de matematik fonksiyonları tanımlıdır.

In [19]:
theta = unumpy.uarray([0,pi/3,pi/4,pi/6,pi/2],[pi/16]*5)
theta
Out[19]:
array([0.0+/-0.19634954084936207,
       1.0471975511965976+/-0.19634954084936207,
       0.7853981633974483+/-0.19634954084936207,
       0.5235987755982988+/-0.19634954084936207,
       1.5707963267948966+/-0.19634954084936207], dtype=object)
In [20]:
print(unumpy.cos(theta))
[1.0+/-0 0.5000000000000001+/-0.1700436903969579
 0.7071067811865476+/-0.13884009181744894
 0.8660254037844387+/-0.09817477042468102
 6.123233995736766e-17+/-0.19634954084936207]

Bu değerleri belirsizlikleri ile birlikte bir dosyaya yazdırmak istediğinizde savetxt fonksiyonunu kullanabilirsiniz.

In [21]:
np.savetxt('cos_theta.txt', unumpy.cos(theta), fmt='%r')

İstendiğinde bu dosya tekrar uncertainties.ufloat_fromstr() kullanılarak okunabilir.

In [22]:
converters = dict.fromkeys(range(1),
                          lambda col_bytes: uncertainties.ufloat_fromstr(col_bytes.decode("utf8")))
cos_th = np.loadtxt('cos_theta.txt', converters=converters, dtype=object)
In [23]:
cos_th
Out[23]:
array([1.0+/-0, 0.5000000000000001+/-0.1700436903969579,
       0.7071067811865476+/-0.13884009181744894,
       0.8660254037844387+/-0.09817477042468102,
       6.123233995736766e-17+/-0.19634954084936207], dtype=object)

umatrix Alt Modülü

Matris işlemleri yapabilmek için tanımlanmış ayrıca bir umatrix alt modülü ve bu modülün fonksiyonları da bulunmaktadır.

In [24]:
mat1 = unumpy.umatrix([1, 2], [0.01, 0.002])
mat2 = unumpy.umatrix([[3],[4]], [0.01, 0.002])
mat1*mat2
Out[24]:
matrix([[11.0+/-0.038262252939417984, 11.0+/-0.03136877428271625]],
       dtype=object)
In [25]:
print(mat1.I)
[[0.19999999999999996+/-0.0012419339757008014]
 [0.3999999999999999+/-0.0016178998732925337]]
In [3]:
import numpy as np
a = np.random.randn(100).reshape(50,2)
print(a)
[[-9.58768977e-01 -2.41936919e+00]
 [ 1.83186250e+00 -1.75614060e-03]
 [-7.39012095e-02 -5.00170979e-01]
 [ 2.62597391e-01  6.96365839e-01]
 [ 1.19807291e+00 -6.81759574e-01]
 [ 6.01467900e-02  9.43985250e-01]
 [-3.77594071e-02  3.94386050e-01]
 [-2.01150054e+00 -2.35749031e-01]
 [ 2.72122257e-01  3.19376704e-01]
 [-5.31655028e-01 -7.35166739e-01]
 [-1.09476475e+00  1.66553531e-01]
 [ 9.39529041e-01  8.47774223e-01]
 [-4.33642222e-02 -4.82351146e-01]
 [-7.66941510e-02  6.54933229e-02]
 [-5.64700881e-01  6.33179991e-01]
 [-3.16627946e-01  1.32846556e+00]
 [-2.11907854e-01 -5.79099109e-01]
 [-4.10227012e-01 -5.59856441e-01]
 [ 4.37026100e-02  1.48536927e-01]
 [-1.06659966e-01 -1.77405548e+00]
 [ 1.08470920e+00  2.42604270e-01]
 [-1.72547731e-01 -5.32229771e-01]
 [ 5.05916208e-01  1.31299466e+00]
 [ 7.17799069e-01  1.23804938e+00]
 [ 1.10465676e+00 -4.98777554e-01]
 [ 2.96792855e+00  1.71986958e+00]
 [-5.31537531e-01 -2.51902773e-01]
 [-1.03751896e+00  9.24722507e-01]
 [ 3.33422805e-01 -1.24050393e+00]
 [ 2.12620929e-01 -1.86312605e-01]
 [ 2.97233029e-01 -8.88261717e-01]
 [-5.99065948e-01  6.48334126e-01]
 [-1.39862310e+00  3.63465565e-01]
 [ 4.96115263e-01  1.52173604e+00]
 [-1.74437629e+00  1.30511292e-01]
 [ 2.71153135e-01  8.88845794e-01]
 [-8.78973609e-01  2.32657197e-01]
 [-2.56681589e+00 -6.93787792e-01]
 [-3.53612842e-01  1.16486674e+00]
 [ 4.11244541e-01 -2.18301379e-01]
 [ 1.89505239e+00  2.86626652e+00]
 [ 1.53322491e+00 -5.97504145e-01]
 [ 8.82892824e-01  1.23005412e-02]
 [-1.67616531e+00  3.20129923e-01]
 [-1.13859973e+00  5.65385674e-01]
 [ 1.19040850e+00  2.86942100e-01]
 [ 2.89481806e-01 -1.22548062e+00]
 [ 4.35432209e-01 -6.29258084e-01]
 [-8.36729687e-01  4.33021680e-01]
 [-3.74726643e-01  7.48202193e-02]]
In [11]:
np.savetxt("a.txt",a)
b = np.loadtxt("a.txt")
a.T[0]
Out[11]:
array([-0.95876898,  1.8318625 , -0.07390121,  0.26259739,  1.19807291,
        0.06014679, -0.03775941, -2.01150054,  0.27212226, -0.53165503,
       -1.09476475,  0.93952904, -0.04336422, -0.07669415, -0.56470088,
       -0.31662795, -0.21190785, -0.41022701,  0.04370261, -0.10665997,
        1.0847092 , -0.17254773,  0.50591621,  0.71779907,  1.10465676,
        2.96792855, -0.53153753, -1.03751896,  0.3334228 ,  0.21262093,
        0.29723303, -0.59906595, -1.3986231 ,  0.49611526, -1.74437629,
        0.27115314, -0.87897361, -2.56681589, -0.35361284,  0.41124454,
        1.89505239,  1.53322491,  0.88289282, -1.67616531, -1.13859973,
        1.1904085 ,  0.28948181,  0.43543221, -0.83672969, -0.37472664])