Search
Curve Fitting

Во подоцнежните тетратки ќе користиме кулмутативни стапки на смрност од ден (како точки во график) и функција (или поточно нашиот модел) што е зависна од некои параметри (смртност, кревети за интензивна нега, респиратори итн.) која истотака ќе ни дава кулмутивни стапки за секој ден, т.е. ќе ги предвидува. Ние потоа ги прилагодуваме кривите (fit the curves), тоа значи дека сакаме да ги најдеме параметрите за нашиот модел што би резултирале со тој да вади предивудвања најблиску до вистински податоци (real-data).

Не знаеме доволно сеуште да навлеземе во математиката зад ова, меѓутоа овој концепт би го објасниле со Гаусова крива која има малку шум, генерираме доста прост пример:

import numpy as np
import pandas as pd
import lmfit
import warnings
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import mpld3
from lmfit.lineshapes import gaussian, lorentzian
from scipy.integrate import odeint
import plotly.graph_objects as go
import plotly.tools as tls
from plotly.offline import plot, iplot, init_notebook_mode
from IPython.core.display import display, HTML
init_notebook_mode(connected = True)
config={'showLink': False, 'displayModeBar': False}
mpld3.enable_notebook()
%matplotlib inline 
pd.options.mode.chained_assignment = None  # default='warn'
warnings.filterwarnings('ignore')

np.random.seed(42)
x = np.linspace(0, 20.0, 1001)

data = (gaussian(x, 21, 6.1, 1.2) + np.random.normal(scale=0.1, size=x.size))  # normal distr. with some noise

figa = go.Figure()

figa.add_trace(go.Scatter(x = x, 
                         y = data, 
                         mode = 'lines',
                         line=dict(color="blue",
                                   width=1),
                         name = "<i>Normal distribution</i>",
                         hovertemplate = "<b>x: %{x}</b><br><b>y: %{y}</b>"))

figa.update_layout(title = 'Нормална распределба;',
                   title_x = 0.5,
                  xaxis_title='x',
                  xaxis=dict(autorange=True, 
                             mirror=False,
                             ticks='outside',
                             showline=True,
                             linecolor='#000',
                             tickvals = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20],
                             tickfont = dict(size=11)),
                  yaxis_title='y',
                  yaxis=dict(autorange=True, 
                             mirror=False,
                             ticks='outside', 
                             showline=True,
                             linecolor='#000',
                             tickvals = [0, 1, 2, 3, 4, 5, 6, 7],
                             tickfont = dict(size=11)),
                  plot_bgcolor='#fff', 
                  width = 500,
                  height = 320,
                  font = dict(size = 10),
                  margin=go.layout.Margin(l=50,
                                         r=50,
                                         b=60,
                                         t=35))

plot(figa, filename = 'fig10_a.html', config = config)
display(HTML('fig10_a.html'))

Следно, нужна ни е функција што зима како аргументи: вредност за $x$ и три параметри што би сакале да ги прилагодиме (во пример тие се p1, p2, p3). Ова функција има интенција да прилагоди податоци (fit the data); библиотеката која ја користиме за fitting (бидејќи сме доста нови и не можеме да правиме потешки fitting) варира еден параметар се додека не најде "адекватно" совпаѓање. Со надеж дека враќа исти параметри со оние кои сме ги генерирале податоците за кривата.

За целото fitting ја користиме библиотеката lmfit. Започнуваме со генерирање на функција (што е само друг збор за модел) кој на влез чита неколку иницијални препоставки за параметрите. Потоа само ги fit-нуваме податоците. Важно е да се напомене дека сите методи за curve-fitting генерално не гарантираат наоѓање на глобален минимум и дека нашите први претпоставки кои одат на влез се премногу важни.

Доколку влезот е векторот: p1: 21, p2: 6.1, p3: 1.2, излезот би бил:

import numpy as np
import pandas as pd
import lmfit
import warnings
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import mpld3
from lmfit.lineshapes import gaussian, lorentzian
from scipy.integrate import odeint
import plotly.graph_objects as go
import plotly.tools as tls
from plotly.offline import plot, iplot, init_notebook_mode
from IPython.core.display import display, HTML
init_notebook_mode(connected = True)
config={'showLink': False, 'displayModeBar': False}
mpld3.enable_notebook()
%matplotlib inline 
pd.options.mode.chained_assignment = None  # default='warn'
warnings.filterwarnings('ignore')

def f(x, p1, p2, p3):
    return gaussian(x, p1, p2, p3)

np.random.seed(42)
x = np.linspace(0, 20.0, 1001)

data = (gaussian(x, 21, 6.1, 1.2) + np.random.normal(scale=0.1, size=x.size))  # normal distr. with some noise

figa = go.Figure()

mod = lmfit.Model(f)
# we set the parameters (and some initial parameter guesses)
mod.set_param_hint("p1", value=10.0, vary=True)
mod.set_param_hint("p2", value=10.0, vary=True)
mod.set_param_hint("p3", value=10.0, vary=True)

params = mod.make_params()

result = mod.fit(data, params, method="leastsq", x=x)  # fitting

figa.add_trace(go.Scatter(x = x, 
                         y = data, 
                         mode = 'lines',
                         line=dict(color="blue",
                                   width=1),
                         name = "<i>Data</i>",
                         hovertemplate = "<b>x: %{x}</b><br><b>y: %{y}</b>"))

figa.add_trace(go.Scatter(x = x, 
                         y = result.best_fit, 
                         mode = 'lines',
                         line=dict(color="orange",
                                   width=2),
                         name = "<i>Best-fit</i>",
                         hovertemplate = "<b>x: %{x}</b><br><b>y: %{y}</b>"))

figa.update_layout(title = 'Нормална распределба; Best-fit',
                   title_x = 0.5,
                  xaxis_title='x',
                  xaxis=dict(autorange=True, 
                             mirror=False,
                             ticks='outside',
                             showline=True,
                             linecolor='#000',
                             tickvals = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20],
                             tickfont = dict(size=11)),
                  yaxis_title='y',
                  yaxis=dict(autorange=True, 
                             mirror=False,
                             ticks='outside', 
                             showline=True,
                             linecolor='#000',
                             tickvals = [0, 1, 2, 3, 4, 5, 6, 7],
                             tickfont = dict(size=11)),
                 legend=dict(x=0.8, 
                              y=0.9,
                              bordercolor="Gray",
                              borderwidth=1),
                  plot_bgcolor='#fff',
                  hovermode = 'x',
                  width = 500,
                  height = 320,
                  font = dict(size = 10),
                  margin=go.layout.Margin(l=50,
                                         r=50,
                                         b=60,
                                         t=35))

plot(figa, filename = 'fig10_b.html', config = config)
display(HTML('fig10_b.html'))

Параметрите добиени: p1: 21.0326, p2: 6.1003, p3: 1.2009.

Поголемиот дел од curve-fitting што го користиме подоцна, беше околку оваа општа идеја. Следно е појаснување на моделот (код + транзиции на состојби) кој го користеме за моделриање на COVID-19.