########################
#                      #
#  Analýza citlivosti  #
#                      #
########################

# Známé analytické řešení (spočítáno mimo R, např. v Pythonu, jak je ve skriptech)
N<-function(t,N0,r,K) {
  return(N0*K/(N0+(K-N0)*exp(-r*t)))
}

t<-0:20
t.s<-3
N0<-100
r<-7/12-1/(11*12)
K<-r/(1/(11*12))

plot(t,N(t,N0,r,K),
     type="l",
     lwd=3,
     col="brown",
     main="Modifikovaný model růstu",
     xlab="Čas",
     ylab="Velikost populace")


###########################
# Analyticky v čase t = 5 #
###########################

# Derivace známého analytického řešení podle N0 (spočítáno mimo R)
diff_N0<-function(t,N0,r,K) {
  K^2*exp(r*t)/((exp(2*r*t)-2*exp(r*t)+1)*N0^2+(2*K*exp(r*t)-2*K)*N0+K^2)
}
S_N0_a<-diff_N0(t.s,N0,r,K)            # index citlivosti pro N0 analyticky
S_N0_a.norm<-S_N0_a*N0/N(t.s,N0,r,K)   # normalizovaný index citlivosti pro N0 analyticky

# Derivace známého analytického řešení podle r (spočítáno mimo R)
diff_r<-function(t,N0,r,K) {
  -((K*N0^2-K^2*N0)*t*exp(r*t))/(N0^2*exp(2*r*t)+(2*K*N0-2*N0^2)*exp(r*t)+N0^2-2*K*N0+K^2)
}
S_r_a<-diff_r(t.s,N0,r,K)             # index citlivosti pro r analyticky
S_r_a.norm<-S_r_a*r/N(t.s,N0,r,K)     # normalizovaný index citlivosti pro r analyticky

# Derivace známého analytického řešení podle K (spočítáno mimo R)
diff_K<-function(t,N0,r,K) {
  (N0^2*exp(2*r*t)-N0^2*exp(r*t))/(K^2+(2*N0*exp(r*t)-2*N0)*K+N0^2*exp(2*r*t)-2*N0^2*exp(r*t)+N0^2)
}
S_K_a<-diff_K(t.s,N0,r,K)             # index citlivosti pro K analyticky
S_K_a.norm<-S_K_a*K/N(t.s,N0,r,K)     # normalizovaný index citlivosti pro K analyticky

barplot(abs(c("N0"=S_N0_a.norm,"r"=S_r_a.norm,"K"=S_K_a.norm)),
        main=paste0("Normalizované lokální koeficienty citlivosti 1. řádu v čase t = ",t.s," let"),
        xlab="Parametry modelu",
        ylab="Výpočet analyticky",
        col="cornflowerblue")


##################################
# Konečné diference v čase t = 5 #
##################################

podminky<-c("N"=N0)
podminky_N0.dolni<-c("N"=N0*0.9)      #N0 snížené o ΔN0 (tedy o 10 %)
podminky_N0.horni<-c("N"=N0*1.1)      #N0 zvýšené o ΔN0 (tedy o 10 %)
t<-0:20
soustava<-function(t,prom,par) {
  return(list(c(par["r"]*prom["N"]*(1-prom["N"]/par["K"]))))
}
parametry<-c("r"=r,"K"=K)
parametry_r.dolni<-c("r"=r*0.9,"K"=K) #r snížené o Δr (tedy o 10 %)
parametry_r.horni<-c("r"=r*1.1,"K"=K) #r snížené o Δr (tedy o 10 %)
parametry_K.dolni<-c("r"=r,"K"=K*0.9) #K snížené o ΔK (tedy o 10 %)
parametry_K.horni<-c("r"=r,"K"=K*1.1) #K snížené o ΔK (tedy o 10 %)
metoda<-"lsoda"

reseni<-as.data.frame(ode(podminky,t,soustava,parametry,metoda))                     # základní řešení beze změn

reseni_n0.dolni<-as.data.frame(ode(podminky_N0.dolni,t,soustava,parametry,metoda))   # řešení se změnou N0 o -ΔN0
reseni_n0.horni<-as.data.frame(ode(podminky_N0.horni,t,soustava,parametry,metoda))   # řešení se změnou N0 o +ΔN0
S_N0_n<-(reseni_n0.horni$N[t.s+1]-reseni_n0.dolni$N[t.s+1])/(N0*1.1-N0*0.9)          # index citlivosti pro N0 numericky
S_N0_n.norm<-S_N0_n*N0/reseni$N[t.s+1]                                               # normalizovaný index citlivosti pro N0 numericky

reseni_r.dolni<-as.data.frame(ode(podminky,t,soustava,parametry_r.dolni,metoda))     # řešení se změnou r o -Δr
reseni_r.horni<-as.data.frame(ode(podminky,t,soustava,parametry_r.horni,metoda))     # řešení se změnou r o +Δr
S_r_n<-(reseni_r.horni$N[t.s+1]-reseni_r.dolni$N[t.s+1])/(r*1.1-r*0.9)               # index citlivosti pro r numericky
S_r_n.norm<-S_r_n*r/reseni$N[t.s+1]                                                  # normalizovaný index citlivosti pro r numericky

reseni_K.dolni<-as.data.frame(ode(podminky,t,soustava,parametry_K.dolni,metoda))     # řešení se změnou K o -ΔK
reseni_K.horni<-as.data.frame(ode(podminky,t,soustava,parametry_K.horni,metoda))     # řešení se změnou K o +ΔK
S_K_n<-(reseni_K.horni$N[t.s+1]-reseni_K.dolni$N[t.s+1])/(K*1.1-K*0.9)               # index citlivosti pro K numericky
S_K_n.norm<-S_K_n*K/reseni$N[t.s+1]                                                  # normalizovaný index citlivosti pro K numericky

barplot(abs(c("N0"=S_N0_n.norm,"r"=S_r_n.norm,"K"=S_K_n.norm)),
        main=paste0("Normalizované lokální koeficienty citlivosti 1. řádu v čase t = ",t.s," let"),
        xlab="Parametry modelu",
        ylab="Výpočet numericky",
        col="brown2")


#################################
# Variance based v čase t = t.s #
#################################

# Stanovení rozptylu parametrů (kovariance jsou nulové) tak aby 90 % hodnot (5. až 95. percentil) leželo v průměru +- 10 %.
N0.mean<-N0
N0.sd<-(N0*1.1-N0*0.9)/(qnorm(0.95,0,1)-qnorm(0.05,0,1))

r.mean<-r
r.sd<-(r*1.1-r*0.9)/(qnorm(0.95,0,1)-qnorm(0.05,0,1))

K.mean<-K
K.sd<-(K*1.1-K*0.9)/(qnorm(0.95,0,1)-qnorm(0.05,0,1))


# Varianta 1 - parametry jsou nezávislé
S_N0_v<-S_N0_n^2*N0.sd^2   # index citlivosti N0 variance based
S_r_v <-S_r_n^2*r.sd^2     # index citlivosti r variance based
S_K_v <-S_K_n^2*K.sd^2     # index citlivosti K variance based

sum<-S_N0_v+S_r_v+S_K_v    # aby byl splněn aditivní předpoklad, že součet relativních variancí je 1, je třeba spočítat sumu

S_N0_v.rel<-S_N0_v/sum     # a podělit
S_r_v.rel <-S_r_v/sum      # jednotlivé
S_K_v.rel <-S_K_v/sum      # hodnoty

barplot(abs(c("N0"=S_N0_v.rel,"r"=S_r_v.rel,"K"=S_K_v.rel)),
        main=paste0("Relativní koeficienty citlivosti 1. řádu v čase t = ",t.s," let"),
        xlab="Parametry modelu",
        ylab="Výpočet variance based",
        col="gold")



# Nyní předpokládáme nenulovou kovarianci mezi r a K 
cov_r.K<--0.15   # kovariance -0.15 znamená, že čím je větší r, je menší K

parametry_r.dolni.K.dolni<-c("r"=r*0.9,"K"=K*0.9)   # r snížené o Δr (tedy o 10 %) a K snížené o ΔK (tedy o 10 %)
parametry_r.dolni.K.horni<-c("r"=r*0.9,"K"=K*1.1)   # r snížené o Δr (tedy o 10 %) a K zvýšené o ΔK (tedy o 10 %)
parametry_r.horni.K.dolni<-c("r"=r*1.1,"K"=K*0.9)   # r zvýšené o Δr (tedy o 10 %) a K snížené o ΔK (tedy o 10 %)
parametry_r.horni.K.horni<-c("r"=r*1.1,"K"=K*1.1)   # r zvýšené o Δr (tedy o 10 %) a K zvýšené o ΔK (tedy o 10 %)

reseni_r.dolni.K.dolni<-as.data.frame(ode(podminky,t,soustava,parametry_r.dolni.K.dolni,metoda))
reseni_r.dolni.K.horni<-as.data.frame(ode(podminky,t,soustava,parametry_r.dolni.K.horni,metoda))
reseni_r.horni.K.dolni<-as.data.frame(ode(podminky,t,soustava,parametry_r.horni.K.dolni,metoda))
reseni_r.horni.K.horni<-as.data.frame(ode(podminky,t,soustava,parametry_r.horni.K.horni,metoda))

S_r.K_n<-(reseni_r.dolni.K.dolni$N[t.s+1]-reseni_r.dolni.K.horni$N[t.s+1]-reseni_r.horni.K.dolni$N[t.s+1]+reseni_r.horni.K.horni$N[t.s+1])/((r*1.1-r*0.9)*(K*1.1-K*0.9))
S_r.K_v<-S_r.K_n^2*abs(cov_r.K)

sum<-S_N0_v+S_r_v+S_K_v+S_r.K_v

S_N0_v.rel <-S_N0_v/sum
S_r_v.rel  <-S_r_v/sum 
S_K_v.rel  <-S_K_v/sum
S_r.K_v.rel<-S_r.K_v/sum

barplot(abs(c("N0"=S_N0_v.rel,"r"=S_r_v.rel,"K"=S_K_v.rel,"r × K"=S_r.K_v.rel)),
        main=paste0("Relativní koeficienty citlivosti 1. a 2. řádu v čase t = ",t.s," let"),
        xlab="Parametry modelu",
        ylab="Výpočet variance based",
        col="forestgreen")


##############################
# Monte-Carlo v čase t = t.s #
##############################

# Již máme definované rozptyly, rocnou uvedeme příklad pro nezávislé parametry

rep<-10000

N0.rand<-rnorm(rep,N0.mean,N0.sd)
r.rand <-rnorm(rep,r.mean,r.sd)
K.rand <-rnorm(rep,K.mean,K.sd)

#N0
reseni.rand.N0<-numeric(0)
for(i in 1:rep) {
  podminky<-c("N"=N0.rand[i])
  t<-0:20
  soustava<-function(t,prom,par) {
    return(list(c(par["r"]*prom["N"]*(1-prom["N"]/par["K"]))))
  }
  parametry<-c("r"=r,"K"=K)
  metoda<-"lsoda"
  
  reseni<-as.data.frame(ode(podminky,t,soustava,parametry,metoda))
  reseni.rand.N0<-c(reseni.rand.N0,reseni$N[t.s+1])
}

S_N0_mc<-(quantile(reseni.rand.N0,0.95)-quantile(reseni.rand.N0,0.05))/(quantile(N0.rand,0.95)-quantile(N0.rand,0.05))
S_N0_mc.norm<-unname(S_N0_mc*N0/reseni$N[t.s+1])

#r
reseni.rand.r<-numeric(0)
for(i in 1:rep) {
  podminky<-c("N"=N0)
  t<-0:20
  soustava<-function(t,prom,par) {
    return(list(c(par["r"]*prom["N"]*(1-prom["N"]/par["K"]))))
  }
  parametry<-c("r"=r.rand[i],"K"=K)
  metoda<-"lsoda"
  
  reseni<-as.data.frame(ode(podminky,t,soustava,parametry,metoda))
  reseni.rand.r<-c(reseni.rand.r,reseni$N[t.s+1])
}

S_r_mc<-(quantile(reseni.rand.r,0.95)-quantile(reseni.rand.r,0.05))/(quantile(r.rand,0.95)-quantile(r.rand,0.05))
S_r_mc.norm<-unname(S_r_mc*r/reseni$N[t.s+1])

#K
reseni.rand.K<-numeric(0)
for(i in 1:rep) {
  podminky<-c("N"=N0)
  t<-0:20
  soustava<-function(t,prom,par) {
    return(list(c(par["r"]*prom["N"]*(1-prom["N"]/par["K"]))))
  }
  parametry<-c("r"=r,"K"=K.rand[i])
  metoda<-"lsoda"
  
  reseni<-as.data.frame(ode(podminky,t,soustava,parametry,metoda))
  reseni.rand.K<-c(reseni.rand.K,reseni$N[t.s+1])
}

S_K_mc<-(quantile(reseni.rand.K,0.95)-quantile(reseni.rand.K,0.05))/(quantile(K.rand,0.95)-quantile(K.rand,0.05))
S_K_mc.norm<-unname(S_K_mc*K/reseni$N[t.s+1])

# Grafy
par(mfrow=c(2,2))
hist(reseni.rand.N0,
     main=paste0("Řešení Monte-Carlo pro N0 v čase t = ",t.s," let"),
     xlab="Velikost populace",
     ylab="Četnost",
     col="pink")
hist(reseni.rand.r,
     main=paste0("Řešení Monte-Carlo pro r v čase t = ",t.s," let"),
     xlab="Velikost populace",
     ylab="Četnost",
     col="pink")
hist(reseni.rand.K,
     main=paste0("Řešení Monte-Carlo pro K v čase t = ",t.s," let"),
     xlab="Velikost populace",
     ylab="Četnost",
     col="pink")
par(mfrow=c(1,1))


barplot(abs(c("N0"=S_N0_mc.norm,"r"=S_r_mc.norm,"K"=S_K_mc.norm)),
        main=paste0("Normalizované lokální koeficienty citlivosti Monte Carlo v čase t = ",t.s," let"),
        xlab="Parametry modelu",
        ylab="Výpočet variance based",
        col="orange")


#########################
#                       #
#  Analýza neurčitosti  #
#                       #
#########################

# Přes Monte-Carlo metody pěkně přejdeme k analýze neurčitosti
reseni.rand<-numeric(0)
for(i in 1:rep) {
  podminky<-c("N"=N0.rand[i])
  t<-0:20
  soustava<-function(t,prom,par) {
    return(list(c(par["r"]*prom["N"]*(1-prom["N"]/par["K"]))))
  }
  parametry<-c("r"=r.rand[i],"K"=K.rand[i])
  metoda<-"lsoda"
  
  reseni<-as.data.frame(ode(podminky,t,soustava,parametry,metoda))
  reseni.rand<-c(reseni.rand,reseni$N[t.s+1])
}

hist(reseni.rand,
     main=paste0("Řešení  Monte-Carlo pro nezávislé parametry v čase t = ",t.s," let"),
     xlab="Velikost populace",
     ylab="Četnost",
     col="pink")
median(reseni.rand)


# Nyní zohledníme kovarianci mezi r a K, tak jak jsme si ji nadefinovali dříve
library("MASS")                # knihovna pro vícedimenzionální pravděpodobnostní rozdělení
cor_r.K<-cov_r.K/(r.sd*K.sd)   # jen pro zajímavost, korelace odpovídající zadané kovarianci

N0.rand<-rnorm(rep,N0.mean,N0.sd)
r.K.rand <-mvrnorm(rep,c(r.mean,K.mean),matrix(c(r.sd^2,cov_r.K,cov_r.K,K.sd^2),2,2))   #mvrnorm vygeneruje dvojrozměrné normální rozdělení se zadanou kovariancí proměnných
r.rand<-r.K.rand[,1]
K.rand<-r.K.rand[,2]
plot(r.rand,K.rand,pch=19)

reseni.rand<-numeric(0)
for(i in 1:rep) {
  podminky<-c("N"=N0.rand[i])
  t<-0:20
  soustava<-function(t,prom,par) {
    return(list(c(par["r"]*prom["N"]*(1-prom["N"]/par["K"]))))
  }
  parametry<-c("r"=r.rand[i],"K"=K.rand[i])
  metoda<-"lsoda"
  
  reseni<-as.data.frame(ode(podminky,t,soustava,parametry,metoda))
  reseni.rand<-c(reseni.rand,reseni$N[t.s+1])
}

hist(reseni.rand,
     main=paste0("Řešení  Monte-Carlo pro nezávislé parametry v čase t = ",t.s," let"),
     xlab="Velikost populace",
     ylab="Četnost",
     col="pink")
median(reseni.rand)


