# Regresné modely - metóda najmenších štvorcov
#
# Programové kódy: SageMath v.9.x(Python 3.0)
# Prerekvizity: štatistické metódy, Python, LaTex.
#
# Použitá literatúra:
# [1] R. E. Mezei, An introduction to SAGE programming, John Wiley,and Sons,2016 
# ISBN: 978-1-119-12280-7.
#
# [2] Gregory V.Bard. Sage for Undergraduates:Second Edition,Compatible with Python3. 
# American Mathematical Society,2022.ISBN978-1-4704-6155-3.
#  URL https://bookstore.ams.org/mbk-143/ - Sage pre pokročilých
#
# [3] Paul, Zimmermann and Alexandre, Casamayou and Nathann, Cohen and Guillaume Con- 
# nan and Thierry Dumont and Laurent Fousse and François Maltey and Matthias Meulien 
# and Marc Mezzarobba and Clément Pernet and Nicolas M. Thiéry and Erik Bray and 
# John Cremona and Marcelo Forets and Alexandru Ghitza and Hugh Thomas: 
# Computational Mathematics with SageMath, 2018, ISBN978-1-611975-45-1.
# https://my.siam.org/Store/Product/viewproduct/?ProductId=30174951,
#  https://www.sagemath.org/sagebook/english.html 
#
# E-kurz: P. Szabó: Numerická matematika, 2022
# Spúšťanie programov: 
#    1) Otvoriť lokalitu https://sagecell.sagemath.org/ 
#    2) Kopírovať obsah príkladu In[n] do SageMath boxu 
#    3) Kliknutie na "Evaluate"
# Príloha publikácie: 
# [4] Peter Szabó a Eva Baranová: Integrálny počet a numerická matematika, 2022, TUKE
# ISBN:978-80-553-4185-9
# Modifikované: 10.12.2022 

# Kreslenie množiny bodov (meraní)
points = [ (0, 1.1), (1, 2.1), (2, 5.1), (3, 9.9), (4, 14) ]
scatter_plot( points )

# Lineárna regresia - nájdenie regresnej priamky
points = [ (0, 1.1), (1, 2.1), (2, 5.1), (3, 9.9), (4, 14) ]
var("a,b")
model(x) = a*x+b
find_fit(points, model)

[a == 3.3600000000051473, b == -0.2800000000027889]

# Lineárna regresia - nájdenie a nakreslenie regresnej priamky
points = [ (0, 1), (1, 2),(2, 5), (3, 9), (4, 14) ]
var("x")
model(x) = a*x+b
#compute the range
xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
p  = scatter_plot( points )
fit =  find_fit(points, model)
f(x) = model(a =fit[0].rhs(), b = fit[1].rhs())
p += plot(f(x), xmin-1, xmax+1 )
show(p)
show(f(x)) 

$$\newcommand{\Bold}[1]{\mathbf{#1}}3.3000000000050163 \, x - 0.4000000000030508$$

# Lineárna regresia - nájdenie a nakreslenie regresnej priamky
# Vstupné body je možné modifikovať
var("a,b")
@interact
def LineFitInteract(
    points = input_box(default = [ (0, 1), (1, 2),\
                        (2, 5), (3, 9), (4, 14) ])):
    var("x")
    model(x) = a*x+b
    #compute the range
    xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
    xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
    p  = scatter_plot( points )
    fit =  find_fit(points, model)
    f(x) = model(a =fit[0].rhs(), b = fit[1].rhs())
    p += plot(f(x), xmin-1, xmax+1 )
    show(p)
    show(f(x)) 

# Lineárna regresia - nájdenie a nakreslenie regresnej priamky
# Vstupné body (100) sú náhodne generované, generátor náhodných čísel 
var("a,b")

points=[(x+random()/10, x+random()/10) for x in [0, 0.01,..,1]]
var("x")
model(x) = a*x+b
#compute the range
xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
p  = scatter_plot( points )
fit =  find_fit(points, model)
f(x) = model(a =fit[0].rhs(), b = fit[1].rhs())
p += plot(f(x), xmin, xmax )
show(p)
show(f(x)) 

$$\newcommand{\Bold}[1]{\mathbf{#1}}0.9804340699532846 \, x + 0.013995647295959479$$

# Regresia - nájdenie a nakreslenie regresnej funkcie
# Vstupné body je možné modifikovať
# Je možné modifikovať aj regresný model
var("a,b,x")
@interact
def LineFitInteract(
    points = input_box(default = [ (0, 1), (1, 2),\
                        (2, 5), (3, 9), (4, 14) ]),
    model = input_box( default = a*x+b )):
    var("x")
    model(x) = model
    #compute the range
    xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
    xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
    p  = scatter_plot( points )
    fit =  find_fit(points, model)
    f(x) = model(a =fit[0].rhs(), b = fit[1].rhs())
    p += plot(f(x), xmin-0.5, xmax+0.5 )
    show(p)
    show(f(x))

# Regresia - grafické porovnanie dvoch regresných modelov
# porovnanie kvadratického a goniometrického regresného modelu
var("a,b,c,d,x")
points = [ (0, 1), (1, 2),(2, 5), (3, 9), (4, 14) ]
model1 = a*x^2+b*x+c
model2 = a*cos(x)+b*sin(x)
    
var("x")
model1(x) = model1
model2(x) = model2
#compute the range
xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
#produces the scatter plot
p  = scatter_plot( points )
fit1 =  find_fit(points, model1)
if( len(model1) == 2 ):
    f1(x) = model1(a =fit1[0].rhs(), b = fit1[1].rhs())
elif( len(model1) == 3 ):
    f1(x) = model1(a =fit1[0].rhs(), b = fit1[1].rhs(),c = fit1[2].rhs())
elif( len(model1) == 4 ):
    f1(x) = model1(a =fit1[0].rhs(), b = fit1[1].rhs(),c = fit1[2].rhs(), d = fit1[3].rhs())
else:
       raise RuntimeError("model 1 should have 1-4 parameters!")
p += plot(f1(x), xmin, xmax, color = "red", legend_label = "model1 =$"+str (f1(x))+"$" )
fit2 =  find_fit(points, model2)
if( len(model2) == 2 ):
       f2(x) = model2(a =fit2[0].rhs(), b = fit2[1].rhs())
elif( len(model2) == 3 ):
       f2(x) = model2(a =fit2[0].rhs(), b = fit2[1].rhs(),c = fit2[2].rhs())
elif( len(model2) == 4 ):
       f2(x) = model2(a =fit2[0].rhs(), b = fit2[1].rhs(),c = fit2[2].rhs(), d = fit2[3].rhs())
else:
       raise RuntimeError("model 2 should have 1-4 parameters!")
p += plot(f2(x), xmin, xmax, color = "blue", legend_label = "model2 =$"+str (f2(x))+"$")
show(p)

var("a,b,c,d,x")
@interact
def LineFitInteract(
    points = input_box(default = [ (0, 1), (1, 2),\
                        (2, 5), (3, 9), (4, 14) ]),
 model1 = input_box( default = a*x^2+b*x+c ),
    model2 = input_box( default = a*cos(x)+b*sin(x) )):
    var("x")
    model1(x) = model1
    model2(x) = model2
    #compute the range
    xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
    xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
    #produces the scatter plot
    p  = scatter_plot( points )
    fit1 =  find_fit(points, model1)
    if( len(model1) == 2 ):
       f1(x) = model1(a =fit1[0].rhs(), b = fit1[1].rhs())
    elif( len(model1) == 3 ):
       f1(x) = model1(a =fit1[0].rhs(), b = fit1[1].rhs(), c = fit1[2].rhs())
    elif( len(model1) == 4 ):
       f1(x) = model1(a =fit1[0].rhs(), b = fit1[1].rhs(),c = fit1[2].rhs(), d = fit1[3].rhs())
    else:
       raise RuntimeError("model 1 should have 1-4 parameters!")
    p += plot(f1(x), xmin, xmax, color = "red", legend_label = "model1 =$"+str (f1(x))+"$" )
    fit2 =  find_fit(points, model2)
    if( len(model2) == 2 ):
       f2(x) = model2(a =fit2[0].rhs(), b = fit2[1].rhs())
    elif( len(model2) == 3 ):
       f2(x) = model2(a =fit2[0].rhs(), b = fit2[1].rhs(),c = fit2[2].rhs())
    elif( len(model2) == 4 ):
       f2(x) = model2(a =fit2[0].rhs(), b = fit2[1].rhs(),c = fit2[2].rhs(), d = fit2[3].rhs())
    else:
       raise RuntimeError("model 2 should have 1-4 parameters!")
    p += plot(f2(x), xmin, xmax, color = "blue", legend_label = "model2 =$"+str (f2(x))+"$")
show(p)

# Zdrojový kód programu v systéme SageMath
# Príklad z kurzu Aplikovaná matematika - 2019/2020 LS
# Program je spustiteľný aj na lokalite: https://sagecell.sagemath.org/
# Pre spustenie - zdrojový kód vložte do boxu Sagecell
# Matematický model šírenia vírusových ochorení
# Model neobsahuje presné výpočty, slúži len na odhad možného šírenia
# epidémie/pandémie. Tento základný model neobsahuje všetky parametre,
# ktoré majú vplyv na šírenie epidémie/pandémie.
# Našim všeobecným cieľom je znížiť "základné reprodukčné číslo" (zrc),
# keď zrc<1 tak pandémia (postupne) zaniká.
#
# 30.marca 2020, Peter Szabó
@interact
def sirenievirusu(
    zrc = slider(vmin=0.7, vmax=9, step_size=0.05,
               default=1.2, label="Základné reprodukčné číslo"),
    idb = input_box(default=6.4, label="Inkubačná doba v dňoch" ),
    pop = input_box(default=5400000, label="Populácia" ),
    cas = slider(vmin=1, vmax=45, step_size=1,
               default=10, label="Deň epidémie")):
    # funkcia šírenia pandémie, x je počet dní od vzniku pandémie
    f(x) = zrc^x/idb
    # Počet dní na infikovanie celej populácie
    dni=log(idb*pop,zrc).n()
    # Nakreslenie grafu šírenia pandémie
    a = 1
    b = cas+2
    if cas>2: a=cas-2

    q = line([(cas,0),(cas,2*f(cas))], color = "red")
    p = plot(f(x), (x, a, b), frame=true, gridlines='minor', axes_labels =[ "Počet dní","Počet infikovaných"])
    (p+q).show()
    print ("Veľkosť populácie                   :",pop)    
    print ("Deň epidémie                        :",cas)     
    if zrc > 1:
        if f(cas) >= pop:
            print ("Nakazená bola celá populácia        :",pop," za ",round(dni)," dní")
        else:
            print ("Počet infikovaných                  :",round(f(cas)))
            print ("Celá populácia môže byť infikovaná o:", round(dni) - cas, "dní")
    else:
        print ("Koniec epidémie")

# Zdrojový kód programu v systéme SageMath
# Príklad z kurzu Aplikovaná matematika - 2019/2020 LS
# Program je spustiteľný aj na lokalite: https://sagecell.sagemath.org/
# Pre spustenie - zdrojový kód vložte do boxu Sagecell
# Matematický model šírenia vírusových ochorení
# Model neobsahuje presné výpočty, slúži len na odhad možného šírenia
# epidémie/pandémie. Tento základný model neobsahuje všetky parametre,
# ktoré majú vplyv na šírenie epidémie/pandémie.
# Našim všeobecným cieľom je znížiť "základné reprodukčné číslo" (zrc),
# keď zrc<1 tak pandémia (postupne) zaniká.
#
# 30.marca 2020, Peter Szabó

zrc = 1.8 
# label="Základné reprodukčné číslo"
idb = 6.4 
# label="Inkubačná doba v dňoch"
pop = 5400000
# label="Populácia" ),
cas = 10
# label="Deň epidémie"

# funkcia šírenia pandémie, x je počet dní od vzniku pandémie
f(x) = zrc^x/idb
# Počet dní na infikovanie celej populácie
dni=log(idb*pop,zrc).n()
# Nakreslenie grafu šírenia pandémie
a = 1
b = cas+2
if cas>2: a=cas-2

q = line([(cas,0),(cas,2*f(cas))], color = "red")
p = plot(f(x), (x, a, b), frame=true, gridlines='minor', axes_labels =[ "Počet dní","Počet infikovaných"])
(p+q).show()
print ("Veľkosť populácie                   :",pop)    
print ("Deň epidémie                        :",cas)     
if zrc > 1:
    if f(cas) >= pop:
        print ("Nakazená bola celá populácia        :",pop," za ",round(dni)," dní")
    else:
        print ("Počet infikovaných                  :",round(f(cas)))
        print ("Celá populácia môže byť infikovaná o:", round(dni) - cas, "dní")
else:
        print ("Koniec epidémie")

Veľkosť populácie                   : 5400000
Deň epidémie                        : 10
Počet infikovaných                  : 56
Celá populácia môže byť infikovaná o: 20 dní

# Numerická a aplikovaná matematika - regresná analýza
# Vývoj priemernej hrubej mzdy na Slovensku podľa údajov štatistického úradu SR.
# Priemerné mzdy v r. 2022 (podľa Eurostat): Euro zóna: 2640 EUR - červená farba
# EU 2576 EUR - čierna farba, Európa: 2155 EUR (EU + ďalšie krajiny Albánsko, Ukrajina atď.) 
# Modrá (regresná) priamka - určuje kedy Slovensko pri súčastnom trende dosiahne úroveň 
# priemerných miezd EU z r. 2022. 
# Odvodové zaťaženie na Slovensku je 48.6% v EU priemer je cca. 38%

points = [(2016,912),(2017,954),(2018,1013),(2019,1092),(2020,1133),(2021,1211),(2022,1291)]
var("a, b, x")
model(x) = a*x+b
year=2045

#compute the range
xmin = min([points[i][0] for i in [0,1,..,len(points)-1]])
xmax = max([points[i][0] for i in [0,1,..,len(points)-1]])
p  = scatter_plot( points )
fit =  find_fit(points, model)
f(x) = model(a =fit[0].rhs(), b = fit[1].rhs())
p += plot(f(x), xmin-1, year,  frame=true, gridlines='minor', axes_labels =[ "Years","2022 average wages in EU"] )

heurope = 2155
heu = 2576
heu80 = 2061
heu60 = 1546
heuro = 2640

p += line([(2016,heurope),(year,heurope)], color = "green")
p += text("Europe 2022 = 2155 EUR", (2016,2200), color = "green", horizontal_alignment='left')
p += line([(2016,heu),(year,heu)], color = "black")
p += text("EU 2022 = 2576 EUR", (2016,2500), color = "black", horizontal_alignment='left')
p += line([(2016,heu80),(year,heu80)], color = "black")
p += text("EU80% 2022", (2016,1990), color = "black", horizontal_alignment='left')
p += line([(2016,heu60),(year,heu60)], color = "black")
p += text("EU60% 2022", (2016,1476), color = "black", horizontal_alignment='left')
p += line([(2016,heuro),(year,heuro)], color = "red")
p += text("Euro 2022 = 2640 EUR", (2016,2700), color = "red", horizontal_alignment='left')
p += text("Slovakia 2022 = 1291 EUR", (2016,1320), rotation=33.0, color = "blue", horizontal_alignment='left')

show(p)
show(f(x)) 

$$\newcommand{\Bold}[1]{\mathbf{#1}}63.24999999999951 \, x - 126615.17857142758$$