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) 

}

animal.risk_function.r

help_animal.risk_.r