animal.risk=function(n,N,p)
  
{
  par(mfrow = c(1,2))
  dist.binom=dbinom(0:N,size=N,prob=p) 
  successes=subset(dist.binom,dist.binom>0.0001) 
  
  plot( dist.binom,type ="h",xlab = "Number of infected animals", ylab = "Probability",
        main = "Distribuição do número de animais infectados",xlim=c(which(dist.binom==successes[1]),
                                                                     which(dist.binom==successes[length(successes)])))
  segments(x0=N*p,y0=0,x1=N*p,y1=dist.binom[(N*p)],col="red",lwd=2)
  
  conf.int1=qbinom(p=.025,size=N,prob=p)
  segments(x0=conf.int1,y0=0,x1=conf.int1,y1=dist.binom[(N*p)],col="blue",lwd=2 ) 
  
  conf.int2=qbinom(p=.975,size=N,prob=p)
  segments(x0=conf.int2,y0=0,x1=conf.int2,y1=dist.binom[(N*p)],col="blue",lwd=2 ) 
  
  cat("\n Most likely number of infected animals =",N*p,"\n Confidence interval of 95% =",
      conf.int1,conf.int2) 
  
  dist.hyper=dhyper(x=0:N,m=p*N,n=(1-p)*N,k=n) 
  successes2=subset(dist.hyper,dist.hyper>0.0001)
  
  plot(dist.hyper,type ="h",xlab = "Number of infected sampled animals", ylab = "Probability",
       main = "Distribuição do número de animais amostrados infectados",xlim=c(which(dist.hyper==successes2[1]),
                                                                               which(dist.hyper==successes2[length(successes2)])))
  segments(x0=n*p,y0=0,x1=n*p,y1=dist.hyper[(n*p)],col="red",lwd=2) 
  
  conf.int3=qhyper(p=.025,m=p*N,n=(1-p)*N,k=n)
  segments(x0=conf.int3,y0=0,x1=conf.int3,y1=dist.hyper[(n*p)],col="blue",lwd=2 ) 
  
  conf.int4=qhyper(p=.975,m=p*N,n=(1-p)*N,k=n)
  segments(x0=conf.int4,y0=0,x1=conf.int4,y1=dist.hyper[(n*p)],col="blue",lwd=2 ) 
  
  cat("\n Most likely number of infected sampled animals =",n*p,"\n Confidence interval of 95% =",
      conf.int3,conf.int4) 
}

{{:bie5782:01_curso_atual:alunos:trabalho_final:isabella.fontana:animal.risk_function.r|}}

{{:bie5782:01_curso_atual:alunos:trabalho_final:isabella.fontana:help_animal.risk_.r|}}