In [1]: 
from pylab import *
%matplotlib inline
import statsmodels.api as sm
from scipy.optimize import curve_fit

# nastavitve za izris grafov (http://matplotlib.org/1.3.1/users/customizing.html)
rc('text', usetex=True)
rc('font', size=12, family='serif', serif=['Computer Modern'])
rc('xtick', labelsize='small')
rc('ytick', labelsize='small')
rc('legend', fontsize='medium')
rc('figure', figsize=(5, 3))
rc('lines', linewidth=2.0)
rc('axes', color_cycle=['k'])
rc('contour', negative_linestyle='solid')

Nelinearna regresija

Pri prejšnji temi smo si pogledali, kako poiščemo linearno kombinacijo funkcij \(y(x)=\sum_i c_i f_i(x)\), ki najbolje opiše podatke \((x_i, y_i, \epsilon_i)\). Podobno postopamo, če parametri \(c_i\) v modelski funkciji nastopajo nelinearno. Minimiziramo namreč izraz

\[\chi^2=\sum\left(\frac{y_i-f(c_1,c_2,c_3\ldots, x_i)}{\epsilon_i}\right)^2=\mathrm{min}.\]

Pri tem je treba paziti, da, za razliko od linearne regresije, pri nelinearni regresiji ni garancije, da ima izraz \(\chi^2\) en sam minimum. Lahko se pojavi več lokalnih minimumov. Za različne začetne približke za parametre \(c_i\), ki jih pri nelinearni regresiji moramo podati, lahko rutina vrne različne krivulje, odvisno pač, v katerega od lokalnih minimumov je postopek skonvergiral.

Poglejmo si razliko med linearno in nelinearno regresijo na primeru. Na sliki (a) gre prilagajanje premice \(y=kx+n\). Na sliki (b) vidimo, da ima \(\chi^2(k, n)\) en sam minimum. Na sliki (c) prilagajamo podatkom funkcijo \(y=A\cos{\left(2\pi kx\right)}\). \(\chi^2(A, k)\), prikazan na sliki (d), ima več lokalnih minimumov. Če rutini za nelinearno regresijo predlagamo, da začne minimum iskati pri \((A, k)=(0.5, 0.8)\), najde lokalni minimum, če pa predlagamo \((A, k)=(1, 1)\), pa najde globalni minimum.

In [55]: 
fig, ((axa, axb), (axc, axd)) = subplots(2, 2, figsize=(10, 7))

podatki=loadtxt('anscombe.dat')
x1=podatki[:, 0]
y1=podatki[:, 1]
x=arange(3, 20.5, 0.5)
fit1=sm.OLS(y1, sm.add_constant(x1, prepend=True)).fit()
f1=lambda x: fit1.params[0] + fit1.params[1] * x
axa.plot(x1, y1, 'ko');
axa.plot(x, f1(x), 'r');
axa.grid(alpha=0.5)
axa.annotate('(a)', xy=(0.03, 0.92), xycoords='axes fraction')
axa.set_xlim([3, 15])
axa.set_ylim([0, 12])
axa.set_xlabel(r'$x$')
axa.set_ylabel(r'$y$')
axa.legend(['podatki', 
            r'$y=kx+n$' + '\n' + r'$k=' + '%.2f' % fit1.params[1] + r'$, ' + r'$n=' + '%.2f' % fit1.params[0] + r'$'], loc=4); 
            
def chi2(k, n):
    return log(sum((y1 - k * x1 - n) ** 2))
k=arange(-0.5, 1.5, 0.02)
n=arange(-2.0, 8.05, 0.05)
K, N=meshgrid(k, n)
chi2=frompyfunc(chi2, 2, 1)
Z=chi2(K, N)
axb.contourf(K, N, Z, cmap="RdBu_r")
axb.annotate('(b)', xy=(0.03, 0.92), xycoords='axes fraction')
p=axb.contour(K, N, Z, colors='k', linewidths=1.0)
clabel(p, inline=1, fmt='%.1f');
axb.set_xlim([-0.5, 1.5])
axb.set_ylim([-2, 8])
axb.set_xlabel(r'$k$')
axb.set_ylabel(r'$n$')

x1=np.linspace(0, 4.0, 100)
y1=cos(2.0 * pi * x1) + 0.5 * cos(1.5 * pi * x1) + np.random.normal(0, 0.2, 100)
axc.plot(x1, y1, 'ko');
axc.grid(alpha=0.5)
axc.annotate('(c)', xy=(0.03, 0.92), xycoords='axes fraction')
axc.set_xlim([0, 4])
axc.set_ylim([-4.0, 1.5])
axc.set_xlabel(r'$x$')
axc.set_ylabel(r'$y$')

def func(x, A, k):
  return A * cos(2.0 * pi * k * x)
params1, cov=curve_fit(func, x1, y1, p0=(1.0, 1.0))
f=lambda x: params1[0] * cos(2.0 * pi * params1[1] * x)
axc.plot(x1, f(x1), 'r')
params2, cov=curve_fit(func, x1, y1, p0=(0.5, 0.8))
f=lambda x: params2[0] * cos(2.0 * pi * params2[1] * x)
axc.plot(x1, f(x1), 'b')
axc.legend(['podatki', 
            r'$y=A\cos{\left(2\pi kx\right)}$' + '\n' + r'$A=' + '%.2f' % params1[0] + r'$, ' + r'$k=' + '%.2f' % params1[1] + r'$',
            r'$y=A\cos{\left(2\pi kx\right)}$' + '\n' + r'$A=' + '%.2f' % params2[0] + r'$, ' + r'$k=' + '%.2f' % params2[1] + r'$'], loc=4);

def chi2(a, k): 
    return log(sum((y1 - a * cos(2.0 * pi * k * x1)) ** 2))
A=arange(0.0, 1.51, 0.02)
k=arange(0.6, 1.11, 0.02)
A1, k1=meshgrid(a, k)
chi2=frompyfunc(chi2, 2, 1)
Z=chi2(A1, k1)
axd.contourf(A1, k1, Z, 15, cmap="RdBu_r")
axd.annotate('(d)', xy=(0.03, 0.92), xycoords='axes fraction')
p=axd.contour(A1, k1, Z, 15, colors='k', linewidths=1.0)
clabel(p, inline=1, fmt='%.1f');
axd.set_ylim([0.6, 1.1])
axd.set_xlim([0.0, 1.5])
axd.set_ylabel(r'$k$')
axd.set_xlabel(r'$A$')

subplots_adjust(hspace=0.25, wspace=0.15)

Poglejmo si še, kako poiščemo najboljšo eksponentno funkcijo \(y=Ae^{-\lambda x}\) za podatke na grafu (a) spodaj. Pri prejšnji temi smo podatke \(y_i\) najprej logaritmirali in nato poiskali najboljšo premico \(\ln{y}=\ln(A)-\lambda x,\) graf (b). Iz koeficientov te premice smo prebrali \(A\) in \(\lambda\) in narisali ustrezno eksponentno funkcijo (c, rdeča črta). S pomočjo nelinearne regresije se lahko vmesnim korakom izognemo in podatkom prilagodimo kar nelinearno odvisnost \(y=Ae^{-\lambda x}\) (c, modra črta).

In [67]: 
fig, ((axa, axb, axc)) = subplots(3, 1, figsize=(5, 9), sharex='col', sharey='row')
axa.set_xlim([0, 10])
axa.set_ylim([0, 200])
axb.set_ylim([1, 6])
axc.set_ylim([0, 200])
axa.set_ylabel(r'$y$')
axb.set_ylabel(r'$\ln{y}$')
axc.set_ylabel(r'$y$')
axc.set_xlabel(r'$x$')

N=50
x=linspace(0, 10, N)
y=exp(log(200.0 * exp(-0.3 * x)) + np.random.normal(0, 0.05, N))

axa.plot(x, y, 'ko')
axa.legend(['podatki'])
axa.annotate('(a)', xy=(0.03, 0.05), xycoords='axes fraction')

axb.plot(x, log(y), 'ko')
fit=sm.OLS(log(y), sm.add_constant(x, prepend=True)).fit()
f=lambda x: fit.params[0] + fit.params[1] * x
axb.plot(x, f(x), 'r')
axb.legend(['podatki', r'$\ln{y}=\ln{A}-\lambda x$'])
axb.annotate('(b)', xy=(0.03, 0.05), xycoords='axes fraction')

axc.plot(x, y, 'ko')
axc.plot(x, exp(f(x)), 'r')
axc.annotate('(c)', xy=(0.03, 0.05), xycoords='axes fraction')

def func(x, A, lam):
  return A * exp(-lam * x)
params, cov=curve_fit(func, x, y, p0=(100.0, -1.0))
f=lambda x: params[0] * exp(-params[1] * x)
axc.plot(x, f(x), 'b--')

axc.legend(['podatki', 
            r'$y=Ae^{-\lambda x}$ (lin.)' + '\n' + r'$A=' + '%.3g' % exp(fit.params[0]) + r'$, ' + r'$\lambda=' + '%.3g' % -fit.params[1] + r'$', 
            r'$y=Ae^{-\lambda x}$ (nelin.)' + '\n' + r'$A=' + '%.3g' % params[0] + r'$, ' + r'$\lambda=' + '%.3g' % params[1] + r'$']);

subplots_adjust(hspace=0.075, wspace=0.05)

Naloge

  1. Teorija kemijske kinetike napove za sigmoidno krivuljo iz podatkov adrenalin.dat naslednjo odvisnost \(F/F_\textrm{max}=c/(a+c)\), kjer pomeni \(a\) koncentracijo s polovičnim maksimalnim učinkom. S pomočjo nelinearne regresije določi koeficienta \(F_\textrm{max}\) in \(a\). Primerjaj z rezultati naloge 7.2.

  2. V datoteki higgs_to_gammagamma.dat je podan histogram števila kandidatov \(\textrm{H}\rightarrow\gamma\gamma\) kot funkcija izmerjene invariantne mase dveh fotonov (v območju med 100 in 160\(\,\)GeV/\(c^2\)). Prvi stolpec podaja zaporedno število predalčka, drugi pa število kandidatov v tem predalčku. Širina predalčka je 2\(\,\)GeV/\(c^2\). Z nelinearno regresijo določi maso Higgsovega bozona. Izmerjene podatke opiši z naslednjo funkcijo:

    \[F(m)=A\,\textrm{ozadje}(m,c)+B\,\textrm{signal}(m,m_0,\sigma),\]

    kjer porazdelitev ozadja opišemo z eksponentno funkcijo,

    \[\textrm{ozadje}(m,c)=e^{mc},\]

    in porazdelitev signala z Gausovo,

    \[\textrm{signal}(m,m_0,\sigma)=e^{-\frac{(m-m_0)^2}{2\sigma^2}}.\]

    Širino Gausove funkcije nastavi na 1.77\(\,\)GeV/\(c^2\). Ostali parametri, \(A\), \(B\), \(c\) in \(m_0\) so neznani.

  3. V datoteki RefractiveIndex.txt so zbrane meritve lomnega količnika nekega materiala v odvisnosti od valovne dolžine pri različnih temperaturah. Odvisnost lomnega količnika \(n\) od valovne dolžine \(\lambda\) opisuje Sellmeierjeva formula

    \[n(\lambda)^2-1=\sum_{i=1}^{N}\frac{S_i}{1-\left(\frac{\lambda_i}{\lambda}\right)^2},\]

    kjer so \(\lambda_i\) valovne dolžine elektronskih prehodov v opazovanem materialu, \(S_i\) pa pripadajoče uteži. Z nelinearno regresijo določi \(\lambda_i\) in \(S_i\) (a) ob predpostavki, da meritev lahko opišemo z enim samim elektronskim prehodom (\(N=1\)) in (b) ob predpostavki, da lahko meritev opišemo z dvema takima prehodoma (\(N=2\)). Nariši temperaturno odvisnost valovnih dolžin \(\lambda_i\) in uteži \(S_i\).