Есть проблема с моим кодом, когда выходит графика, и я не могу понять, где ошибки?
simulation_epidemie_SIR<-function(taille1,suspect1,infect1,remov1,alpha1,beta1)
{
taille<-taille1 #population
suscept<-rep(0,taille) #suspect
infect<-rep(0,taille) #infected persons
remov<-rep(0,taille) #non infected
time<-rep(0,taille) #time
nvcas<-rep(0,taille) #new cases
suscept[1]<-suspect1 #my nombre of suspect
infect[1]<-infect1
remov[1]<-remov1
beta<-beta1 #infectious rate/day
alpha<-alpha1 #recovery rate/day
for(i in 1:(taille-1)) #from here it's my teacher code
{
if(infect[i]!=0) # i did'nt change any thing here
{
s<--log(runif(1))/(beta*suscept[i]*infect[i]+alpha*
infect[i])
time[i+1]<-time[i]+s
rand<-runif(1) #
proba<-beta*suscept[i]*infect[i]/(beta*suscept[i]*infect[i]+alpha*infect[i])
if (rand<=proba)
{
suscept[i+1]<-suscept[i]-1 #
infect[i+1]<-infect[i]+1 #
remov[i+1]<-remov[i]
nvcas[i+1]<-1 #
}
if (rand>proba)
{
suscept[i+1]<-suscept[i] #the method of my teacher i did'nt change any thing here
infect[i+1]<-infect[i]-1
remov[i+1]<-remov[i]+1 #
nvcas[i+1]<-0 #
}
} ##
} #
suscept_red<-suscept[suspect=!0] #from here i write it
infect_red<-infect[suscept=!0] # for starting
remov_red<-remov[suscept=!0] # i could do better here
time_red<-time[suspect=!0] # i think here ther eis a proplem
par(mfrow=c(2,3)) #with this 4 codes
plot(time_red,suscept_red,type = "l",col= "red",xlab = "Temps",ylab = "Nb de cas susceptible",main = "Evolution des patients suspects") #i have a problem with this graphic
plot(time_red,infect_red,type = "l",col= "blue",xlab = "Temps",ylab = "Nb de cas infectés",main = "Evolution des patients infectés") #and this one
plot(time_red,remov_red,type = "l",col= "green",xlab = "Temps",ylab = "Nb de cas retirés",main = "Evolution des patients suspects") #also this one
tnvcas<-rep(0,365) #for one year from day one
for (i in 1:365)
tnvcas[i]<-sum(nvcas[(i-1)<time & time<=i])
plot(1:365,tnvcas,type = "l",col = "purple",xlab = "Jours", ylab = "Nb de nouveau cas"
,main = "Evolution du Nb de nouveau cas par jours", font.main = 1.5,font.lab = 3
,font.sub = 3,cex.main = 1.8, cex.lab = 1.5, cex.sub = 1.5) #plot of new cases
sum(nvcas) #to sum new cases
sum(tnvcas)
incidence<- tnvcas/taille1 #here to calculte the impact = new cases/ population
incidence_red<-incidence[incidence=!0] #here to calculate the decrease in incidence
return(list(taille=sum(tnvcas),incidence=incidence,durée=length(incidence_red))) #
#
}
simulation_epidemie_SIR2<-simulation_epidemie_SIR(100000,10000,5,90000,0.4,0.0000524) #
#these are the chiffres i have to use them
