#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< margibibe <- function(n,a,b) #TITLE returns the marginal probabilities of a binomial * beta #DESCRIPTION # returns the vector of the \code{n+1} probabilities # of a discrete random variable defined as a binomial of # size \code{n} with a parameter \code{p} following # a beta distribution with parameters \code{a} and \code{b}. #DETAILS #KEYWORDS #PKEYWORDS #INPUTS #{n} <> #{a} <> #{b} <> #[INPUTS] #VALUE # a vector with \code{n+1} probabilities summing one, # associated to the values \code{0:n}. #EXAMPLE # margibibe(10,2,5); #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 11_09_20 #REVISED 11_09_20 #-------------------------------------------- { # checking check4tyle(n, "integer",1, c(0,Inf),message = "'n' is not valid"); check4tyle(a, "numeric",1, c(0,Inf),message = "'a' is not valid"); check4tyle(b, "numeric",1, c(0,Inf),message = "'b' is not valid"); # protection if(n+a+b > 120) { erreur(c(n,a,b),"Arguments with too high values"); } # special case if (n==0) { return(1);} # the general constant kk <- gamma(a+b)/gamma(a)/gamma(b)*gamma(n+1)/gamma(n+a+b); # the varying part vv <- 0:n; kv <- gamma(a+vv)/gamma(1+vv)*gamma(b+n-vv)/gamma(1+n-vv); # returning kk*kv; } #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8histo <- function(ou, hist, echelle=1, orientation="N", grafo=vector("list",0)) #TITLE draws an histogram within a plot #DESCRIPTION # draws the density of any discrete distribution # with rectangles (possibly colored) in any of the # four cardinal directions #DETAILS # The plot is supposed already open. #KEYWORDS #PKEYWORDS #INPUTS #{ou} <> #{hist} <> #[INPUTS] #{echelle} <> #{orientation} <> #{grafo} <> #VALUE # Nothing, some histogram is added to the current plot. #EXAMPLE # plot(0,0,xlim=c(-1,6),ylim=c(-0.1,1.1),type="n"); # draw8histo(ou=c(0,5,0,1),hist=list(n=1:5,f=1:5)); #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 11_09_20 #REVISED 11_09_20 #-------------------------------------------- { # checking check4tyle(ou, "numeric", 4, message = "'ou' must be a numeric(4)"); check4tyle(hist,"list",2, message = "'hist' must be a list(2)"); check4tyle(orientation,"character",1,c("N","W","S","E"), message = "'orientation' is not valid"); # decoding nb <- length(hist$n); xl <- range(ou[1:2]); yl <- range(ou[3:4]); if (is.null(grafo$lty)) { grafo$lty <- rep(1,nb);} if (is.null(grafo$lwd)) { grafo$lwd <- rep(1,nb);} if (is.null(grafo$wid)) { grafo$wid <- rep(1,nb);} if (is.null(grafo$col)) { grafo$col <- rep(NA,nb);} # just drawing the asked bars for (bar in bc(nb)) { # getting the assocated rectangle if (orientation=="N") { # abscissa xx <- rep(( 0.5 + (hist$n[bar]-min(hist$n))),4); xx <- xx + grafo$wid[bar]*c(-1,1,1,-1)/2; xx <- xl[1] + xx * diff(range(xl)) / (diff(range(hist$n))+1); # ordinates yy <- rep(yl[1],4); yy[3:4] <- yy[3:4] + hist$f[bar]/max(hist$f) * rep(diff(yl),2) * echelle; } if (orientation=="S") { # abscissa xx <- rep(( 0.5 + (hist$n[bar]-min(hist$n))),4); xx <- xx + grafo$wid[bar]*c(-1,1,1,-1)/2; xx <- xl[1] + xx * diff(range(xl)) / (diff(range(hist$n))+1); # ordinates yy <- rep(yl[2],4); yy[3:4] <- yy[3:4] - hist$f[bar]/max(hist$f) * rep(diff(yl),2) * echelle; } if (orientation %in% c("E","W")) { pause(paste("Sorry, this orientation (",orientation,") is not yet implemented")); } # drawing it polygon(xx,yy,lty=grafo$lty[bar], lwd=grafo$lwd[bar],col=grafo$col[bar] ); } # returning invisible(); } #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8beta <- function(ou, a,b, pp=c(0.25,0.50,0.75), ci=-1, echelle=1, orientation="E", grafo=c("1","2","NA"), ccol=c("red","green") ) #TITLE draws a beta density within a plot #DESCRIPTION # draws the density of a beta distribution # with parameters \code{a} and \code{b} in any of the # four cardinal directions #DETAILS # The plot is supposed already open. #KEYWORDS #PKEYWORDS #INPUTS #{ou} <> #{a} <> #{b} <> #[INPUTS] #{pp} <> #{ci} <> #{echelle} <> #{orientation} <> #{grafo} <> #{ccol} <> #VALUE # Nothing, some density is added to the current plot. #EXAMPLE # plot(0,0,xlim=c(-1,6),ylim=c(-0.1,1.1),type="n"); # draw8histo(ou=c(0,5,0,1),3,5); #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 11_09_20 #REVISED 11_09_21 #-------------------------------------------- { # constant nb <- 123; # checking check4tyle(ou, "numeric", 4, message = "'ou' must be a numeric(4)"); check4tyle(a,"numeric",1,c(0,Inf), message = "'a' must be a positive value"); check4tyle(b,"numeric",1,c(0,Inf), message = "'b' must be a positive value"); check4tyle(orientation,"character",1,c("N","W","S","E"), message = "'orientation' is not valid"); check4tyle(grafo,"character",3, message = "'grafo' is not valid"); # decoding if (is.null(pp)) { pp <- numeric(0);} xl <- range(ou[1:2]); yl <- range(ou[3:4]); lty <- as.numeric(grafo[1]); lwd <- as.numeric(grafo[2]); col <- grafo[3]; if (col=="NA") { col <- NA;} vx <- seq(0,1,length=nb)[-c(1,nb)]; credi <- ((ci > 0) & (ci < 1)); vv <- dbeta(vx,a,b); if (credi) { cy <- qbeta(c((1-ci)/2,(1+ci)/2),a,b); quels <- (vx>cy[1])&(vx>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8bino <- function(ou,n,p,ci=0.95,echelle=1,ccol=c("red","orange","green")) #TITLE draws a binomial #DESCRIPTION # draws the density of a binomial distribution # with vertical bars with the help of \code{draw8histo}. #DETAILS # The plot is supposed already open. #KEYWORDS #PKEYWORDS #INPUTS #{ou} <> #{n} <> #{p} <> #[INPUTS] #{ci} <> #{echelle} <> #{ccol} <> #VALUE # Nothing, some vertical bars are added to the current plot. #EXAMPLE # plot(0,0,xlim=c(-1,6),ylim=c(-0.1,1.1),type="n"); # draw8bino(ou=c(0,1,0,1),n=5,p=0.4); #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 11_09_19 #REVISED 13_06_11 #-------------------------------------------- { # preparing the call to draw8histo hist <- list(n=0:n,f=dbinom(0:n,n,p)); # if ((ci > 0) & (ci < 1)) { colo <- rep(3,n+1); alpha <- (1 - ci) / 2; dd <- dbinom(0:n,n,p); cu1 <- cumsum(dd); cu2 <- cumsum(dd[(n+1):1])[(n+1):1]; colo[(cu1 < alpha) | (cu2 < alpha)] <- 1; colo[sum(cu1 < alpha) + 1] <- 2; colo[n+1-sum(cu2 < alpha)] <- 2; } else { colo <- rep(NA,n+1); } form <- list(col=ccol[colo],wid=rep(1,n+1)); hist <- list(n=0:n,f=dbinom(0:n,n,p)); # drawing the histogram draw8histo(ou,hist,"N",form,echelle=echelle); # returning invisible(); } #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8binos <- function(n,pp,ci=-1, mti=NULL,sti=NULL, xlimi=NULL, main=NULL,sub=NULL, xlab=NULL,ylab=NULL, ccol=c("red","orange","green"), ...) #TITLE draws a binomial for different values of p #DESCRIPTION # draws the density of binomial distributions # with vertical bars #DETAILS #KEYWORDS #PKEYWORDS #INPUTS #{n} <> #{pp} <> #[INPUTS] #{ci} <> #{mti} <
> #{sti} <> #{xlimi} <> #{main} << Main title of the plot.>> #{sub} << Secondary title for the plot.>> #{xlab} << x-axis label for the plot.>> #{ylab} << y-axis label for the plot.>> #{ccol} <> #{\dots} << to be given to the \code{plot} call.>> #VALUE # A complete plot is issued #EXAMPLE # draw8binos(5,c(0.4,0.5,0.7)); #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 11_09_19 #REVISED 13_06_11 #-------------------------------------------- { ## separating constant ysep <- 0.05; ## values to be considered val <- 0:n; ## computing the range for abscissaes xlim <- c(-1,1) + range(val); ## computing the position for each histogram # height of each histo hm <- rep(NA,length(pp)); for (ip in bf(pp)) { hm[ip] <- max(dbinom(val,n,pp[ip])); } # starting position of each histo ou <- matrix(NA,length(pp),4); for (ip in bf(pp)) { ou[ip,1:2] <- xlim; ou[ip,3:4] <- c(sum(hm[bc(ip-1)]) + ysep * sum(hm) * (ip-1), sum(hm[bc(ip )]) + ysep * sum(hm) * (ip-1)); } ## computing the range for ordinates plus <- 1/5; ylim <- c(-ysep-plus,1+length(pp)*ysep)*sum(hm); ## preparing the labels lala <- unique(round(seq(0,n,le=7))); ## preparing the titles if (is.null(mti)) { mti <- paste(length(pp),"binomial distribution(s)");} if (is.null(sti)) { sti <- paste("size of",n);} ## opening the graph if (is.null(xlimi)) { xlimi <- xlim;} if (is.null(main)) { main <- mti;} if (is.null(sub)) { sub <- sti;} if (is.null(xlab)) { xlab <- "Binomial density";} if (is.null(ylab)) { ylab <- paste("For p =",paste(pp,collapse=" / "));} if (is.null(xlab)) { xlab <- "Binomial density";} if (is.null(ylab)) { ylab <- paste("For p =",paste(pp,collapse=" / "));} if (is.null(xlab)) { xlab <- "Binomial density";} if (is.null(ylab)) { ylab <- paste("For p =",paste(pp,collapse=" / "));} cat("xlim =",xlim,"ylim =",ylim,"\n"); plot(0,0,type="n", yaxt="n",xaxt="n", xlab=xlab, ylab=ylab, xlim=xlimi,ylim=ylim, main=main,sub=sub, ... ); ## writing the labels for (la in bf(lala)) { lal <- lala[la]; lax <- 0.5 + ou[1,1] + lal * diff(ou[1,1:2]) / (diff(range(val))+1); text(lax,ylim[1]/2,lal,cex=1.5,pos=1); } ## drawing each density for (ip in bf(pp)) { draw8bino(ou[ip,],n,pp[ip],ci,ccol=ccol); abline(h=ou[ip,2],lwd=1); } ## done invisible(); } #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8betabino <- function(n,a,b, pp=c(0.25,0.50,0.75), rr=NULL, ci=-1, type="joint", grafo="complete", mti=NULL,sti=NULL, xlab=NULL,ylab=NULL, ccol=c("red","orange","green"), ...) #TITLE draws a beta-binomial joint distribution #DESCRIPTION # draws on the same plot, the marginal of the beta for \code{p} # in ordinates, the conditional distributions for \code{r|p} # and the marginal of \code{r}. #DETAILS #KEYWORDS #PKEYWORDS #INPUTS #{n} <> #{a} <> #{b} <> #[INPUTS] #{pp} <> #{rr} <> #{ci} <> #{type} <> #{grafo} <> #{mti} <
> #{sti} <> #{xlab} <