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

library("EnvStats")


# Rovnice modelu
model<-function(t,N0,a,b) {
  return(N0*exp((a-b)*t))
}

# Definice parametrů

# věk dožití: 10-12 let
# porodnost: 3-4 krát ročně 1-7 mláďat (tj. 3-28 mláďat ročně)
N0<-100
a <-7/12
b <-1/(11*12)

# Graf
plot(c(0:12),model(c(0:12),N0,a,b),type="l",col="navyblue",lwd=3,main="Neomezený vývoj populace v čase",xlab="Čas [měsíc]",ylab="Počet jedinců [1]")

#########################################
# Intervalová aritmetika - pomocí cyklů #
#########################################

# Definice parametrů
a.interval<-c(0.125,0.6475,1.17)
b.interval<-c(0.007,0.00765,0.0083)

plot(0,0,cex=0,xlim=c(0,12),ylim=c(0,1e5),main="Neomezený vývoj populace v čase",xlab="Čas [měsíc]",ylab="Počet jedinců [1]")

for(a in a.interval) {
  for(b in b.interval) {
    lines(c(0:12),model(c(0:12),N0,a,b),col="navyblue",lwd=3)
  }
}


###########################################
# Pravděpodobnostní analýza - Monte Carlo #
###########################################

# Definice parametrů
a.rozdeleni<-rtri(1e5,0.125,1.17,0.6475)
b.rozdeleni<-rtri(1e5,0.007,0.0083,0.0076)

hist(a.rozdeleni)
hist(b.rozdeleni)

plot(0,0,cex=0,xlim=c(0,12),ylim=c(0,1e5),main="Neomezený vývoj populace v čase",xlab="Čas [měsíc]",ylab="Počet jedinců [1]")

for(i in 1:1e5) {
  a<-a.rozdeleni[i]
  b<-b.rozdeleni[i]
  lines(c(0:12),model(c(0:12),N0,a,b),col="navyblue",lwd=3)
}


# A nyní s průhlednými čarami
col2rgb("navyblue")

plot(0,0,cex=0,xlim=c(0,12),ylim=c(0,1e5),main="Neomezený vývoj populace v čase",xlab="Čas [měsíc]",ylab="Počet jedinců [1]")

for(i in 1:1e5) {
  a<-a.rozdeleni[i]
  b<-b.rozdeleni[i]
  lines(c(0:12),model(c(0:12),N0,a,b),col=rgb(0,0,128/255,0.01),lwd=3)
}

# Histogramy výsledků
reseni<-as.data.frame(matrix(NA,13,0))

for(i in 1:1e5) {
  a<-a.rozdeleni[i]
  b<-b.rozdeleni[i]
  reseni<-cbind(reseni,model(c(0:12),N0,a,b))
}

hist(as.numeric(reseni[2,]))
hist(as.numeric(reseni[5,]))


