#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8cproba <- function(rectangle, parameters,varange, dproba,qproba, vax=NULL, ci=-1, echelle=1, orientation="E", grafo=c("1","2","NA") ) #TITLE draws the continuous density of any proba within a plot #DESCRIPTION # draws the density of any continuous distribution # with \code{parameters} on the \code{varange} in any of the # four cardinal directions when the density (\code{dproba}) # and quantile (\code{qproba}) can be provided as arguments. #DETAILS # The plot is supposed already open. #KEYWORDS #PKEYWORDS #INPUTS #{rectangle} <> #{parameters} <> #{varange} <> #{dproba} <> #{qproba} <> #[INPUTS] #{vax} <> #{ci} <> #{echelle} <> #{orientation} <> #{grafo} <> #VALUE # Some density is added to the current plot , # and the maximum height of the curve is returned # (without taking into account the \code{echelle} coefficient). #EXAMPLE # plot(0,0,xlim=c(-1,6),ylim=c(-0.1,1.1),type="n"); # draw8histo(rectangle=c(0,5,0,1),3,5); #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE The use of echelle/kechelle is to be explicited... #AUTHOR J.-B. Denis #CREATED 13_03_19 #REVISED 13_03_21 #-------------------------------------------- { # the number of points to discretize the curve into a polygon nb <- 123; # checking check4tyle(rectangle, "numeric", 4, message = "'rectangle' must be a numeric(4)"); check4tyle(parameters,"nnumeric",-1, message = "'parameters' must be a named numeric vector"); check4tyle(varange,"numeric",2, message = "'range' must be a 2-lengthed vector"); check4tyle(orientation,"character",1,c("N","W","S","E"), message = "'orientation' is not valid"); check4tyle(grafo,"character",c(3,4), message = "'grafo' is not valid"); # decoding if (is.null(vax)) { vax <- numeric(0);} xl <- range(rectangle[1:2]); yl <- range(rectangle[3:4]); lty <- as.numeric(grafo[1]); lwd <- as.numeric(grafo[2]); if (is.na(grafo[3])) { grafo[3] <- "green";} if (grafo[3] == "NA") { grafo[3] <- "green";} if (is.na(grafo[4])) { grafo[4] <- "red";} if (grafo[4] == "NA") { grafo[4] <- "red";} coldens <- grafo[3]; coltail <- grafo[4]; # vx <- seq(varange[1],varange[2],length=nb)[-c(1,nb)]; credi <- ((ci > 0) & (ci < 1)); paraval <- paste(names(parameters)[1],"=",parameters[1],sep=""); for (ii in bd(2,length(parameters))) { paraval <- paste(paraval, paste(names(parameters)[ii],"=",parameters[ii],sep=""), sep=","); } afaire <- paste("vv <- dproba(vx,",paraval,")",sep=""); eval(parse(text=afaire)); if (credi) { afaire <- paste("cy <- qproba(c((1-ci)/2,(1+ci)/2),",paraval,")",sep=""); eval(parse(text=afaire)); quels <- which((vx>cy[1])&(vx>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8cprobas <- function(parameters,varange, dproba,qproba, pval=1:nrow(parameters), palab=NULL,sipalab=0.05, ci=-1, mti=NULL,sti=NULL, xlab=NULL,ylab=NULL, kechelle=1,grafo=c("1","2","green","red"), orientation="N", rplot=NULL, ...) #TITLE draws a continuous proba for different values of parameters #DESCRIPTION # draws the density of a continuous distribution # with different values of its parameters. #DETAILS # Just calling \code{draw8cproba} and assembling the results.\cr # All densities are drawn in an identical scale which is chosen # to fit into the intervals left by \code{pval}. To do so the upper # density is not taken into account since there is no upper bound # for it. #KEYWORDS #PKEYWORDS #INPUTS #{parameters} <> #{varange} <> #{dproba} <> #{qproba} <> #[INPUTS] #{pval} << y-values to place the distributions.>> #{palab} << Labels to associate to each distribution.>> #{sipalab} << length to give to \code{palab} in proportion of varange.>> #{ci} <> #{mti} <
> #{sti} <> #{xlab} << x-axis label for the plot.>> #{ylab} << y-axis label for the plot.>> #{kechelle} << Addotional coefficient to multiplicatively apply to the # \code{echelle} argument when calling the # \code{draw8cproba} function.>> #{grafo} <> #{orientation} << Either 'N' for the horizontal display or not for a # vertical display.>> #{rplot} << When \code{NULL} a new plot is open and the densities are drawn within it, # if not must be the definition of a recangle \code{c(xA,xB,yA,yB)} in # which the densities have to be drawn in an already opened plot.>> #{\dots} << to be given to the \code{plot} call.>> #VALUE # A complete plot is issued or an already plot is completed. # Returns \code{cbind(xlim,ylim)} in which the plot was created or modified. #EXAMPLE # draw8cprobas(as.data.frame(matrix(c(1,2,3,1,2,0.5),3,dimnames=list(NULL,c("mean","sd")))), # c(0,4),dnorm,qnorm) #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 13_03_19 #REVISED 13_03_22 #-------------------------------------------- { # some checking nbd <- nrow(parameters); if (nbd < 2) { stop("At least two distributions are required");} if (nbd != length(pval)) { stop("'parameters' and 'pval' are not consistent"); } if (!is.null(palab)) { if (nbd != length(palab)) { stop("'parameters' and 'palab' are not consistent"); }} # ordering the distribution according to pval ooo <- order(pval); parameters <- parameters[ooo,]; pval <- pval[ooo]; # computing the scale to draw the densities occupied <- 0.98; nbp <- 200; vax <- seq(varange[1],varange[2],length=nbp+2)[-c(1,nbp+2)]; deltay <- rep(0,nbd); for (dd in bc(nbd)) { paraval <- paste(names(parameters)[1],"=",parameters[dd,1],sep=""); for (ii in bd(2,length(parameters))) { paraval <- paste(paraval, paste(names(parameters)[ii],"=", parameters[dd,ii],sep=""), sep=","); } afaire <- paste("vay <- dproba(vax,",paraval,")",sep=""); eval(parse(text=afaire)); deltay[dd] <- max(vay); } echelle <- min(diff(pval)/deltay[-nbd]); pval <- c(pval,pval[nbd]+deltay[nbd]*echelle); # computing the window for each density rectangles <- matrix(NA,nbd,4); dimnames(rectangles) <- list(NULL,c("xm","XM","ym","YM")); for (ip in bc(nbd)) { rectangles[ip,1:2] <- varange; rectangles[ip,3:4] <- c(pval[ip],pval[ip+1]); } # computing the global window of the plot if (is.null(palab)) { sipalab <- 0;} xlim <- c(varange[1] - sipalab*diff(varange),varange[2]); ylim <- range(pval); if (is.null(rplot)) { # opening the plot if (is.null(xlab)) { xlab="";} if (is.null(ylab)) { ylab="";} if (is.null(mti)) { mti="";} if (is.null(sti)) { sti="";} nvarange <- varange; } else { # modifying the rectangle changer <- function(values,ancien,nouveau) { #print(rbind(ancien,nouveau)); # returns values from ancien to nouveau nouveau[1] + (values - ancien[1]) / diff(ancien) * diff(nouveau); } xlin <- range(rplot[1:2]); ylin <- range(rplot[3:4]); nvarange <- changer(varange, xlim,xlin); rectangles[,1:2] <- changer(rectangles[,1:2],xlim,xlin); rectangles[,3:4] <- changer(rectangles[,3:4],ylim,ylin); xlim <- xlin; ylim <- ylin; } if (orientation == "N") { ### horizontal disposal if (is.null(rplot)) { # opening the plot plot(xlim,ylim,type="n", yaxt="n", xlab=xlab, ylab=ylab, xlim=xlim,ylim=ylim, main=mti,sub=sti, ... ); } # writing the labels for (la in bf(palab)) { xpos <- 0.05; ypos <- 0.5; lax <- xlim[1] # lay <- mean(rectangles[la,3:4]); text(lax,lay,labels=palab[la],cex=1.5,pos=4); } # drawing each density for (ip in bc(nbd)) { mespara <- as.numeric(parameters[ip,]); names(mespara) <- names(parameters); draw8cproba(as.numeric(rectangles[ip,]), mespara, varange, dproba,qproba, vax=NULL, ci=ci, echelle=echelle/kechelle, orientation="N", grafo=grafo ); res <- cbind(xlim,ylim); } } else { ### vertical disposal if (is.null(rplot)) { # opening the plot plot(ylim,xlim,type="n", xaxt="n", xlab=xlab, ylab=ylab, xlim=ylim,ylim=xlim, main=mti,sub=sti, ... ); } # writing the labels for (la in bf(palab)) { #ypos <- 0.05; xpos <- 0.5; lay <- xlim[1] #+ xpos*(nvarange[1]-xlim[1]); lax <- mean(rectangles[la,3:4]); text(lax,lay,labels=palab[la],cex=1.5,pos=3); } # drawing each density for (ip in bc(nbd)) { mespara <- as.numeric(parameters[ip,]); names(mespara) <- names(parameters); recto <- as.numeric(rectangles[ip,c(3,4,1,2)]); draw8cproba(recto, mespara, varange, dproba,qproba, vax=NULL, ci=ci, echelle=echelle/kechelle, orientation="E", grafo=grafo ); res <- cbind(ylim,xlim); } } ## done and returning res; } #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< draw8nornor <- function(mean.x=0,sd.x=1,mean.y=function(x){x;},sd.y=1, xval=seq(mean.x-3*sd.x,mean.x+3*sd.x,length=5), xran=NULL,yran=NULL, style="I", ci=0.75, type="joint",marx=0.10,mary=0.10, colo=c("green","red"), spa=0.03, xlab=NULL,ylab=NULL,main=NULL, ... ) #TITLE draws a joint (or conditional) bi-normal distribution (not yet finished) #DESCRIPTION # draws the density of a couple of random variables (X,Y) following # a bi-normal distribution expressed by the marginal of X and # the conditional of Y|x.\cr # Still difficulty with type == "E". #DETAILS # Just calling \code{draw8cproba} and \code{draw8cprobas} with # some assembling.\cr # The difference between the 'joint' or 'conditional' is the # weighting applied to the non marginal distributions. #KEYWORDS #PKEYWORDS #INPUTS #[INPUTS] #{mean.x} << mean value for X>> #{sd.x} << standard deviation value for X.>> #{mean.y} << function with one argument representing an X # value returning the mean of Y conditioned by this value. # The argument can be vectorial as well.>> #{sd.y} << standard deviation value for Y.>> #{xval} <> #{xran} <> #{yran} <> #{style} <<"I" for isocontours, then the type is "joint" and "xval" not used, # "H" for horizontal cuttings and "V" for vertical cuttings.>> #{ci} << level of the credibility interval to draw (-1 means no)>> #{type} << either 'joint' for a joint plot, otherwise a conditional plot.>> #{marx} << if positive means that the X marginal must be drawn occupying # this proportion of xran.>> #{mary} << same for Y marginal.>> #{colo} << vector of the color for the density and the tails (when # a credibility interval is asked for).>> #{spa} << % of \code{xran} and \code{yran} to introduced between the joint # densities and the marginal ones.>> #{xlab} << X label>> #{ylab} << X label>> #{main} << main title>> #{\dots} << to be given to the \code{plot} call.>> #VALUE # A complete plot is issued #EXAMPLE # draw8nornor() #REFERENCE #SEE ALSO #CALLING #COMMENT #FUTURE #AUTHOR J.-B. Denis #CREATED 13_03_19 #REVISED 13_06_11 #-------------------------------------------- { # some checking to be done # some preparation margx <- c(mean=mean.x,sd=sd.x); margy <- c(mean=mean.y(mean.x), sd=sqrt((mean.y(mean.x+1)-mean.y(mean.x))^2*sd.x^2+sd.y^2)); joint <- data.frame(mean.y(xval),sd.y); names(joint) <- c("mean","sd"); if (is.null(xran)) { xran <- margx[1] + c(-3,3)*margx[2];} if (is.null(yran)) { margym <- mean.y(min(xval)); margyM <- mean.y(max(xval)); yran <- c(margym - 3*margy[2],margyM + 3*margy[2]); names(yran) <- NULL; } xli <- xran; yli <- yran; spax <- spay <- 0; if (mary > 0) { spax <- spa*diff(xran); xli[1] <- xli[1] - spax - mary * diff(xran);} if (marx > 0) { spay <- spa*diff(xran); yli[1] <- yli[1] - spay - marx * diff(yran);} jrect <- c(yran,xran); jrect <- c(xran,yran); jrecy <- c(xli[1], xran[1]-spax, yran); jrecx <- c( xran, yli[1],yran[1]-spax); # opening the graph if (is.null(xlab)) { xlab <- "";} if (is.null(ylab)) { ylab <- "";} if (is.null(main)) { main <- "";} plot(xli,yli,type="n",xlab=xlab,ylab=ylab, xlim=xli,ylim=yli, #yaxt="n",xaxt="n", main=main,...); # drawing the Y marginal if (mary > 0) { draw8cproba(jrecy,margy,yran,dnorm,qnorm, grafo=c(1,2,colo), vax=NULL,ci=ci,orientation="W"); } # drawing the X marginal if (mary > 0) { draw8cproba(jrecx,margx,xran,dnorm,qnorm, grafo=c(1,2,colo), vax=NULL,ci=-1,orientation="S"); } # drawing the joint fait <- FALSE; if (style == "H") { draw8cprobas(joint,xran,dnorm,qnorm,xval, grafo=c(1,2,colo), ci=ci, orientation="N",rplot=jrect); fait <- TRUE; } if (style == "V") { draw8cprobas(joint,xran,dnorm,qnorm,xval, grafo=c(1,2,colo), ci=ci, orientation="E",rplot=jrect); fait <- TRUE; } if (style == "I") { nbva <- 100; xva <- seq(xran[1],xran[2],length=nbva+2)[-c(1,nbva+2)]; yva <- seq(yran[1],yran[2],length=nbva+2)[-c(1,nbva+2)]; zva <- matrix(NA,nbva,nbva); sdy <- sqrt((mean.y(mean.x+1)-mean.y(mean.x))^2*sd.x^2+sd.y^2); for (ix in bf(xva)) { for (iy in bf(yva)) { xx <- xva[ix]; yy <- yva[iy]; muy <- mean.y(xx); zva[ix,iy] <- dnorm(xx,mean.x,sd.x) * dnorm(yy,muy,sdy); }} contour(xva,yva,zva,add=TRUE,lwd=2); fait <- TRUE; } if (!fait) { stop("'style' argument not understood."); } # returning invisible(); } #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>