{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# AST416 Astronomide Sayısal Çözümleme - II #\n", "## Ders - 04 Eğri Uyumlama ##" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Doç. Dr. Özgür Baştürk
\n", "Ankara Üniversitesi, Astronomi ve Uzay Bilimleri Bölümü
\n", "obasturk at ankara.edu.tr
\n", "http://ozgur.astrotux.org" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Bu derste neler öğreneceksiniz?#\n", "## Eğri Uyumlama ##\n", "\n", "* [Eğri Uyumlama](#Eğri-Uyumlama)\n", "* [Lineer Regresyon](#Lineer-Regresyon)\n", " * [Örnek: Lineer Regresyon](#Örnek:-Lineer-Regresyon)\n", " * [Regresyon Katsayısı ve Tahmin Üzerindeki Standart Hata](#Regresyon-Katsayısı-ve-Tahmin-Üzerindeki-Standart-Hata)\n", "* [Lineer Olmayan İlişkilerin Lineerleştirilmesi](#Lineer-Olmayan-İlişkilerin-Lineerleştirilmesi)\n", " * [Üstel Fonksiyonların Lineerleştirilmesi](#Üstel-Fonksiyonların-Lineerleştirilmesi)\n", " * [Kuvvet Fonksiyonlarının Lineerleştirilmesi](#Kuvvet-Fonksiyonlarının-Lineerleştirilmesi)\n", " * [Rasyonel Bazı Fonksiyonların Lineerleştirilmesi](#Rasyonel-Bazı-Fonksiyonların-Lineerleştirilmesi)\n", "* [Örnek: Lineer Olmayan İlişkilerin Lineerleştirilmesi](#Örnek:-Lineer-Olmayan-İlişkilerin-Lineerleştirilmesi)\n", "* [Polinom Regresyonu](#Polinom-Regresyonu)\n", " * [Örnek: Polinom Regresyonu](#Örnek:-Polinom-Regresyonu)\n", " * [numpy.polyfit İle Parabol Uyumlaması Örneği ve Uyumlamanın Başarısı](#numpy.polyfit-İle-Parabol-Uyumlaması-Örneği-ve-Uyumlamanın-Başarısı)\n", "* [Çok Değişkenli Lineer Regresyon](#Çok-Değişkenli-Lineer-Regresyon)\n", " * [Örnek: Çok Değişkenli Lineer Regresyon](#Örnek:-Çok-Değişkenli-Lineer-Regresyon)\n", " * [Çok Değişkenli Lineer Regresyon: Genel İfade](#Çok-Değişkenli-Lineer-Regresyon:-Genel-İfade)\n", "* [En Küçük Kareler Minimizasyonu Yöntemi Genel İfade](#En-Küçük-Kareler-Minimizasyonu-Yöntemi-Genel-İfade)\n", "* [Doğrusal Olmayan En Küçük Kareler Yöntemi](#Doğrusal-Olmayan-En-Küçük-Kareler-Yöntemi)\n", " * [Levenberg Marquardt Algoritması](Levenberg-Marquardt-Algoritması)\n", " * [scipy.optimize.leastsq ile Uyumlama](#scipy.optimize.leastsq-ile-Uyumlama)\n", " * [scipy.optimize.curve_fit ile Uyumlama](#scipy.optimize.curve_fit-ile-Uyumlama)\n", "* [LMFIT Paketiyle Modelleme](#LMFIT-Paketiyle-Modelleme)\n", "* [Eğri Uyumlamanın İnterpolasyon ve Ekstrapolasyon Amacıyla Kullanımı](#Eğri-Uyumlamanın-İnterpolasyon-ve-Ekstrapolasyon-Amacıyla-Kullanımı)\n", "* [Kaynaklar](#Kaynaklar)\n", "* [Ödev 4](#Ödev-4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Eğri Uyumlama #\n", "\n", "Sadece astronomide değil, deneysel veriye bağlı herhangi bir bilim dalında gözlem veriye bir matematiksel modelin uyumlanması en sık başvurulan sayısal çözümleme tekniklerinden biridir. Uyumlanan matematiksel model, çoğu zaman analitik ifadesi bilinen bir eğridir. Bu eğriden yola çıkılarak i) verinin zaman içinde nasıl değiştiğini anlamak (trend analizi), ii) interpolasyon ve ekstrapolasyon yapmak; iii) herhangi bir model önerisinin başarısını ve öngörülen matematiksel modelin karşılık geldiği hipotezi sınamak olanaklı hale gelir. \n", "\n", "# Lineer Regresyon #\n", "\n", "Bir veriye yapılabilecek en basit model bir doğru uyumlamaktır. Eğri uyumlamada izlenen klasik istatistiksel (frekansçı; ing. frequentist) yaklaşımları anlamak için temel olarak lineer regresyonu ele alalım ve konuya en basit eğri (doğru) uyumlama yöntemiyle başlayalım. Lineer regresyonda eldeki veriye\n", "\n", "$ y = a_0 + a_1 x $\n", "\n", "modeli uyumlanmaya çalışılır. Model, deneysel veriyle ölçüm olayının kaçınılmaz bir şekilde hata barındırması nedeniyle birebir aynı olmayacağı için onun da üzerinde bir hata bulunacaktır. Lineer regresyonda amaç bu hatayı minimize etmektir. Hata $y_{deney}$ deneysel (ya da gözlemsel) ölçümleri, $y_{model}$ lineer modeli göstermek üzere, aşağıdaki şekilde ifade edilebilir. \n", "\n", "$ e = y_{deney} - y_{model} = y_{deney} - a_0 + a_1 x $\n", "\n", "Hatayı minimize etmek için çeşitli stratejiler düşünülebilir.\n", "\n", "(a) Doğrudan bu farkın minimize edilmesi hedeflenebilir. Bu durumda her bir deney sonucu ile model arasındaki farklar alınır ve bu farkı minimum yapan model (doğru denklemi) aranır. $n$, nokta sayısını göstermek üzere;\n", "\n", "$$ \\sum_{i=1}^{n} e_i = \\sum_{i=1}^{n} (y_{i} - (a_0 + a_1 x_i)) $$\n", "\n", "Bu durumda birden fazla farklı model aynı sonucu verebilir. Aşağıdaki grafikte gösterildiği gibi doğru model (sürekli doğru) ile açıkça yanlış olduğu görülen ikinci bir model (süreksiz doğru) modelden farklar toplandığında aynı $\\Sigma e_i$ değerini verir. Birinci modelde gözlemsel veriler doğru üzerinde olduğu için modelden farkların her biri 0 olduğundan farklar toplamı da 0'dır. Aynı şey ikinci modelde her bir noktanın modelden farkı 0 olmasa da geçerlidir; zira ilk noktanın modelden farkı negatif olup, son nokta ile onun pozitif olması dışında büyüklüğü aynıdır. Ortadaki nokta da modelin üzerinde olduğu için modelden farklar toplamı yine 0'dır.\n", "\n", "
\n", "\"dogru\n", "
\n", "\n", "(b) Artıkların pozitif ya da negatif olabilmesi sebebiyle birbirlerini götürmelerinden kaçınmak için mutlak değer alınabilir. Bu stratejide minimize edilmesi istenen nicelik aşağıdaki şekilde ifade edilebilir.\n", "\n", "$$ \\sum_{i=1}^{n} |e_i| = \\sum_{i=1}^{n} |y_{i} - (a_0 + a_1 x_i)| $$\n", "\n", "Bu durumda da birden fazla model aynı sonucu verebilir. Aşağıdaki dört gözlemsel nokta iki ayrı doğruyla $\\Sigma |e_i| $ değerini verecek şekilde temsil edilebilir.\n", "\n", "
\n", "\"dogru\n", "
\n", "\n", "(c) Her iki problemden de kaçınmak için fark karelerin toplamının minimize edilmesi yoluna gidilir ki. Bu gayet iyi bilinen en küçük kareler minimizasyonu (ing. least squares minimization) tekniğidir. Bu teknikte aşağıda $S_r$ ile verilen niceliği minimize edecek model (burada doğru) aranır.\n", "\n", "$$ S_r = \\sum_{i=1}^{n} e_i^2 = \\sum_{i=1}^{n} (y_{i} - (a_0 + a_1 x_i))^2 $$\n", "\n", "Bu klasik bir minimizasyon problemidir. Çok iyi bildiğiniz gibi bir fonksiyonun uç değerleri (örneğin minimumu) birinci türevinin 0 olduğu noktada aranır. Ancak burada fonksiyonun türevinin alınacağı değişkenler $(x,y)$ değişkenleri değildir. Zira bu değerler değişmeyen, gözlemsel veriyi $(x_i,y_i)$ oluşturmaktadırlar. Değştirerek bu fonksiyonu minimize edecek değerlerin arandığı parametreler ise lineer modelin parametreleridir $(a_0,a_1)$. Dolayısıyla fonksiyonun türevinin alınacağı değişkenler de bu değişkenler olmalıdır. Daha sonra bu türevler 0'a eşitlenerek $a_0$ ve $a_1$'in bu farkları minimize edecek değerleri aranacak ve bu değerlerin işaret ettiği doğru \"en iyi doğru modeli\" olacaktır. Bu model tektir (ing. unique) ve bu stratejide ulaşılan en iyi modeldir.\n", "\n", "$$ \\frac{\\partial S_r}{\\partial a_0} = -2 \\sum_{i=1}^{n} (y_{i} - a_0 - a_1 x_i) $$\n", "$$ \\frac{\\partial S_r}{\\partial a_1} = -2 \\sum_{i=1}^{n} [(y_{i} - a_0 - a_1 x_i) x_i] $$\n", "\n", "Bu türevler 0'a eşitlenip, toplam da bileşenlerine ayrılacak olursa;\n", "\n", "$$ \\sum_{i=1}^{n} y_{i} - \\sum_{i=1}^{n} a_0- \\sum_{i=1}^{n} x_i a_1 = 0 $$\n", "$$ \\sum_{i=1}^{n} x_{i} y_{i} - \\sum_{i=1}^{n} x_i a_0 - \\sum_{i=1}^{n} x_i^2 a_1 = 0 $$\n", "\n", "İfadeler aşağıdaki şekilde düzenlenebliir. $n$ tane $a_0$ toplanınca $n~a_0$ elde edileceğine ve $a_0$ ile $a_1$'in toplam ifadelerinin dışına alınabileceğine dikkat ediniz.\n", "\n", "$$ (n)~a_0 + (\\sum_{i=1}^{n} x_i) a_1 = \\sum_{i=1}^{n} y_{i} $$\n", "$$ (\\sum_{i=1}^{n} x_i) a_0 + (\\sum_{i=1}^{n} x_i^2) a_1 = \\sum_{i=1}^{n} x_{i} y_{i}$$\n", "\n", "Bu basit iki bilinmeyenli ($a_0$ ve $a_1$) birinci dereceden iki denklem içeren bir denklem takımıdır, çözümü de tektir ve aşağıdaki şekilde ifade edilebilir (gösteriniz!).\n", "\n", "$$ a_1 = \\frac{n~\\Sigma x_{i} y_{i} - \\Sigma x_i \\Sigma y_i}{n~\\Sigma x_{i}^2 - (\\Sigma x_{i})^2} $$\n", "\n", "$$ a_0 = \\bar{y} - a_1~ \\bar{x}, \\bar{x} = \\frac{\\Sigma x_i}{n}, \\bar{y} = \\frac{\\Sigma y_i}{n} $$\n", "\n", "Parametrelerin hataları;\n", "\n", "$$ \\epsilon_{a_1} = \\sqrt{\\frac{(n~\\Sigma y_i^2 - (\\Sigma y_i)^2 - a_1^2~(n \\Sigma x_i^2 - (\\Sigma x_i)^2))}{(n-2)(n~\\Sigma x_i^2 - (\\Sigma x_i)^2)}}$$\n", "\n", "$$ \\epsilon_{a_0} = \\sqrt{\\epsilon_{a_1}^2 \\frac{\\Sigma x_i^2}{n}} $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Örnek: Lineer Regresyon ##\n", "\n", "Bir roketin hızına karşılık maruz kaldığı hava direnci aşağıdaki çizelgede verilmiştir. Bu veriye en küçük kareler yöntemiyle en uygun doğru uyumlamasını yapınız.\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
i $x_i$ $y_i$
010 25
12070
230380
340550
450610
5601220
670830
7801450
\n", "\n", "Bu problemi öncelikle $a_0$ ve $a_1$ parametrelerini kendi fonksiyonlarımızla hesaplayacak şekilde çözelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import numpy as np\n", "def dogru_denklemi(x,y):\n", " \"\"\"\n", " x ve y dizileri ile verilen veriye en iyi uyan dogrunun\n", " y-eksenindeki kayma (a0) ve egim (a1) parametrelerini\n", " donduren fonksiyon\n", " \"\"\"\n", " n = len(x)\n", " a1 = (n*sum(x*y) - (sum(x)*sum(y))) / (n*sum(x**2) - (sum(x)**2))\n", " a0 = (sum(y)/n) - a1*(sum(x)/n)\n", " return(a0,a1)\n", "\n", "x = np.arange(10,81,10.)\n", "y = np.array([25.,70.,380.,550., 610.,1220.,830.,1450.])\n", "n,m = dogru_denklemi(x,y)\n", "print(\"Verilen veriye uyumlanan en iyi dogru: y = {:.4f}x + {:.4f}\".format(m,n))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Şimdi veriyi ve en küçük kareler yöntemiyle veriye uyumladığımız en iyi doğru modelinin bir grafiğini çizdirelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "# oncelikle x'in uc degerleri arasinda\n", "# 50 noktadan olusan yeni bir x ekseni tanimlayalim\n", "xyeni = np.linspace(x[0],x[-1])\n", "# Bu ekseni dogru denkleminde yerine koyalim\n", "yyeni = m*xyeni + n\n", "plt.plot(x, y, 'o', xyeni, yyeni, '-')\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(['veri', 'model'], loc = 'best')\n", "plt.plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Şimdi aynı soruyu bu kez numpy.polfyit() fonksiyonu yardımıyla Python'a çözdürelim. [numpy.polyfit](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html) istenen dereceden bir polinomu $x$ ve $y$ parametreleriyle sağlanan veriye uyumlayan bir fonksiyondur. 1. dereceden bir polinomun doğru olduğu bilgisinden hareketle $deg$ parametresini 1'e eşitleyerek aranan çözüme ulaşılabilir. Fonksiyon uyumlanan doğrunun (polinomun) katsayılarını ve istendiğinde parametre hatalarının türetilebileceği kovaryans matrisini verir." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "import numpy as np\n", "%matplotlib inline\n", "katsayilar = np.polyfit(x, y, deg = 1) # Sondaki 1 dogru uyumlamasıini gosterir\n", "print(katsayilar)\n", "print(\"uyumlanan dogrunun katsayilari: \", katsayilar)\n", "# Istenirdse bu katsayilar tek tek kullanilarak dogrunun \n", "# y degerleri hesaplabilecegi gibi,\n", "# bu amaacla poly1d fonksiyonu da kullanilabilir.\n", "dogru = np.poly1d(katsayilar)\n", "ydogru = dogru(xyeni)\n", "plt.plot(x, y, 'o', xyeni, ydogru, '-')\n", "plt.xlabel(\"x\")\n", "plt.xlabel(\"y\")\n", "plt.legend(['veri', 'model'], loc = 'best')\n", "plt.plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Regresyon Katsayısı ve Tahmin Üzerindeki Standart Hata ##\n", "\n", "Bu basit örnek sonrası uyumlanan doğrunun veri setini ne derecede temsil edebildiğine yönelik parametrelere geçelim. Regresyon katsayısı, bir doğrusal modelin eldeki veri setini aritmetik ortalamaya oranla ne ölçüde temsil edebildiğinin bir ölçüsüdür. Bu parametre aşağıdaki şekilde hesaplannır ve 1'e ne kadar yakınsa doğru modeli eldeki veri setine ortalamadan o denli daha iyi bir yaklaşımken, 0'a ne kadar yakınsa aritmetik ortalamaya göre temsil niteliği de o derece azdır.\n", "\n", "$S_r$ verinnin modelden, $S_t$ aritmetik ortalamadan fark kareleri toplamını ifade etmek üzere aşağıdaki şekilde tanımlanır:\n", "\n", "$$ S_r = \\sum_{i=1}^{n} (y_{i} - (a_0 + a_1 x_i))^2 $$\n", "$$ S_t = \\sum_{i=1}^{n} (y_{i} - \\bar{y_i})^2, \\bar{y_i} = \\frac{\\Sigma yi}{n} $$\n", "\n", "Korelasyon katsayısı $r$ ise,\n", "\n", "$$ r^2 = \\frac{|S_t - S_r|}{S_t} $$\n", "\n", "şeklinde tanımlandığı için açıkça modelin ortalamaya göre başarısını ifade etmektedir. Aşağıdaki grafiklerde aynı veri setinin aritmetik ortalamayla (solda) ve doğru modeliyle (sağda) hangi başarımda temsil edildiği model verinin modell etrafındaki dağılımıyla gösterilmiştir.\n", "\n", "
\n", "\"lineer\n", "
\n", "\n", "Nümerik bir örnek için tekrar roket hızına karşılık hava direnci örneğimize dönelim ve bu örnek için kendi yazdığımız dogru_uyumlama fonksiyonu korelasyon katsayısını da hazırlamak üzere biraz geliştirelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def lineer_regresyon(x,y):\n", " \"\"\"\n", " x ve y dizileri ile verilen veriye en iyi uyan dogrunun\n", " y-eksenindeki kayma (a0) ve egim (a1) parametrelerini\n", " ve lineer korelasyon katsayisini hesaplayarak\n", " donduren fonksiyon\n", " \"\"\"\n", " n = len(x)\n", " a1 = (n*sum(x*y) - (sum(x)*sum(y))) / (n*sum(x**2) - (sum(x)**2))\n", " a0 = (sum(y)/n) - a1*(sum(x)/n)\n", " Sr = sum((y - a0 - a1*x)**2)\n", " St = sum((y - y.mean())**2)\n", " r2 = (np.abs(St - Sr) / St)\n", " return(a0,a1,r2)\n", "\n", "x = np.arange(10,81,10.)\n", "y = np.array([25.,70.,380.,550., 610.,1220.,830.,1450.])\n", "n,m,r2 = lineer_regresyon(x,y)\n", "print(\"Verilen veriye uyumlanan en iyi dogru: y = {:.4f}x + {:.4f}\".format(m,n))\n", "print(\"Korelasyon katsayisi r = {:.2f}\".format(np.sqrt(r2)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Korelasyon katsayısı r'nin yerine genellikle karesi $r^2 = 0.94^2 = 0.88$ kullanılır. Zira $r^2$, örneğimizde lineer modelin, verinin üzerindeki belirsizliğin %88'ini açıklayabildiğini, ortalamaya oranla getirdiği iyileştirmenin bu düzeyde olduğunu göstermektedir.\n", "\n", "Bu doğruya dayanarak yapılacak tahminlerin üzerindeki standart hata, tahmin üzerindeki standart hata parametresi ($S_{y/x}$) (ing. standard estimation error) ile verilir. \n", "\n", "$$ S_{y/x} = \\sqrt{\\frac{S_r}{n - 2}} $$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Sr = sum((y - n - m*x)**2)\n", "Syx = np.sqrt(Sr / (len(y) - 2))\n", "print(\"Tahmin üzerindeki standart hata: {:g}\".format(Syx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ortalamanın standart sapması ($S_y$) ise \n", "\n", "$$ S_{y} = \\sqrt{\\frac{S_t}{n - 1}} $$\n", "\n", "ile verilir. Standart sapma, sadece ortalama üzerinden hesaplandığı için serbestlik derecesi (ing. degrees of freedom) $n - 1$ iken, tahmin üzerindeki standart hata ise $a_0$ ve $a_1$ şeklinde iki parametre içeren bir model üzerinden hesaplandığı için serbestlik derecesi, $n - 2$'dir. Dolayısı ile tahmin üzerindeki standart hata aslında lineer modelin standart sapmasıdır. Bu durumu verinin ortalama ve lineer model etrafında nasıl dağıldığına bakarak görebiliriz. Tıpkı aşağıdaki şekilde olduğu gibi örneğimizde de verinin lineer model etrafındaki dağılımının düzeyini ifade eden tahmin üzerindeki standart hatanın (sağdaki dağılımın standart sapmasının), ortalama etrafındaki dağılımını ifade eden standart sapmadan küçük olmasını ($S_{y/x} < S_y$) bekleriz; zira korelasyon katsayısının ($r^2 = 0.88$) 1'e yakın olması bize bunu söylemektedir.\n", "\n", "
\n", "\"lineer\n", "
\n", "\n", "Ortalama etrafındaki dağılımın standart sapmasını hesaplayacak olursak:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Sy = np.sqrt(np.sum((y-y.mean())**2) / (len(y) - 2))\n", "print(\"Standart sapma: {:g}\".format(Sy))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lineer modelin, veriyi temsil etmekteki ortalamaya göre başarısı buradan da görülebliir: $S_{y/x} = 189.79 < S_y = 548.98$.\n", "\n", "Regresyon katsayısı aslında tek başına verinin lineer modelle ortalamay göre ne denli başarıyla temsil edildiğini vermektedir. Dolayısyla verinin kendi içindeki değişimin doğasından bağımsız olarak değerlendirilmemesi gerekir. Aşağıda verilen dört grafikte verilen dört farklı veri seti için $y = 3 + 0.5~x$ doğrusunun lineer korelasyon katsayıları ($r^2$) eşittir ([Anscombe 1973](http://www.sjsu.edu/faculty/gerstman/StatPrimer/anscombe1973.pdf)). Bu grafikten hareketle lineer korelasyon katsayısına bakmadan önce verinin değişiminin doğasına bakılması gerektiği görülmektedir. \n", "\n", "
\n", "\"Anscombe\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lineer Olmayan İlişkilerin Lineerleştirilmesi #\n", "\n", "Bazı doğrusal olmayan ilişkileri logaritmanın ve üs alma işleminin özelliklerini kullanarak doğrusal hale getirilmesi (lineerleştirilmesi) mümkündür. Bu özellikler, bu ilişkilerin analitik formlarının belirlenmesinde büyük kolaylık sağlar, zira lineer model basit bir modeldir ve uyumlama işlemi de görece kolaydır.\n", "\n", "## Üstel Fonksiyonların Lineerleştirilmesi ##\n", "\n", "Üstel fonksiyonlar logaritma alınarak aşağıdaki şekilde doğrusallaştırılabilir.\n", "\n", "$$ y = \\alpha_1 e^{\\beta_1 x} $$\n", "\n", "olmak üzere, her iki tarafın logaritması alınacak olursa;\n", "\n", "$$ ln(y) = ln(\\alpha_1) + \\beta_1~x $$\n", "\n", "şeklinde ifade edilebilir. Bu durumda ilişkiyi lineerleştirmek için yapılması gereken tek şey her iki tarafın logaritmasını almak olduğundan; veri setindeki $x$ dizisi olduğu gibi bırakılırken $y$ ölçümlerinin logartiması alınır ve yukarıda anlatıldığı şekilde en küçük kareler yöntemiyle (ya da daha sofistike bir teknikle) lineer regresyon yapılarak doğrunun parametreleri ($a_0 = ln(\\alpha_1)$ ve $a_1 = \\beta_1$) hesaplanır. Bu parametreler hesaplandıktan sonra istenirse ilişki doğrusal form ya da ($ln(y) = ln(\\alpha_1) + \\beta_1~x$) üstel formda ($ y = \\alpha_1 e^{\\beta_1 x} $) ifade edilebilir.\n", "\n", "## Kuvvet Fonksiyonlarının Lineerleştirilmesi ##\n", "\n", "Kuvvet fonksiyonları da logaritma alınarak aşağıdaki şekilde doğrusallaştırılabilir.\n", "\n", "$$ y = \\alpha_2 x^{\\beta_2} $$\n", "\n", "olmak üzere, her iki tarafın logaritması alınacak olursa;\n", "\n", "$$ log(y) = log(\\alpha_2) + \\beta_2~log(x) $$\n", "\n", "şeklinde ifade edilebilir. Bu durumda ilişkiyi lineerleştirmek için yapılması gereken tek şey her iki tarafın logaritmasını almak olduğundan; veri setindeki $x$ ve $y$ dizilerinin logartimaları alınır ve yukarıda anlatıldığı şekilde en küçük kareler yöntemiyle (ya da daha sofistike bir teknikle) lineer regresyon yapılarak doğrunun parametreleri ($a_0 = log(\\alpha_2)$ ve $a_1 = \\beta_2$) hesaplanır. Bu parametreler hesaplandıktan sonra istenirse ilişki doğrusal form ($log(y) = log(\\alpha_2) + \\beta_2~log(x)$) veya kuvvet yasası formunda ($ y = \\alpha_2 x^{\\beta_2} $) ifade edilebilir.\n", "\n", "## Rasyonel Bazı Fonksiyonların Lineerleştirilmesi ##\n", "\n", "Aşağıdaki yapıda verilen rasyonel bir fonksiyon eşitliğin her iki tarafının 1'e bölünmesiyle doğrusallaştırılabilir.\n", "\n", "$$ y = \\alpha_3 \\frac{x}{\\beta_3 + x} $$\n", "\n", "olmak üzere, her iki taraf 1 ile bölünecek olursa;\n", "\n", "$$ \\frac{1}{y} = \\frac{1}{\\alpha_3} + \\frac{\\beta_3}{\\alpha_3}~\\frac{1}{x} $$\n", "\n", "şeklinde ifade edilebilir. Bu durumda ilişkiyi lineerleştirmek için yapılması gereken tek şey her iki tarafı 1 ile bölmek olduğundan; veri setindeki $x$ ve $y$ dizileri de 1 ile bölünür ($1/x$ ve $1/y$) ve yukarıda anlatıldığı şekilde en küçük kareler yöntemiyle (ya da daha sofistike bir teknikle) lineer regresyon yapılarak doğrunun parametreleri ($a_0 = 1 / \\alpha_3$ ve $a_1 = \\beta_3 / \\alpha_3$) hesaplanır. Bu parametreler hesaplandıktan sonra istenirse ilişki istenen formda ifade edilebilir.\n", "\n", "Aşağıda her üç durum da grafiklerle gösterilmiştir.\n", "\n", "
\n", "\"Anscombe\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Örnek: Lineer Olmayan İlişkilerin Lineerleştirilmesi #\n", "\n", "Bir roketin hızına karşılık maruz kaldığı hava direnci örneği üzerinden devam ederek, bu örnekteki veriye daha iyi bir model bulunup bulunamayacağına bakalım. Hava direncinin roket hızının kuvvetiyle değiştiğini varsayarak, kuvvet yasasıyla bir ilişki aramayı deneyelim. Bu durumda,\n", "\n", "$$ y = \\alpha x^{\\beta} $$\n", "\n", "olmak üzere, her iki tarafın logaritması alınacak olursa;\n", "\n", "$$ log(y) = log(\\alpha) + \\beta log(x) $$\n", "\n", "doğrusal ilişkisi için veri ile en küçük fark kareleri verecek modelin $\\alpha$ ve $\\beta$ parametrelerini arayalım. Bun amaçla daha önce yazdığımız lineer_regresyon fonksiyonunu kullanabiiriz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "x = np.arange(10,81,10.)\n", "y = np.array([25.,70.,380.,550., 610.,1220.,830.,1450.])\n", "logx = np.log10(x)\n", "logy = np.log10(y)\n", "n,m,r2 = lineer_regresyon(logx,logy)\n", "print(\"Verilen veriye uyumlanan en iyi dogru: y = {:.4f}x + {:.4f}\".format(m,n))\n", "print(\"Korelasyon katsayisi r = {:.2f}\".format(np.sqrt(r2)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Görüldüğü gibi bu modelin korelasyon katsayısı daha da yüksek ($r^2 = 0.94$) çıkmıştır. Bu durumda bu modeli diğer modele göre daha üstün kabul edebiliriz. Her iki modeli de verini üzerine çizdirelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "xyeni = np.linspace(x[0],x[-1])\n", "ykuvvet = (10**n)*(xyeni**m)\n", "plt.plot(x, y, 'o', xyeni, ydogru, '--', xyeni, ykuvvet, '-')\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(['veri', 'lineer', 'kuvvet'], loc = 'best')\n", "plt.plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bu yeni model her iki ekseni de logaritmik olarak ayarlanmak suretiyle doğru şeklinde de çizdirilebilir." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(x, y, 'o', xyeni, ykuvvet, '-')\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(['veri', 'kuvvet'], loc = 'best')\n", "plt.plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Polinom Regresyonu #\n", "\n", "Eğri uyumlamadaki temel yaklaşımı en küçük kareler minimizasyonu yöntemiyle gördükten sonra aynı yaklaşımı polinom uyumlama için de kullanabiliriz. Doğada iki değişken arasındaki pek çok ilişki lineer ya da lineerleştirilebilir değildir. Bir polinomla temsil edilebilecek veri setleri için de uyumlanabilecek en iyi polinomun derecesini belirldikten sonra katsayılarını belirlemek üzere en küçük kareler minimizasyonu tercih edilebilir. Örneğin ikinci dereceden bir polinom (parabol) uyumlanmak istendiğinde model;\n", "\n", "$$ y = a_0 + a_1~x + a_2~x^2 $$\n", "\n", "şeklinde ifade edilir. Yine minimize edilmek istenen parametre ($S_r$) verinin modelden fark karelerinin toplamıdır ve modelin hatasını (başarımının derecesini) sayısallaştırır. \n", "\n", "$$ S_r = \\sum_{i=1}^{n} (y_{i} - (a_0 + a_1 x_i + a_2 x_i^2))^2 $$\n", "\n", "Yine tipik bir minimizasyon problemi olan bu problemin çözümünde minimize edilmek istenen fonksiyonun ($S_r$), değişkenlerine ($a_0$, $a_1$, $a_2$) göre birinci türevi alınıp sıfıra eşitlenerek minimum değeri aranır. Fonksiyon ifadesindeki $x_i$, $y_i$ değerlerinin değişmez ölçümler olduğunu hatırlatmakta fayda vardır.\n", "\n", "$$ \\frac{\\partial S_r}{\\partial a_0} = -2 \\sum_{i=1}^{n} (y_{i} - a_0 - a_1 x_i - a_2 x_i^2) $$\n", "$$ \\frac{\\partial S_r}{\\partial a_1} = -2 \\sum_{i=1}^{n} [(y_{i} - a_0 - a_1 x_i - a_2 x_i^2) x_i] $$\n", "$$ \\frac{\\partial S_r}{\\partial a_2} = -2 \\sum_{i=1}^{n} [(y_{i} - a_0 - a_1 x_i - a_2 x_i^2) x_i^2] $$\n", "\n", "ifadeleri 0'a eşitlenip düzenlendiğinde;\n", "\n", "$$ (n)~a_0 + (\\sum_{i=1}^{n} x_i) a_1 + (\\sum_{i=1}^{n} x_i^2) a_2 = \\sum_{i=1}^{n} y_{i} $$\n", "$$ (\\sum_{i=1}^{n}x_i) a_0 + (\\sum_{i=1}^{n} x_i^2) a_1 + (\\sum_{i=1}^{n} x_i^3) a_2 = \\sum_{i=1}^{n} x_{i} y_{i}$$\n", "$$ (\\sum_{i=1}^{n}x_i^2) a_0 + (\\sum_{i=1}^{n} x_i^2) a_1 + (\\sum_{i=1}^{n} x_i^3) a_2 = \\sum_{i=1}^{n} x_{i} y_{i}$$\n", "\n", "Çözülmesi gereken denklem sistemi budur. Bu denklem sistemini matris formunda ifade edecek olursak;\n", "\n", "\\begin{equation*}\n", "\\begin{bmatrix}\n", "n & \\Sigma x_i & \\Sigma x_i^2 \\\\\n", "\\Sigma x_i & \\Sigma x_i^2 & \\Sigma x_i^3 \\\\\n", "\\Sigma x_i^2 & \\Sigma x_i^3 & \\Sigma x_i^4\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "a_0 \\\\\n", "a_1 \\\\\n", "a_2\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "\\Sigma y_i \\\\\n", "\\Sigma x_i y_i \\\\\n", "\\Sigma x_i^2 y_i\n", "\\end{bmatrix}\n", "\\end{equation*}\n", "\n", "Bu denklem $ R x A = M $ şeklinde bir lineer denklemdir ve çözümü oldukça kolaydır. Yapılması gereken, bilinmeyen katsayıları içeren $A = [a_0, a_1, a_2]$ matrisini yalnız bırakmak üzere denklemin her iki tarafını $R$ matrisinin tersi $R^{-1}$ ile çarpmaktır. \n", "\n", "$$ R x A = M \\Rightarrow R^{-1} x R x A = R^{-1} x M \\Rightarrow I x A = R^{-1} x M \\Rightarrow A = R^{-1} x M $$\n", "\n", "## Genelleştirilmiş Çözüm: ##\n", "\n", "Bu ifade genelleştirilerek herhangi bir dereceden polinom uyumlaması için aşağıdaki çözüm elde edilebilir.\n", "\n", "$$ y = a_0 + a_1~x + a_2~x^2 + ... + a_m~x^m $$\n", "$$ S_r = \\sum_{i=1}^{n} (y_{i} - (a_0 + a_1 x_i + a_2 x_i^2 + ... + a_m x^m))^2 $$\n", "\n", "olmak üzere;\n", "\n", "\\begin{equation*}\n", "\\begin{bmatrix}\n", "n & \\Sigma x_i & \\Sigma x_i^2 & \\cdots & \\Sigma x_i^m \\\\\n", "\\Sigma x_i & \\Sigma x_i^2 & \\Sigma x_i^3 & \\cdots & \\Sigma x_i^{m+1} \\\\\n", "\\Sigma x_i^2 & \\Sigma x_i^3 & \\Sigma x_i^4 & \\cdots & \\Sigma x_i^{m+2} \\\\\n", "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "\\Sigma x_i^m & \\Sigma x_i^{m+1} & \\Sigma x_i^{m+2} & \\cdots & \\Sigma x_i^{2m} \\\\\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "a_0 \\\\\n", "a_1 \\\\\n", "a_2 \\\\\n", "\\vdots \\\\\n", "a_m\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "\\Sigma y_i \\\\\n", "\\Sigma x_i y_i \\\\\n", "\\Sigma x_i^2 y_i \\\\\n", "\\vdots \\\\\n", "\\Sigma x_i^m y_i\n", "\\end{bmatrix}\n", "\\end{equation*}\n", "\n", "Daha sonra tıpkı kuadratik polinom (parabol) çözümünde olduğu gibi matris formunda ifade edilen bu lineer denklem aşağıdaki şekilde çözülür ve en iyi polinomun katsayıları elde edilir.\n", "\n", "$$ R x A = M \\Rightarrow R^{-1} x R x A = R^{-1} x M \\Rightarrow I x A = R^{-1} x M \\Rightarrow A = R^{-1} x M $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Örnek: Polinom Regresyonu ##\n", "\n", "Aşağıdaki tabloda verilen veri setine uyan en iyi polinomu en küçük kareler minimzasyonu yöntemiyle belirleyerek veri seti üzerine bu polinomu çizdiriniz.\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
$x_i$ $y_i$
02.1
17.7
213.6
327.2
440.9
561.1
\n", "\n", "Öncelikle problemi yukarıdaki eşitliklerden yararlanarak kendi fonksiyonlarımızı yazmak suretiyle kendimiz çözelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def parabol_uyumlama(x,y):\n", " n = len(y)\n", " Sxi = np.sum(x)\n", " Sxi2 = np.sum(x**2)\n", " Sxi3 = np.sum(x**3)\n", " Sxi4 = np.sum(x**4)\n", " Syi = np.sum(y)\n", " Sxiyi = np.sum(x*y)\n", " Sxi2yi = np.sum(x**2*y)\n", " R = np.matrix([[n, Sxi, Sxi2], [Sxi, Sxi2, Sxi3], [Sxi2, Sxi3, Sxi4]])\n", " M = np.matrix([[Syi],[Sxiyi],[Sxi2yi]])\n", " Rinv = np.linalg.inv(R)\n", " A = np.dot(Rinv,M)\n", " return(A)\n", "x = np.linspace(0,5,6)\n", "y = np.array([2.1,7.7,13.6,27.2,40.9,61.1])\n", "katsayilar = parabol_uyumlama(x,y)\n", "print(katsayilar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Görüldüğü gibi denklem takımımızı matris formuna getirmek üzere [`numpy.matrix`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.matrix.html) fonksiyonundan faydalandık ve her bir satırı birer liste olan listeleri hazırlayarak matris formuna dönüştürmüş olduk. $\\Sigma x_i$, $\\Sigma x_i^2$ gibi toplamları birden çok kez kullanmamız gerektiği için değişkenlere atadık ve bu değişkenleri `numpy.matrix` fonksiyonundaki listelere matris formundaki lineer denklemimizde verildiği şekliyle aldık. Daha sonra $R$ matrisinin tersini ($R^{-1}$) `numpy.linalg.inv` fonksiyonuyla aldık ve sonrasında $R^{-1}~M$ matris çarpımını `numpy.dot` fonksiyonuyla yaptık. Burada iki önemli noktayı hatırlatmakta fayda var; birincisi iki matrisin doğrudan çarpımında ($R * M$) matrislerdeki karşılıklı elemanlar çarpıldığı için bunun bir matris çarpımı olmadığı, bunu ancak `numpy.dot` fonksiyonuyla yapabildiğimiz; ikincisi de burada sadece pedagojik olması açısından matris işlemlerini tek tek kendimizin yaptığı. Oysa ki `numpy.linalg.solve()` fonksiyonunu kullanarak bunu doğrudan da yapabilirdik. [`numpy.linalg`](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html) modülü lineer cebir uygulamaları için kullanılabilecek pek çok fonksiyon sağlaması bakımından oldukça kullanışlıdır. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def parabol_uyumlama(x,y):\n", " n = len(y)\n", " Sxi = np.sum(x)\n", " Sxi2 = np.sum(x**2)\n", " Sxi3 = np.sum(x**3)\n", " Sxi4 = np.sum(x**4)\n", " Syi = np.sum(y)\n", " Sxiyi = np.sum(x*y)\n", " Sxi2yi = np.sum(x**2*y)\n", " R = np.matrix([[n, Sxi, Sxi2], [Sxi, Sxi2, Sxi3], [Sxi2, Sxi3, Sxi4]])\n", " M = np.matrix([[Syi],[Sxiyi],[Sxi2yi]])\n", " A = np.linalg.solve(R,M)\n", " return(A)\n", "x = np.linspace(0,5,6)\n", "y = np.array([2.1,7.7,13.6,27.2,40.9,61.1])\n", "katsayilar = parabol_uyumlama(x,y)\n", "print(katsayilar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bu kolaylık fonksiyonumuzu sadece bir satır kısalttı ama bir miktar gereksiz iş yükünden de kurtardı. Bu tür durumlara programlarımız ve betiklerimizde dikkat etmeli, Python ya da `numpy`, `scipy` gibi güvenilir paketlerdeki uygun fonksiyonları ile yapabildiklerimizi öğrenip mümkün olduğunca ilgili fonksiyonları kullanmalıyız. Bu yaklaşım hem kodlarımızı hızlandıracak, hem de yazma, okuma ve dokümantasyon kolaylığı sağlayacaktır. Özünde $numpy$ ve $scipy$ fonksiyonlarını kullanarak parabol de uyumlayabiliriz ama şimdilik öğretici olması açısından katsayıları kendimiz hesaplamyı tercih ettik. Şimdi uyumladığımız parabolü ve verimizi çizdirelim ve bir sonraki aşamada bu örneği bir de `numpy.polyfit` fonksiyonuyla çözelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "xyeni = np.linspace(0,5,50)\n", "a0, a1, a2 = katsayilar\n", "yparabol = float(a0) + float(a1)*xyeni + float(a2)*xyeni**2\n", "plt.plot(x,y,'ro',label=\"veri\")\n", "plt.plot(xyeni,yparabol,'b-',label=\"model\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Korelasyon katsayısı ($r^2$) ve tahmin üzerindeki standart hata ($S_{y/x}$) terimlerini hesaplamak için de bir fonksiyon yazalım ve uyumlamanın başarısını da sayısallaştıralım." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def parabol_uyumlama_basarim(x,y,katsayilar):\n", " n = len(y)\n", " a0, a1, a2 = katsayilar\n", " St = np.sum((y - y.mean())**2)\n", " Sr = np.sum((y - (float(a0)+float(a1)*x+float(a2)*x**2))**2)\n", " r2 = np.abs((St - Sr) / St)\n", " Syx = np.sqrt(Sr / (n-2))\n", " return(r2,Syx)\n", "r2,Syx = parabol_uyumlama_basarim(x,y,katsayilar)\n", "print(\"Verilen veriye uyumlanan en iyi parabol: y = {:.4f}x^2 + {:.4f}x + {:.4f}\".\\\n", " format(float(a2),float(a1),float(a0)))\n", "print(\"Korelasyon katsayisi r^2 = {:.4f}\".format(r2))\n", "print(\"Tahmin üzerindeki standart hata: {:.4f}\".format(Syx))\n", "Sy = np.sqrt(np.sum((y - y.mean())**2) / (len(y)-2))\n", "print(\"Ortalama modelin standart sapması: {:.4f}\".format(Sy))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bu fonksiyonu parabol_uyumlama fonksiyonuyla birleştirmek ve uyumlamanın parametreleri (parabolün katsayıları) ile birlikte korelasyon katsayısı ve tahmin üzerindeki standart hatayı da bu fonksiyondan döndürmekte fayda vardır (deneyiniz!). Grafiği gördükten sonra anlamlı olmasa dahi, parabolik modelin ortalamaya göre ne oranda başarılı olduğunu göstermek üzere ortalama modelin standart sapması ile tahmin üzerindeki standart hatayı karşılaştıabileceğiniz gibi $r^2 = 0.0085$ değerinin 1'e yakınlığına da bakabilirsiniz." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## numpy.polyfit İle Parabol Uyumlaması Örneği ve Uyumlamanın Başarısı ##\n", "\n", "Aynı problemi daha önce de gördüğümüz `numpy.polyfit` fonksiyonundan yararlanarak çözelim. Bu kez fonksiyonu kullanırken $full$ parametresini $True$ değerine ayarlayarak çalıştıralım. Bu durumda sadece uyumlanan polinomun katsayılarını değil, aynı zamanda modelden artıkları ($S_r$), katsayılar matrisinin rank'ini, singüler olduğu değerleri, singüler değerlerin dikkate alınmayacağı limit değerini de verir. Bunlardan son üçü modelin hatasını hesaplarken işe yaramaz. Bu nedenle bu üç parametreyi \"_\" 'ye ayarlayarak kullanmayacağımızı belirtmiş olduk. Python'da birden fazla değer döndüren fonksiyonlarda kullanılmayacak değerler genellikle \"_\"'e atanarak gözardı edilir. Ayrıca istenirse y değerlerinin ölçüm hataları da fonksiyona $w$ parametresi ile geçirilerek her bir noktanın hatası ile ($1 / e_i^2$) ağırlıklandırılmak suretiyle uyumlama yapılır. Bu durumda hatası büyük olan noktaya daha az ağırlık veirlmiş olur. Bu özellikten faydalanmak üzere veriye bir de hata sütunu ekleyelim.\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
$x_i$ $y_i$ $e_i$
02.10.5
17.71.0
213.63.5
327.20.5
440.92.5
561.10.5
\n", "\n", "Ayrıca uyumlamanın başarısını göstermek için indirgenmiş $\\chi^2$ parametresini de hesaplayalım ve çıktı olarak sağlayalım. İndirgenmiş $\\chi^2$ (\"ki-kare\" şeklinde okunur) parametresi aşağıdaki şekilde tanımlanır.\n", "\n", "$$ \\chi^2_\\nu = \\frac{\\chi^2}{\\nu} $$\n", "\n", "Burada sistemin $\\nu$ serbestlik derecesi (dof, ing. degree of freedom), veri noktası sayısı ile model paramatreleri sayısı arasındaki farktır ($\\nu = n - m$). $\\chi^2$ parametresi ise gözlenen noktalar ile model arasındaki fark karelerinin varyansa (ya da noktanın hatasının karesine) oranlarının toplamıdır ve aşağıdaki şekilde verilir.\n", "\n", "$$ \\chi^2 = \\sum_{i=1}^{n} \\frac{(y - y_i)^2}{e_i^2} $$\n", "\n", "`numpy.polyfit` tarafından döndürülen artıklar, her bir noktanın hatasının $w$ parametresinden sağlanması durumunda aslında doğrudan $\\chi^2$'ye eşittir. Zira verinin modelden fark kareler toplamı, noktanıın hatasının karesine bölünerek toplanması suretiyle hesaplanır. İndirgenmiş $\\chi^2$ değerini bulmak için yapılması gereken sadece bu değeri serbestlik derecesine bölmektir. \n", "\n", "$ \\chi^2 >> 1 $ olduğu durumda uyumlama oldukça başarısız ve veriyi temsil etmekten uzaktır. Genel olarak $ \\chi^2_\\nu > 1 $ olması, fitin veri üzerindeki belirsizliği tam olarak modelleyemediği anlamına gelir. Tüm noktalar uyumlanan eğri üzerinde değilse durum zaten bu olmalıdır. İndirgenmiş $ \\chi^2 = 1 $ ise gözlemler ile modelin, varyansla (ya da hatalarla) uyumlu bir şekilde örtüştüğü anlamına gelir. Bu durumda $\\chi^2_\\nu$'nin 1'e yakın olması, uyumlamanın (ya da modelin) da başarılı olduğu şeklinde yorumlanır. $\\chi^2_\\nu < 1$ durumu ise ancak modelden artıkların serbestlik derecesinden küçük olmasıyla mümkündür. Bu durum, verinin göstermediği kadar yüksek dereceden bir fonksiyonla uyumlandığı (ing. overfitting) ya da üzerindeki varyansın olduğundan fazla tahmin edildiği (ing. overestimated) şeklinde yorumlanır. \n", "\n", "Örnek problemimizde nokta sayısı $n = 6$, modelin parametre sayısı $m = 3$ olduğundan (polinomun derecesi 2 iken parametre sayısı m = 3'tür), serbestlik derece $\\nu = n - m = 6 - 3 = 3$'tür. Elde edilen indirgenmiş $\\chi^2_\\nu$ değeri 1'e yakın olup, parabolik modelin bu veri seti için kabul edilebilir, iyi bir model olduğunu göstermektedir.\n", "\n", "Not: Model başarımı ve karşılaştırması için kullanılan diğer parametre ve testleri model karşılaştırması konusunda ayrıntılı olarak göreceğiz. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import math\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "x = np.linspace(0,5,6)\n", "y = np.array([2.1,7.7,13.6,27.2,40.9,61.1])\n", "yerr = np.array([0.5,1.0,3.5,0.5,2.5,0.5])\n", "katsayi,Sr,_,_,_ = np.polyfit(x,y,deg=2,w=yerr,full=True) #sondaki 2, 2. dereceden polinomu gosterir\n", "parabol = np.poly1d(katsayi)\n", "xyeni = np.linspace(0,5) #varsayilan deger 50 noktadir\n", "yparabolnumpy = parabol(xyeni)\n", "plt.errorbar(x,y,yerr=yerr,c='r',fmt='.',label=\"veri\")\n", "plt.plot(xyeni,yparabolnumpy,'b-',label=\"model\")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(loc=\"best\")\n", "plt.show()\n", "Sr = Sr[0] # dondurulen artiklar bir dizidir\n", "St = np.sum((y - y.mean())**2)\n", "r2 = np.abs((St - Sr) / St)\n", "Syx = np.sqrt(Sr / (len(y) - 2))\n", "nu = len(y) - len(katsayi) #serbestlik derecesi\n", "indChi2 = Sr / nu\n", "print(\"Korelasyon katsayisi r^2 = {:.4f}\".format(r2))\n", "print(\"Tahmin üzerindeki standart hata: {:.4f}\".format(Syx))\n", "print(\"İndirgenmiş ki-kare: {:.4f}\".format(indChi2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Çok Değişkenli Lineer Regresyon #\n", "\n", "Astronomide ilgilenilen pek çok problemin tek bir değişkeni yoktur. Dolayısı ile modellemede birden fazla değişkene bağlılığın uyumlanması gerekir. Birden fazla değişkenin doğrusal bir modelinin elde edilmesi özünde tek değişkenli lineer regresyondan farklı değildir. Örneğin üç boyutta düşünecek olursak $x_1$ ve $x_2$ bağımsız değişkenleri göstermek üzere lineer bir model\n", "\n", "$$ y = a_0 + a_1~x_1 + a_2~x_2 $$\n", "\n", "şeklinde ifade edilebilir. Yine minimize edilmek istenen parametre ($S_r$) verinin modelden fark karelerinin toplamıdır ve modelin hatasını (başarımının derecesini) sayısallaştırır.\n", "\n", "$$ S_r = \\sum_{i=1}^{n} (y_{i} - (a_0 + a_1 x_{1,i} + a_2 x_{2,i}))^2 $$\n", "\n", "Yine tipik bir minimizasyon problemi olan bu problemin çözümünde minimize edilmek istenen fonksiyonun ($S_r$), değişkenlerine ($a_0$, $a_1$, $a_2$) göre birinci türevi alınıp sıfıra eşitlenerek minimum değeri aranır. Fonksiyon ifadesindeki $x_{1,i}$, $x_{2,i}$ ve $y_i$ değerleri yine değişmez ölçümler olduğundan çözüm parabol uyumlamaya çok benzerdir.\n", "\n", "$$ \\frac{\\partial S_r}{\\partial a_0} = -2 \\sum_{i=1}^{n} (y_{i} - a_0 - a_1 x_{1,i} - a_2 x_{2,i}) $$\n", "$$ \\frac{\\partial S_r}{\\partial a_1} = -2 \\sum_{i=1}^{n} [(y_{i} - a_0 - a_1 x_{1,i} - a_2 x_{2,i}) x_{1,i}] $$\n", "$$ \\frac{\\partial S_r}{\\partial a_2} = -2 \\sum_{i=1}^{n} [(y_{i} - a_0 - a_1 x_{1,i} - a_2 x_{2,i}) x_{2,i}] $$\n", "\n", "ifadeleri 0'a eşitlenip düzenlendiğinde;\n", "\n", "\\begin{equation*}\n", "\\begin{bmatrix}\n", "n & \\Sigma x_{1,i} & \\Sigma x_{2,i} \\\\\n", "\\Sigma x_{1,i} & \\Sigma x_{1,i}^2 & \\Sigma x_{1,i} x_{2,i} \\\\\n", "\\Sigma x_{2,i} & \\Sigma x_{1,i} x_{2,i} & \\Sigma x_{2,i}^2\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "a_0 \\\\\n", "a_1 \\\\\n", "a_2\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "\\Sigma y_i \\\\\n", "\\Sigma x_{1,i} y_i \\\\\n", "\\Sigma x_{2,i} y_i\n", "\\end{bmatrix}\n", "\\end{equation*}\n", "\n", "Denklemin çözümü, birden fazla değişken olduğu için bir eğri (lineer modelde doğru) olmayıp, bir düzlem verecektir (aşağıdaki şekil). Şekilde görülen beyaz noktalar uyumlanan en iyi düzlemin altında, siyah noktalar ise üstünde kalan veri noktalarıdır. Dik doğrular modele (uyumlanan en iyi düzleme) uzaklıkları göstermektedir.\n", "\n", "
\n", "\"Cok\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Örnek: Çok Değişkenli Lineer Regresyon #\n", "\n", "Diyelim ki aşağıdaki tabloda verilen veriye lineer bir model arıyor olalım.\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
$x_{1,i}$ $x_{2,i}$ $y_i$
005
2110
2.529
130
463
7227
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def duzlem_uyumlama(x1,x2,y):\n", " n = len(y)\n", " Sx1i = np.sum(x1)\n", " Sx2i = np.sum(x2)\n", " Sx1ix2i = np.sum(x1*x2)\n", " Sx1i2 = np.sum(x1**2)\n", " Sx2i2 = np.sum(x2**2)\n", " Syi = np.sum(y)\n", " Sx1iyi = np.sum(x1*y)\n", " Sx2iyi = np.sum(x2*y)\n", " R = np.matrix([[n, Sx1i, Sx2i], [Sx1i, Sx1i2, Sx1ix2i], [Sx2i, Sx1ix2i, Sx2i2]])\n", " M = np.matrix([[Syi],[Sx1iyi],[Sx2iyi]])\n", " A = np.linalg.solve(R,M)\n", " return(A)\n", "x_1 = np.array([0.0,2.0,2.5,1.0,4.0,7.0])\n", "x_2 = np.array([0.0,1.0,2.0,3.0,6.0,2.0])\n", "y = np.array([5.0,10.0,9.0,0.0,3.0,27.0])\n", "katsayilar = duzlem_uyumlama(x_1,x_2,y)\n", "print(katsayilar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Şimdi uyumladığımız düzlemi ve verimizi 3-boyutlu bir grafik üzerinde görelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "%matplotlib inline\n", "fig = plt.figure(figsize=(12,8))\n", "ax = fig.gca(projection='3d')\n", "# Veri noktalarini cizdirme\n", "ax.scatter3D(x_1, x_2, y, color=\"r\", marker=\"o\", s=50)\n", "# Duzlem cizdirme\n", "X1, X2 = np.meshgrid(x_1, x_2)\n", "a0,a1,a2 = katsayilar\n", "Y = float(a0) + float(a1)*X1 + float(a2)*X2\n", "ax.plot_surface(X1, X2, Y, color=\"gray\", alpha=0.5, rstride=1)\n", "ax.set_xlabel(\"$x_1$\")\n", "ax.set_ylabel(\"$x_2$\")\n", "ax.set_zlabel(\"$y$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Çok Değişkenli Lineer Regresyon: Genel İfade #\n", "\n", "Çok değişkenli doğrusal regresyonu iki değişken örneği üzerinden gördük. İfadeyi $m$ sayıda değişkene genelleştirmek istersek yukarıdaki ifadeleri aşağıda verildiği şekliyle değiştirmemiz gerekir.\n", "\n", "$$ y = a_0 + a_1~x_1 + a_2~x_2 + ... + a_m~x_m $$\n", "\n", "şeklinde ifade edilebilir. Yine minimize edilmek istenen parametre ($S_r$) verinin modelden fark karelerinin toplamıdır ve modelin hatasını (başarımının derecesini) sayısallaştırır. \n", "\n", "$$ S_r = \\sum_{i=1}^{n} (y_{i} - a_0 - a_1 x_{1,i} - a_2 x_{2,i} - ... - a_m x_{m,i})^2 $$\n", "\n", "Yine tipik bir minimizasyon problemi olan bu problemin çözümünde minimize edilmek istenen fonksiyonun ($S_r$), değişkenlerine ($a_0$, $a_1$, $a_2$, ... , $a_m$) göre birinci türevi alınıp sıfıra eşitlenerek minimum değeri aranır. Fonksiyon ifadesindeki $x_{1,i}$, $x_{2,i}$, ..., $x_{m,i}$ ve $y_i$ değerleri yine değişmez ölçümler olduğundan çözüm parabol uyumlamaya çok benzerdir. Sonuçta aşağıda matris formunda ifade edilen lineer denklem sistemi elde edilmiş olur ve bu sistem çözülerek en iyi modelin parametreleri elde edilir.\n", "\n", "\\begin{equation*}\n", "\\begin{bmatrix}\n", "n & \\Sigma x_{1,i} & \\Sigma x_{2,i} & \\cdots & \\Sigma x_{m,i} \\\\\n", "\\Sigma x_{1,i} & \\Sigma x_{1,i}^2 & \\Sigma x_{1,i} x_{2,i} & \\cdots & \\Sigma x_{1,i} x_{m,i} \\\\\n", "\\Sigma x_{2,i} & \\Sigma x_{1,i} x_{2,i} & \\Sigma x_{2,i}^2 & \\cdots & \\Sigma x_{2,i} x_{m,i} \\\\\n", "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "\\Sigma x_{m,i} & \\Sigma x_{1,i} x_{m,i} & \\Sigma x_{2,i} x_{m,i} & \\cdots & \\Sigma x_{m,i}^2 \\\\\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "a_0 \\\\\n", "a_1 \\\\\n", "a_2 \\\\\n", "\\vdots \\\\\n", "a_m\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "\\Sigma y_i \\\\\n", "\\Sigma x_{1,i} y_i \\\\\n", "\\Sigma x_{2,i} y_i \\\\\n", "\\vdots \\\\\n", "\\Sigma x_{m,i} y_i \\\\\n", "\\end{bmatrix}\n", "\\end{equation*}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# En Küçük Kareler Minimizasyonu Yöntemi Genel İfade #\n", "\n", "En küçük kareler minimizasyonu yöntemi gördüğünüz gibi değişken sayısı ve uyumlanan fonksiyonun derecesi ya da parametre sayısından bağımsız olarak benzer bir yol izlemektedir. Bu benzer yol aşağıdaki genel ifadenin elde edilmesini ve bu şekilde en küçük kareler minimizasyonunun tüm durumalr için tek bir ifadeyle verilmesini sağlar.\n", "\n", "$$ y = a_0 ~z_0 + a_1~z_1 + a_2~z_2 + ... + a_m~z_m $$\n", "\n", "Bu ifade \n", "\n", "$ z_0 = 1, z_1 = x, z_2 = .. = z_m = 0 $ için 1-boyutlu lineer regresyon,\n", "$ z_0 = 1, z_1 = x_1, z_2 = x_2, ... , z_m = x_m $ için m-boyutlu lineer regresyon,\n", "$ z_0 = 1, z_1 = x, z_2 = x^2, ... , z_m = x^m $ için m. dereceden polinom regresyonuna\n", " \n", "karşılık gelir. Çözüm her zaman aynı şekilde aranır.\n", "\n", "$$ Z x A = y $$ \n", "\n", "şeklinde matris formunda verilen denklem sisteminde \n", "\n", "\n", "\\begin{equation*}\n", "[Z] = \n", "\\begin{bmatrix}\n", "z_{01} & z_{11} & \\cdots & z_{m1} \\\\\n", "z_{02} & z_{12} & \\cdots & z_{m2} \\\\\n", "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "z_{0n} & z_{1n} & \\cdots & z_{mn}\n", "\\end{bmatrix}\n", "\\begin{bmatrix}\n", "a_0 \\\\\n", "a_1 \\\\\n", "\\vdots \\\\\n", "a_m\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "y_1 \\\\\n", "y_2 \\\\\n", "\\vdots \\\\\n", "y_n \\\\\n", "\\end{bmatrix}\n", "\\end{equation*}\n", "\n", "$A$ katsayılar matrisinin bulunabilmesi için $Z$ matrisinin tersi ($Z^{-1}$) alınır ve her iki taraf sol taraftan $Z^{-1}$ matrisi ile çarpılır. Sonuç $A = Z^{-1} Y$ matrisidir.\n", "\n", "Korelasyon katsayısı: \n", "$$ r^2 = \\frac{S_t - S_r}{S_t} = 1 - \\frac{S_t}{S_r} $$\n", "\n", "$S_t$ aritmetik ortalama ($\\bar{y} = \\frac{\\Sigma y_i}{n}$) model olarak kabul edildiğinde ondan fark karelerin toplamı, $S_r$, en küçük kareler minimizasyonu tekniğiyle elde edilen modelden ($\\hat{y}$) fark karelerin toplamını göstermek üzere;\n", "\n", "$$ r^2 = 1 - \\frac{\\Sigma (y_i - \\hat{y_i})^2}{\\Sigma (y_i - \\bar{y_i})^2} $$\n", "\n", "ile verilir.\n", "\n", "Tahmin üzerindeki standart hata: \n", "\n", "$$ S_{y/x} = \\sqrt{\\frac{\\Sigma (y_i - \\hat{y_i})^2}{n - (m+1)}} $$\n", "\n", "İndirgenmiş $\\chi^2$:\n", "\n", "$$ \\chi^2 = \\sum_{i=1}^{n} \\frac{(\\hat{y_i} - y_i)^2}{e_i^2} $$\n", "\n", "olmak üzere\n", "\n", "$$ \\chi^2_\\nu = \\frac{\\chi^2}{(n - m)} $$\n", "\n", "ile verilir." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Doğrusal Olmayan En Küçük Kareler Yöntemi #\n", "## Levenberg Marquardt Algoritması ##\n", "\n", "Matematik ve bilgisayar bilimlerinde Levenberg-Marquardt Algoritması (LMA) , sönümlenmiş en küçük kareler yöntemi (Damped Least Squares, DLS) olarak da bilinir. m adet gözlemsel veriyi n tane bilinmeyenle modellemek üzere başlangıçta doğrusal bir model üretilip, modelin iterasyonlar sonucu iyileştirilerek doğrusal olmayan bir modelle değiştirilmesi prensibine dayanır.\n", "\n", "$(x_0, y_0)$, $(x_1, y_1)$, $(x_2, y_2)$, … , $(x_m, y_m)$ gözlemsel noktalarımız olsun. Bu noktaları $f(x, \\beta)$ gibi bir fonksiyonla uyumlamak (fit etmek) istiyor olalım. Burada $m \\ge n$ olmak üzere $\\beta = (\\beta_0, \\beta_1, … , \\beta_{n-1}$) fonksiyonun parametreleridir ve fonksiyon $x$ ile birlikte bu parametrelere de bağlıdır. Algoritmada amaç en küçük kareler yöntemine benzer şekilde en küçük artıkları verecek $S_r = \\Sigma (y_i - f(x_i, \\beta))^2$ minimizasyonunun gerçekleşeceği $\\beta$ parametrelerini elde etmektir.\n", "\n", "Tıpkı en küçük kareler yönteminde olduğu gibi yine $S_r$'nin bu kez $\\beta_j$ ($j = 0, 1, 2, … , n$) adını alan fit fonksiyonu parametrelerine göre türevleri alınır ve 0'a eşitlenir. Doğrusal olmayan bir sistemde bu türevler hem $x$'e hem de $\\beta$'ya bağlı olacağından doğrudan çözüm elde etmek yerine probleme uygun seçilmiş başlangıç parametreleri ile başlanır ($\\beta$'lara başlangıç değerleri verilir) ve iteratif bir yaklaşımla her seferinde $S_r$ hesaplanıp, daha öncekilerle ile karşılaştırılarak en küçük $S_r$'yi verecek parametre seti aranır. $\\beta$ her seferinde $\\beta + \\delta$ ile değiştirilir. $\\delta$'yı bulmak için $f(x_i, \\beta + \\delta)$ doğrusallaştırılmak üzere Taylor serisine açılır ($f(x_i, \\beta + \\delta) ≈ f(x_i, \\beta) + J_i\\delta$, $J_i$: Jakobiyen). $S_r$'nin $\\beta$'ya görevi türevi 0 olduğunda $S_r(\\beta + \\delta) = \\Sigma (y_i – f(x_i, \\beta) – J_i\\delta)^2$ toplamının $\\delta$'ya göre türevi de 0 olmalıdır ki minimizasyon sağlanmış olsun. Vektörel notasyonla ($J^T J) \\delta = J^T |(y_i – f(\\beta)|$ denklemi $\\delta$ için çözülür ve $\\beta$ parametreleri bulunur. \n", "\n", "Levenberg'in algoritmaya getirdiği yenilik bir sönümleme (ing. damping) parametresi eklemiş olmasıdır. Böylece denklem ($J^T J + \\lambda I) \\delta = J^T |(y_i – f(\\beta)|$ şeklini alır. Sönümleme parametresi $\\lambda$, her bir iterasyon adımında değiştirilir ve $\\beta + \\delta$ 'nın belirenen bir tolerans değerinin altına indiği anda iterasyon durdurulur. En son $\\beta$ vektörü fonksiyonun parametrelerini verir.\n", "\n", "Doğrusal olmayan (non-lineer) en küçük kareler problemini çözmek üzere geliştirilmiş tek yöntem Levemberg-Marquardt algoritması değildir. [Trust Region Reflective Algorithm](https://nmayorov.wordpress.com/2015/06/19/trust-region-reflective-algorithm/) ya da bu algoritmanın varyantları olan [dogleg ve Steihaug](http://www.numerical.rl.ac.uk/people/nimg/course/lectures/raphael/lectures/lec7slides.pdf) algoritmaları da doğrusal olmayan en küçük kareler probleminin çözümünde kullanılmaktadır.\n", "\n", "## Levenberg-Marquardt Algoritmasına Yönelik Bazı Uyarılar: ##\n", "\n", "Her ne kadar sık kullanılsa da doğrusal olmayan en küçük kareler minimizasyonu probleminin çözümünde Levenberg-Marquardt algoritmasına yönelik olarak aşağıdaki uyarılar dikkate alınmalıdır:\n", "\n", "* Algoritma verilen başlangıç parametrelerine en yakın yerel minimumu bulması nedeniyle dezavantajlıdır.\n", "\n", "* Bir süre benzer (ya da yakın) değerler alan fonksiyonlarda minimizasyon konusunda problemlidir (bu noktalarda zig-zag yapar!).\n", "\n", "* Bazı durumlarda (başlangıç parametrelerinin çözüme yakınlığı) oldukça hızlı çalışırken, bazı durumlarda en iyi fiti verecek parametreleri bulmak konusunda oldukça yavaş olabilir." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# scipy.optimize.leastsq ile Uyumlama #\n", "\n", "`scipy.optimize` modülü eğri uyumlama için pek çok fonksiyon ve alt modül sağlar. [`leastsq`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html) bu fonksiyonlardan biridir ve doğrusal olmayan (non-lineer) en küçük kareler problemini çözerek verilen veriye uyan en iyi eğriyi bulur. Buna ek olarak `scipy.optimize` paketi sadece Levenberg-Marquardt algoritmasıyla değil, Trust Region Reflective algoritması ve dogleg algoritmasıyla da çözüm arayabilen [`least_squares`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) adlı bir başka fonksiyonu daha içermektedir. Bu fonksiyonda çözüm için istenen algoritma, $method$ parametresi sırasıyla $lm$, $trf$ ya da $dogbox$ değerlerine ayarlanarak seçilebilir. Klasik en küçük kareler yönteminden farklı olarak daha hızlı çalışır ve yerel minimumlara saplanma olasılığı daha azdır. Ayrıca her iki fonksiyonla sadece doğru, parabol ya da istenen dereceden başka bir polinom değil, istenen her türden fonksiyon uyumlanabilir. \n", "\n", "Bir örnekle `scipy.optimize.leastsq` fonksiyonunun kullanılışını görelim. Aşağıdaki kodun başında $x$ ve $y$ `numpy` dizileriyle verilen veri setine Levenberg-Marquardt algoritmasıyla bir doğru bir de Parabol fit etmek istiyor olalım. Yapmamız gereken leastsq fonksiyonuna gözlemsel $(x,y)$ değerleri ile hesaplananlar arasındaki hatayı veren bir $HataFonksiyonu$, birer başlangıç parametre seti ($a_0 + a_1 x$ doğrusu için $\\beta = (\\beta_0, \\beta_1) = (a_0, a_1)$ parametreleri, $a_0 + a_1 x + a_2 x^2$ parabolü için $\\beta = (\\beta_0, \\beta_1, \\beta_2) = (a_0, a_1, a_2)$ parametreleri olarak sırasıyla $(1.0, 2.0)$ ve $(1.0, 2.0, 3.0)$, ve veri setimizi $(args = (x,y))$ içeren fonksiyon argümanlarını vermektir. Levenberg-Marquardt algoritması bu parametreler için en iyi değerleri hesaplar ve bu parametreleri doğru ve parabol denklemlerinde yerine koyarak fitlerimizi elde etmiş oluruz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import leastsq\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "# veri\n", "x=np.array([1.0,2.5,3.5,4.0,1.1,1.8,2.2,3.7])\n", "y=np.array([6.008,15.722,27.130,33.772,5.257,9.549,11.098,28.828])\n", "\n", "# Dogru ve Parabol fitleri icin birer lamda fonksionu yaratalim\n", "# tpl burada fitin parametrelerini (a0, a1, ) tutan bir demet degisken\n", "fonkDogru = lambda beta,x : beta[0]*x+beta[1]\n", "fonkParabol = lambda beta,x : beta[0]*x**2+beta[1]*x+beta[2]\n", "# fonk fit etmek istedigimiz fonksiyonun yerini alsin\n", "# Bunu istedigimiz gibi degistirebiliriz.\n", "#------------------------------------------------\n", "# Dogrusal Fit\n", "#-------------------------------------------------\n", "fonk = fonkDogru\n", "# HataFonksiyonu, fit fonksiyonumuzla hesaplanan y ile gozlemsel y arasindaki\n", "# farki versin (toplami S)\n", "HataFonksiyonu = lambda beta,x,y: fonk(beta,x)-y\n", "# tplBaslangicD Levenberg-Marquard algoritmasini baslatmak icin baslangic\n", "# parametre setimiz (beta), dogrusal oldugunu gostermek uzere D kullandim\n", "# beta parametre seti a0 (sabit) ve a1 (egim)'den olusuyor --> a0 + a1*x\n", "beta_baslangicD = (1.0,2.0)\n", "# scipy.optimize.leastsq baslangicD ile baslatilan beta parametre setinin\n", "# HataFonksiyonu toplamini (Sr) veren degerlerini bulur\n", "# HataFonksiyonu = yfit - yIterasyon\n", "beta_sonD, basarim = leastsq(HataFonksiyonu,beta_baslangicD[:],args=(x,y))\n", "print(\" Linear Fit \",beta_sonD)\n", "xxD=np.linspace(x.min(),x.max(),50)\n", "yyD=fonk(beta_sonD,xxD)\n", "#------------------------------------------------\n", "# Parabol Fiti\n", "#-------------------------------------------------\n", "fonk = fonkParabol\n", "beta_baslangicP = (1.0,2.0,3.0)\n", "beta_sonP, basarim = leastsq(HataFonksiyonu,beta_baslangicP[:],args=(x,y))\n", "print(\"Parabol Fiti: \" ,beta_sonP)\n", "xxP = xxD\n", "yyP = fonk(beta_sonP,xxP)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.plot(x,y,'ro',label=\"veri\")\n", "plt.plot(xxD,yyD,'b-',label=\"dogru\")\n", "plt.plot(xxP,yyP,'g-',label=\"parabol\")\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# scipy.optimize.curve_fit ile Uyumlama #\n", "\n", "Tıpkı `scipy.optimize.least_squares` gibi [`curve_fit`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) fonksiyonu da doğrusal olmayan en küçük kareler minimizasyonu problemini çeşitli yöntemlerle çözer. Yapısı ve uygulama şekli ona oranla biraz farklı ve daha esnektir.\n", "\n", "Bu fonksiyonun kullanımını örneklemek için bu kez doğru ya da parabol uyumlamadan farklı bir fonksiyon (üstel fonksiyon) uyumlayalım. Bu amaçla öncelikle üstel fonksiyonla uyumlanabilecek sentetik bir veri seti oluşturalım. Aşağıdaki fonksiyon bu amaçla yazılmıştır. Öncelikle $[0,4]$ arasında 50 nokta içeren $x$ dizisi $y = 2.5 e^{-1.3x} + 0.5$ fonksiyonuna (veri_olusturma) sağlanarak bir $y$ dizisi oluşturulmakta, daha sonra rastgele gürültü eklemek üzere bu dizi `numpy.random.normal` fonksiyonu ile normal (Gaussyen) bir dağılımdan rastgele örneklemeyle oluşturulan y_gurultu dizisiyle toplanmaktadır." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "veri_olusturma = lambda x,a,b,c: a * np.exp(-b * x) + c\n", "x = np.linspace(0,4)\n", "y = veri_olusturma(x,2.5, 1.3, 0.5)\n", "# Veri setimize biraz Gaussyen (normal dagilan) gurultu ekleyelim\n", "np.random.seed(1729)\n", "y_gurultu = 0.2 * np.random.normal(size=x.size)\n", "y += y_gurultu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Şimdi bu oluşturduğumuz sentetik veri setini `curve_fit` fonkisyonuyla uyumlamaya çalışalım. Bu amaçla `curve_fit` fonksiyonuna uyumlayacağı fonksiyonun türünü ve verimizi argüman olarak sağlamamız gerekmektedir. Uyumlamak istediğimiz fonksiyon üstel bir fonksiyon olduğu için aslında sentetik veri oluşturmak için kullandığımız $lambda$ fonksiyonuyla (veri_olusturma) aynı olacaktır. O nedenle yeni bir fonksiyon oluşturmak yerine bu fonksiyonu sağlayalım ve `curve_fit` 'in sentetik verimiz için en iyi $a$, $b$ ve $c$ değerlerini bulmasını sağlayalım. Biz bu değerlerin aslında sırasıyla $2.5$, $1.3$ ve $0.5$'e yakın olmasını bekliyoruz. Ancak veriye gürültü de eklediğimiz için `curve_fit` tarafından bulunan değerler doğal olarak aynı çıkmamalıdır." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import curve_fit\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "param, kovar = curve_fit(veri_olusturma, x, y)\n", "print(\"Veriye uyumlamanan ustel fonksiyonun parametreleri: \", param)\n", "a = float(param[0])\n", "b = float(param[1])\n", "c = float(param[2])\n", "print(\"Model: {:.2f} e^(-{:.2f}x) + {:.2f}\".\\\n", " format(a,b,c))\n", "plt.plot(x,y,'ro',label=\"veri\")\n", "y_model = a*np.exp(-b*x) + c\n", "plt.plot(x, y_model,'b-',label='model')\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# LMFIT Paketiyle Modelleme #\n", "\n", "Levenberg-Marquardt algoritmasıyla daha sofistike ve esnek uyumlamalar için geliştirilmiş olan `LMFIT` paketi, scipy.optimize modülü üzerine kurulu ve onun optimizasyon ve eğri uyumlama yeteneklerini daha ileriye taşıyan bir paket. İndirmek, kurmak ve kullanımı için [bkz.](https://lmfit.github.io/lmfit-py/).\n", "\n", "Aşağıda önce `scipy.optimize.curve_fit` fonksiyonu kullanılarak çözülen bir uyumlama problemi örneklenmiştir. Öncelikle $[-10,10]$ aralığında 50 nokta içeren bir $x$ dizisi `numpy.linspace` fonksiyonuyla oluşturulmaktadır. Daha sonra sentetik veriyi oluşturmak için bu dizi $gauss$ fonksiyonuna $A = 2.33$ genlik, $\\mu = 0.21$ ortalama değer, $\\sigma = 1.51$ standart sapma değerleriyle gönderilmekte ve $y = A e^{-(x-\\mu)^2 / (2\\sigma)}$ fonksiyonuyla sentetik y verisi oluşturulmaktadır. Daha sonra bu veriye ortalma değeri 0.0, standart sapması 0.2 olan bir Gauss dağılımından rastgele seçilen 50 değer eklenmekte, bu şekilde veri üzerinde sentetik bir gürültü (saçılma) oluşturulmaktadır.\n", "\n", "Bu veri `scipy.optimize.curve_fit` fonksiyonuyla yine daha önceki gibi veriyi oluşturmak için kodlanan $gauss$ fonksiyonu ile birlikte sağlanarak modellenmek istenmektedir. Gauss fonksiyonunun başlangıç parametreleri $A = 1.0$, $\\mu = 0.0$, $\\sigma = 1.0$ olarak seçilmiştir. Biz bu parametrelerin fit sonunda $A = 2.33$, $\\mu = 0.21$ ve $\\sigma = 1.51$ değerlerine yakınsamasını bekliyoruz; zira bu veriyi bu değerleri kullanarak fit için de kullandığımız $gauss$ fonksiyonuyla kendimiz oluşturduk." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.optimize import curve_fit\n", "def gauss(x, genlik, ortalama, stsapma):\n", " return genlik*np.exp(-(x-ortalama)**2/(2*stsapma))\n", "x = np.linspace(-10, 10)\n", "y = gauss(x, 2.33, 0.21, 1.51) + np.random.normal(0, 0.2, len(x))\n", "baslangicPar = [1, 0, 1]\n", "param, kovar = curve_fit(gauss, x, y, p0 = baslangicPar)\n", "print(\"Veriye uyumlamanan ustel fonksiyonun parametreleri: \", param)\n", "A = float(param[0])\n", "mu = float(param[1])\n", "sigma = float(param[2])\n", "print(\"Model: {:.2f} e^(-(x-{:.2f}) / 2*{:.2f})\".\\\n", " format(A,mu,sigma))\n", "x2 = np.linspace(-10,10,100)\n", "yfit = param[0]*np.exp(-(x2 - param[1])**2/param[2])\n", "plt.plot(x,y,'ro',label=\"veri\")\n", "plt.plot(x2,yfit,'b-',label=\"model\")\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gördüğümüz gibi bulduğumuz değerler veriyi oluşturmak için kullandığmız değerlere yakın. Yani `scipy.optimize.curve_fit` iyi bir sonuç vermiş görünüyor. Ancak bu durumda `LMFIT` paketinin getirdiği bazı özelliklerden mahrum kalırız. Bunlardan biri `LMFIT`'in sağladığı $Model$ metodudur. Bu metodun ürettiği modelin parametrelerinin sırasını bilmeniz gerekmez, adlarını bilmeniz yeterlidir. Şimdi aynı soruyu `LMFIT` modülü fonksiyonlarından yararlanarak çözelim. $Model$ fonksiyonuna, veriyi hem oluşturmak hem de modellemek üzere kullandığımız $gauss$ fonksiyonunu sağlayalım ve öncelikle modelimizin parametrelerini görelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import lmfit as lm\n", "gmodel = lm.Model(gauss)\n", "print(gmodel.param_names)\n", "print(gmodel.independent_vars)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gördüğünüz gibi modelimizin parametreleri kendi yazdığmız $gauss$ fonksiyonunda kendi tanımladığımız parametreler; bağımsız değişkeni ise yine bizim verdiğimiz isimle $x$. Şimdi modelimizi veriye uyumlayalım (fit edelim) ve fitimizin parametrelerini fit_report fonksiyonuyla inceleyelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sonuc = gmodel.fit(y, x=x, genlik=1, ortalama=0, stsapma=1)\n", "print(sonuc.fit_report())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`LMFIT` [dokümantasyoundan](https://lmfit.github.io/lmfit-py/) da anlaşılacağı üzere oldukça gelişmiş bir uyumlama ortamı sunar. Uyumlama temel parametrelerle yapıldığında görece kolay olmakla birlikte `LMFIT` pek çok başka girdi parametresine daha sahiptir ve bu opsiyonel parametreler uyumlama için oldukça esnek bir yapı sunar. Tüm bu detayları dokümantasyon yoluyla öğrenip, yerli yerince kullanabilirsiniz. Ancak bu aşamada bu detay düzeyiyle yetineceğiz. fit_report() fonksiyonunun sağladğı rapor parametrelerine baktığınızda model parametrelerinin değerleri ve hatalarının (yüzde düzeyleriyle birlikte) yanı sıra $\\chi^2$, indirgenmiş $\\chi^2_\\nu$ (reduced chi-square) gibi uyumlamanın başarısını gösteren istatistiklerle birlikte Akaike ve Bayesian ölçütleri gibi model karşılaştırması konusunda öğreneceğiniz diğer parametrelerin değerlerini de görebilirsiniz. \n", "\n", "Son olarak yapılan uyumlamayı grafik üzerinde gösterelim. Sonuçtaki modele daha önce report_fit metodunu gördüğünüz $sonuc$ nesnesinin best_fit, başlangıç parametreleriyle oluşan modele ise init_fit metoduyla ulşabilir; uyumlama parametrelerinin değerlerini params sözlüğünün ilgili parametreye karşılık gelen anahtarının value metoduyla öğrenebilirsiniz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "plt.plot(x, y, 'ro', label=\"veri\") # veri\n", "plt.plot(x, sonuc.init_fit, 'b--', label=\"başlangic modeli\") # A=1, mu=0, sigma=1 icin fit\n", "plt.plot(x, sonuc.best_fit, 'g-', label=\"sonuc model\") # Veri setine uyumlanan en iyi fit\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(loc=\"best\")\n", "plt.show()\n", "# Fitimizin parametrelerini ekrana yazdiralim\n", "print(\"En iyi uyumlamanin genligi = {:.2f}, ortalama degeri = {:.2f}, st.sapma = {:.2f}\".\\\n", " format(sonuc.params['genlik'].value,sonuc.params['ortalama'].value,sonuc.params['stsapma'].value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Örnek: Küresel Isınmanın Ciddiyeti ##\n", "\n", "Karbondioksit seviyesinde 1960'dan bu yanaki değişim Keeling Eğrisi adı verilen bir eğriyle ifade edililr. Söz konusu değişim küresel ısınmaya neden olduğu gerekçesiyle endişe konusudur. Aşağıda 1960'dan bu yana geçen süre içerisinde onar yıllık ölçüm sonuçları milyon parçacık başına verilmiştir. Keeling Eğrisinin bir parabolle (2. dereceden bir polinomla) temsil edilebileceğini düşünerek ve eldeki veriyi kullanarak eğriye en uygun formülü bulunuz. Keeling Eğrisi yapısını korursa 2030 yılı sonunda parçacık başına kaç karbondioksit molekülü düşer hesaplayınız.\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Yıl $CO_2$
0317
10326
20338
30354
40369
50390
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Öncelikle verimizi bir grafiğe döküp nasıl bir modelle veriyi en iyi uyumlayabileceğimizi düşünelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "t = np.arange(0,51,10.)\n", "CO2 = np.array([317,326,338,354,369,390])\n", "plt.plot(t,CO2,'ro')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Değişim lineer gibi gözükse de bir miktar eğiklik kendini hissettirdiğinden öncelikle bir parabolik model ($a_0 + a_1~x + a_2~x^2$) varsayalım. Daha önce olduğu gibi modelimizin parametrelerini bir görelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import lmfit as lm\n", "parabol = lambda x,a0,a1,a2: a0 + a1*x + a2*x**2\n", "parmodel = lm.Model(parabol)\n", "print(parmodel.param_names)\n", "print(parmodel.independent_vars)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Modelimizin ($parmodel$) üç parametresi ($a_0$, $a_1$ ve $a_2$) ve bir de bağımsız değişkeni ($x$) var. Şimdi bu modeli verimize fit etmek üzere modelimiz ($parmodel$) üzerinde tanımlı `fit` metodunu kullanalım. Metoda verimizi önce $y = CO2$, sonra $x = t$ ile verdiğimize ve başlangıç parametresi olarak da göz kararı belirlediğimiz $a_0 = 300$, $a_1 = 2$ ve $a_2 = 1$ değerlerini argüman olarak verdiğimize dikkat ediniz. Uyumlama sonucunu $sonucpar$ değişkenine alıyoruz. Bu değişken bir $lmfit$ model sonucu nesnesi olarak böylece tanımlanmış oluyor. Bu nesnenin uyumlamaya ilişkin istatistiksel bilgileri alabileceğiniz bir `fit_report()` metodu ve parametrelerini almanız için de `params` adında bir sözlük değişkeni tanımlı. Öncelikle uyumlamaya ilişkin istatistiksel bilgileri görelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sonucpar = parmodel.fit(CO2, x=t, a0 =300, a1=2, a2=1)\n", "print(sonucpar.fit_report())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bu bilgilerden özellikle $\\chi^2$ ve indirgenmiş $\\chi^2_\\nu$ istatistiğinin uyumlama başarısını göstermekte kullanıldığını hatırlayalım. Akaike ve Bayes parametreleri de uyum başarısını verir ancak bu parametrelerin nasıl tanımlandığını daha sonra model karşılaştırmasına yönelik dersimizde göreceğiz. Ayrıca parametrelerimize ilişkin değerler, korelasyon katsayıları ve hataları da raporda sunulmuş durumdadır. Şimdi modelimizi verimizin üzerine çizdirelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "plt.plot(t, CO2, 'ro', label=\"veri\") # veri\n", "plt.plot(t, sonucpar.init_fit, 'b--', label=\"başlangic modeli\") # a0=300, a1=2, a2=1\n", "plt.plot(t, sonucpar.best_fit, 'g-', label=\"sonuc model\") # Veri setine uyumlanan en iyi fit\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(loc=\"best\")\n", "plt.show()\n", "# Fitimizin parametrelerini ekrana yazdiralim\n", "print(\"En iyi uyumlama: {:.2f} + {:.2f}x + {:.2f}x^2\".\\\n", " format(sonucpar.params['a0'].value,sonucpar.params['a1'].value,sonucpar.params['a2'].value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Görüldüğü üzere başlangıç parametreleriyle oluşturulan eğri (mavi kesikli) verimizi uyumlamaktan uzakken, `lmfit` 'in bu başlangıç parametrelerinden hareketle uyumladığı parabol (yeşil sürekli) başarılı görünmektedir. $a_2$ katsayısının değerinin $0.01$ olarak bulunmuş olması bize doğrusal bir modelin de belki başarılı olacağını düşündürtebilir. Bunun için verimize bir de doğru modeli uyumlamaya çalışalım. Bu kez iki parametreli bir doğru fonksiyonu ($dogru$) tanimlayacağız ve bunu bir `lmfit` modeline `Model` fonkisyonuyla dönüştüreceğiz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dogru = lambda x,a0,a1: a0 + a1*x\n", "dogrusalmodel = lm.Model(dogru)\n", "print(dogrusalmodel.param_names)\n", "print(dogrusalmodel.independent_vars)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ardından modelimizin doğrusal olmayan en küçük kareler yöntemlerinden Levenberg-Marquardt algoritmasıyla uyumlanması sonucunda oluşan en iyi modelimize ilişkin parametreleri sonuc_dogrusal model sonucu değişkenine alalım ve sonuç modelimizin parametrelerini görelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sonuc_dogrusal = dogrusalmodel.fit(CO2, x=t, a0 =300, a1=1.0)\n", "print(sonuc_dogrusal.fit_report())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gördüğünüz gibi $\\chi^2_\\nu >> 1$ durumu oluşmakta, doğru modelin parabole göre verinin üzerindeki belirsizliği açıklamakta daha başarısız olduğu görünmektedir. Ancak her iki modeli bir de grafik üzerinde karşılaştıralım. Gelecek derste model karşılaştırmasına yönelik olarak yapılması gerekenleri ayrıntılı olarak göreceğiz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "plt.plot(t, CO2, 'ro', label=\"veri\") # veri\n", "plt.plot(t, sonucpar.best_fit, 'b--', label=\"parabol modeli\") # a0=300, a1=2, a2=1\n", "plt.plot(t, sonuc_dogrusal.best_fit, 'g-', label=\"dogrusal model\") # Veri setine uyumlanan en iyi fit\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "plt.legend(loc=\"best\")\n", "plt.show()\n", "# Fitimizin parametrelerini ekrana yazdiralim\n", "print(\"En iyi uyumlama: {:.2f} + {:.2f}x\".\\\n", "\n", " format(sonuc_dogrusal.params['a0'].value,sonuc_dogrusal.params['a1'].value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Eğri Uyumlamanın İnterpolasyon ve Ekstrapolasyon Amacıyla Kullanımı #\n", "\n", "Bir matematiksel ya da fiziksel model, bağımsız bir değişkenin herhangi bir değeri için bağımlı değişkene ilişkin bir sonuç vereceği için interpolasyon ya da ekstrapolasyon amacıyla da kullanılabilir. Örneğimizde atmosferimizdeki milyon parçacık başına $CO_2$ molekülü sayısının Keeling eğrisiyle artması koşuluyla gelecekte alabileceği değerleri ya da 1960'dan önceki değerlerini hesaplamak isteyebiliriz. Eldeki veriden yola çıkarak gelecek ve geçmişe yapılan bu projeksiyon ekstrapolasyon olarak bilinir. Bilimde \"tehlikeli\" bir işlem olmakta birlikte sıklıkla uygulanır. Tehliklesi verinin alındğı aralıktaki fonksiyonun gelecek ve geçmiş için de geçerli olup olmayacağının bilinmemesidir. Diğer taraftan verinin alındığı aralıkta bulunan ancak veri noktasının bulunmadığı bir nokta için de hesap yapılabilir. Daha önce çeşitli yöntemlerini gördüğünüz bu yöntem interpolasyon olarak bilinir ve eğri uyumlamayla elde edilen bir matematiksel model interpolasyon amacıyla da kullanılabilir.\n", "\n", "Örneğimizden hareketle 2030 ve 2005 yılları için atmosferimizdeki milyon parçacık başına, parabolik modelimizin kaç $CO_2$ molekülü öngördüğünü hesaplayarak sırasıyla ekstrapolasyon ve interpolasyon kavramlarını örnekleyelim." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# 2030 yıilinda modelimiz bize ne kadar bir CO2 parcacik sayisi\n", "# olacagini ongoruyor\n", "# Ekstrapolasyon\n", "x2030 = 2030 - 1960 # 1960dan beri gecen yil sayisi\n", "a0 = sonucpar.params['a0'].value\n", "a1 = sonucpar.params['a1'].value\n", "a2 = sonucpar.params['a2'].value\n", "CO2_2030 = a0 + a1*x2030 + a2*x2030**2\n", "print(\"Parabol modelimiz {:d} yilinda milyon parcacik basına {:.2f} CO2 parcacigi olacagini ongormektedir.\".\\\n", " format(x2030+1960,CO2_2030))\n", "# Interapolasyon\n", "x2005 = 2005 - 1960 # 1960dan beri gecen yil sayisi\n", "CO2_2005 = a0 + a1*x2005 + a2*x2005**2\n", "print(\"Parabol modelimiz {:d} yilinda milyon parcacik basına {:.2f} CO2 parcacigi oldugunu ongormektedir.\".\\\n", " format(x2005+1960,CO2_2005))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Kaynaklar #\n", "\n", "* Measurements and their Uncertainties, Ifan G. Hughes & Thomas P.A. Hase, Oxford University Press, 2010\n", "\n", "* Data Reduction and Error Analysis for the Physical Sciences, 3rd ed., Philip R.Bevington & D. Keith Robinson, 2003, McGraw Hill Higher Education\n", "\n", "* Numerical Analysis Using Matlab and Excel 3rd ed., Steven T. Karris, Orchard Publications, 2007\n", "\n", "* Numerical Methods for Engineers 6th ed., Steven C. Chapra, Raymond P. Canale, McGraw Hill, 2010\n", "\n", "* Numerical Methods, Rao V. Dukkipati, New Age International, 2010\n", "\n", "* Ascombe, F.J., 1973, The American Statistician, Vol. 27, No. 1. (Feb., 1973), pp. 17-21.\n", "\n", "* [numpy.polyfit Dokümantasyonu](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html)\n", "\n", "* [numpy.linalg Lineer Cebir Modülü Dokümantasyonu](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html)\n", "\n", "* [numpy.matrix Fonksiyonu Dokümantasyonu](https://docs.scipy.org/doc/numpy/reference/generated/numpy.matrix.html)\n", "\n", "* [scipy.optimize.curve_fit Dokümantasyonu](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) \n", "\n", "* [scipy.optimize.leastsq Dokümantasyonu](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html)\n", "\n", "* [scipy.optimize.least_squares Dokümantasyonu](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html)\n", "\n", "* [LMFIT Dokümantasyonu](https://lmfit.github.io/lmfit-py/)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Başa Dön](#Bu-derste-neler-öğreneceksiniz?)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.9" } }, "nbformat": 4, "nbformat_minor": 2 }