pocatecni.pocet<-10
a<-0.35
b<-0.25
max.cas<-10


###################################
#                                 #
#  Spojitý deterministický model  #
#                                 #
###################################

library("deSolve")

podminky<-c(N=pocatecni.pocet)
casy<-0:max.cas
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]-param[2]))))
}
parametry<-c(a,b)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="brown",lwd=3,which="N",main="Spojitý deterministický model",xlab="čas",ylab="počet jedinců")


##################################
#                                #
#  Diskértní stochastický model  #
#                                #
##################################

N<-pocatecni.pocet

for(cas in 1:max.cas) {
  N[cas+1]<-N[cas]
  for(jedinec in 1:N[cas]) {
    
    # rozmnozeni
    if(runif(1,0,1)<a) {
      N[cas+1]<-N[cas+1]+1
    }
    
    # umrti
    if(runif(1,0,1)<b) {
      N[cas+1]<-N[cas+1]-1
    }
  }
}

plot(0:max.cas,N,pch=19,col="brown",main="diskrétní stochastický model",xlab="čas",yalb="počet jedinců")
lines(0:max.cas,N,lwd=2,col="gold")
points(0:max.cas,N,pch=19,col="brown")


#########################################
#                                       #
#  Inverzní problém (pro stoch. model)  #
#                                       #
#########################################

nove.pocty<-N[-1]
predchozi.pocty<-N[-length(N)]

nove.pocty/predchozi.pocty-1

plot(nove.pocty~predchozi.pocty,pch=19,col="navyblue")
model<-lm(nove.pocty~predchozi.pocty+0)
summary(model)
model$coefficients[1]
confint(model,level=0.95)


########################
#                      #
#  Číslo podmíněnosti  #
#                      #
########################

N1<-1000
N2<-1000
r<-0.1

e<-0.05

for(rok in 1:max.cas) {
  N1[rok+1]<-N1[rok]
  for(jedinec in 1:N1[rok]) {
    reprodukujese<-runif(1,0,1)
    if(reprodukujese<r) {
      N1[rok+1]<-N1[rok+1]+1
    }
  }
}

for(rok in 1:max.cas) {
  N2[rok+1]<-N2[rok]
  for(jedinec in 1:N2[rok]) {
    reprodukujese<-runif(1,0,1)
    if(reprodukujese<(r+e)) {
      N2[rok+1]<-N2[rok+1]+1
    }
  }
}

plot(N1,main="Vývoj počtu jedinců v čase",xlab="Čas [rok]",ylab="Počet jedinců [1]",pch=19,col="brown",ylim=c(0,max(N2)))
points(1:11,N2,pch=19,col="chocolate")

max((abs(N2-N1)/N1)/(e/r))