inst/original/rbsb4.code.r

 
###########################################
###########################################
########
#((((((( NEW S4 CLASS empirical
########
###########################################
###########################################

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
valid8empirical <- function(object)
#TITLE  checks an /empirical/
#DESCRIPTION
#   This function checks an /empirical/ object
#DETAILS
# It is the validity method for /empirical/ objects.
#KEYWORDS classes
#INPUTS
#{object} <<The empirical object to be validated.>>
#[INPUTS]
#VALUE
# TRUE when the object seems acceptable
# else a character describing the error(s)
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# valid8empirical(rbsb.epi1);
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_13
#--------------------------------------------
{
    if (class(object)!="empirical") {
        erreur(NULL,paste("This object is not of class 'empirical' but '",class(object),"'",sep=""));
    }
    #
    rrr <- valid8des(object@des);
    if (identical(rrr,TRUE)) {
      res <- character(0);
    } else {
      res <- rrr;
    }
    #
    res <- c(res,check4tyle(object@repa,"numeric",2,c(0,Inf),message="@repa not accepted"));
    if (sum(object@repa) <= 0) {
      res <- c(res,"@repa cannot be c(0,0), a 'true' probability is wanted");
    }
    #
    res <- c(res,check4tyle(object@sup,"numeric",2,message="@repa not accepted"));
    if (diff(object@sup)<0) {
      res <- c(res,"The @sup must defined an interval");
    }
    #
    res <- c(res,check4tyle(object@xw,"matrix",-1,message="@xw not accepted"));
    if (is.matrix(object@xw)) {
      if (object@repa[1] > 0) {
        if (nrow(object@xw) == 0) {
          res <- c(res,"When @repa[1]>0, @xw cannot be empty");
        }
      }
      if (ncol(object@xw)!=2) {
        res <- c(res,"@xw must have two columns");
      }
      if (any(object@xw[,2]<0)) {
        res <- c(res,"@xw[,2] must be non negative weights");
      }
    }
    if (any(is.na(object@xw)) |
        any(is.infinite(object@xw))) {
      res <- c(res,"All values of @xw must be defined");
    }
    #
    res <- c(res,check4tyle(object@xy,"matrix",-1,message="@xy not accepted"));
    if (is.matrix(object@xy)) {
      if (object@repa[2] > 0) {
        if (nrow(object@xy) < 3) {
          res <- c(res,"When @repa[2]>0, @xy must have at least 3 columns");
        }
      }
      if (ncol(object@xy)!=2) {
        res <- c(res,"@xy must have two columns");
      }
      if (any(diff(object@xy[,1])<0)) {
        res <- c(res,"@xy[,1] must be non decreasing");
      }
      if (any(object@xy[,2]<0)) {
        res <- c(res,"@xy[,2] must be non negative weights");
      }
      if (nrow(object@xy)>1) { if(any(object@xy[c(1,nrow(object@xy)),2]!=0)) {
        res <- c(res,"First and Last element of @xy must be null");
      }}
    }
    if (any(is.na(object@xy)) |
        any(is.infinite(object@xy))) {
      res <- c(res,"All values of @xy must be defined");
    }
    #
    if (length(res)== 0) { res <- TRUE;
    } else { erreur(res,w=rbsb.mwa);}
    res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>


###########################################
# description
setClass("empirical", representation(
    des ="des",        # description of the object
    sup ="numeric",    # the support for the distribution
                       # when some values of 'xw' and 'xy'
                       # are outside the support (the bounds
                       #  are included), they are eliminated.
    repa="numeric",    # repartition between the discrete [1]
                       #    and continuous [2] parts
                       #    irrespectively of the support, this
                       #    means that the repartition is applied
                       #    before truncating.
    xw="matrix",     # matrix with two columns (x,p) describing
                       #    the discrete density
    xy="matrix"      # matrix of two columns (x,y) describing
                       #    the continuous density
                         ),
                      prototype=list(
    des=new("des",name="empi"),
    sup=c(0,10),
    repa=c(0.3,0.7),
    xw=matrix(c(1:3,rep(1,3)),3),
    xy=matrix(c(0,0,1,0,1,0),3)
                         ),
                      validity=valid8empirical
        );
#

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
print8empirical <- function(x,...,what="rswc",empha=0)
#TITLE  prints an /empirical/
#DESCRIPTION
#   This function prints in an /empirical/ object.
#DETAILS
# Global constant \code{rbsb.mwi} is used to justify the paragraphes.
# It is the specific print for /empirical/ objects.
#KEYWORDS classes
#INPUTS
#{x} <<The \code{empirical} object to print.>>
#[INPUTS]
#{\dots} <<Further arguments to be passed to the print function.>>
#{what} <<the fields to print:
#          \code{a}=all fields,
#          \code{d}=description,
#          \code{s}=support,
#          \code{r}=repartition,
#          \code{w}=discrete part,
#          \code{c}=continuous part.>>
#{empha} <<Emphasize level of printing;
#           between 0 and 10>>
#VALUE
# nothing but a print is performed
#EXAMPLE
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
# For the moment, only empha==0 is done
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_10
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  che <- valid8empirical(x);
  if (!identical(che,TRUE)) {
      erreur(che,"/empirical/ is not valid");
  }
}
empha <- round(max(0,min(10,empha)));
#
empha <- 0;
what <- tolower(what);
if (expr3present("a",what)) { what <- "dsrwc";}
# preparing the constant according to the emphasis
if (empha == 0) {
}
dimnames(x@xw) <- list(NULL,c("X","W"));
dimnames(x@xy) <- list(NULL,c("X","Y"));
#
# printing
if (expr3present("d",what)) {
  print8des(x@des,empha=empha);
}
#
if (expr3present("s",what)) {
  form3titre("Support:",empha);
  form3line(0,wid=1);
  cat("       From",x@sup[1],"to",x@sup[2],"\n");
}
#
if (expr3present("r",what)) {
  form3titre("Repartition:",empha);
  form3line(0,wid=1);
  cat("      ",x@repa[1],"for discrete and",x@repa[2],"for continuous\n");
}
#
if (expr3present("w",what)) {
  if (x@repa[1] > 0) {
    form3titre("The discrete part:",empha);
    form3line(0,wid=1);
    print(x@xw);
  } else {
    form3titre("No discrete part",empha);
    form3line(0,wid=1);
  }
}
#
if (expr3present("c",what)) {
  if (x@repa[2] > 0) {
    form3titre("The continuous part:",empha);
    form3line(0,wid=1);
    print(x@xy);
  } else {
    form3titre("No continuous part",empha);
    form3line(0,wid=1);
  }
}
# returning
invisible();
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
plot8empirical <- function(x,y="useless",...,
                           what="dc",
                           ticks="s",
                           xli=NULL,ylid=NULL,ylic=NULL,main=NULL)
#TITLE  plots an /empirical/
#DESCRIPTION
#   This function plots an /empirical/ object.
#DETAILS
#KEYWORDS classes
#INPUTS
#{x} <<The \code{empirical} object to plot.>>
#[INPUTS]
#{y} << For compatibility with the generic plot function. Useless here.>
#{\dots} <<Further arguments to be passed to the plot function.>>
#{what} <<What to plot? \code{d} for density and \code{c} for cumulative.>>
#{ticks} <<Markers of the discontinuity of the function. When numerical,
#          0 means no ticks if not the height of the ticks.
#          When character, \code{s} means until the curve,
#          \code{S} the same but ticks will have a minimum length, if not nothing.>>
#{xli} <<possible imposed \code{xlim} for the density and the cumulative>>
#{ylid} <<possible imposed \code{ylim} for the density>>
#{ylic} <<possible imposed \code{ylim} for the cumulative>>
#{main} <<Some prepared title for the plot.>>
#VALUE
# nothing but a plot is performed
#EXAMPLE
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_21
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  che <- valid8empirical(x);
  if (!identical(che,TRUE)) {
      erreur(che,"/empirical/ is not valid");
  }
}
# normalizing
x <- normalize8empirical(x);
# looking for the frame
if (is.null(xli)) {
  xli <- numeric(0);
  if (x@repa[1]>0) { xli <- c(xli,x@xw[,1]);}
  if (x@repa[2]>0) { xli <- c(xli,x@xy[,1]);}
  if (length(xli)==0) { return();}
  xli <- geom3lmargin(xli);
} else {
  check4tyle(xli,"numeric",2,message="Bad 'xli'");
}
if (is.null(ylid)) {
  ylid <- numeric(0);
  if (x@repa[1]>0) { ylid <- c(ylid,x@xw[,2]);}
  if (x@repa[2]>0) { ylid <- c(ylid,x@xy[,2]);}
  ylid <- range(0,ylid);
} else {
  check4tyle(ylid,"numeric",2,message="Bad 'ylid'");
}
if (is.null(ylic)) {
  ylic <- 0:1;
} else {
  check4tyle(ylic,"numeric",2,message="Bad 'ylic'");
}
if (is.null(main)) {
  titd <- paste(x@des@name,"(density)");
  titc <- paste(x@des@name,"(cumulative)");
} else {
  titd <- main;
  titc <- main;
}
#
#
### THE DENSITY
#
if (expr3present("d",what)) {
  plot(0,0,
       xlab="values",ylab="density",
       xlim=xli,ylim=ylid,main=titd,type="n");
  abline(h=0);
  #
  # drawing the continuous part
  if (x@repa[2] > 0) {
    NN <- 1000;
    xxx <- seq(xli[1],xli[2],length=NN);
    yyy <- continuous2density(x@xy,xxx);
    polygon(xxx,yyy,density=NA,col="lightgrey")
    lines(xxx,yyy);
  }                        
  #
  # drawing the discrete part
  if (x@repa[1] > 0) {
    for (ww in bc(nrow(x@xw))) {
      xx <- x@xw[ww,1];
      yy <- x@xw[ww,2]
      lines(c(xx,xx),c(0,yy),lwd=5);
    }
  }
}
#
#
### THE CUMULATIVE
#
if (expr3present("c",what)) {
  plot(0,0,
       xlab="values",ylab="probability",
       xlim=xli,ylim=ylic,main=titc,type="n");
  abline(h=c(0,1));
  #
  # computing the cumulative
  N <- 1000;
  ran <- range5empirical(x)$remp;
  if (diff(ran)==0) {
    # to tackle properly Dirac distributions
    ran <- (xli+ran)/2;
  }
  xxx <- seq(ran[1],ran[2],length=N);
  cumu <- cbind(xxx,pempirical(xxx,x));
  # adding a starting point when the discrete part
  # has got the lowest value
  if (cumu[1,2] > 0) {
    cumu <- rbind(c(cumu[1,1],0),cumu);
  }
  #
  # drawing it
  lines(cumu,lwd=2);
  #
  # adding the ticks
  xt <- unique(c(x@xw[,1],x@xy[,1]));
  if (is.numeric(ticks)) {
    if (ticks > 0) {
      for (xx in xt) {
        lines(rep(xx,2),c(0,min(1,ticks)),lty=3);
      }
    }   
  }
  if (is.character(ticks)) {
    if (tolower(ticks)=="s") {
      for (xx in xt) {
        ix <- which(min(abs(xx-cumu[,1]))==
                    abs(xx-cumu[,1])
                   )[1];
        xxx <- cumu[ix,1]; yyy <- cumu[ix,2];
        if (ticks=="S") { yyy <- max(0.05,yyy);}
        lines(rep(xxx,2),c(0,yyy),lty=3);
      }
    }   
  }
}
# returning
invisible();
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#
#

setMethod("print",signature(x = "empirical"), print8empirical);

setMethod( "plot",signature(x = "empirical"), plot8empirical);

#
#
#

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
normalize8empirical <- function(empi)
#TITLE  normalizes an /empirical/
#DESCRIPTION
#   This function normalizes an /empirical/ object,
# that is returns the object with weight transformed  
# into probabilities (taking care of the two parts)
#DETAILS
# Also a sorting of the discrete values is performed
# and multiple identical values are merged. Zero values
# are eliminated as much as possible. Also are values outside
# the proposed support.
#KEYWORDS classes
#INPUTS
#{empi} <<The \code{empirical} object to normalize.>>
#[INPUTS]
#VALUE
# the transformed /empirical/
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# print(normalize8empirical(new("empirical")));
# e1 <- new("empirical",des=char2des("e1"),sup=c(0,10),repa=c(10,0),xw=matrix(c(1,2,1,2,1,1:5),ncol=2));
# print(e1);
# print(normalize8empirical(e1));
# e2 <- new("empirical",des=char2des("e2"),sup=c(2.5,3.5),repa=c(0,5),xy=matrix(c(1:5,0,1,1,1,0),ncol=2));
# print(e2);
# print(normalize8empirical(e2));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_14
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
    erreur(che,"/empirical/ is not valid");
  }
}
res <- empi;
# the repartition
res@repa <- res@repa / sum(res@repa);
#
## the discrete part
if (res@repa[1] == 0) {
  res@xw <- matrix(0,0,2);
} else {
  # eliminating the zero y values
  ok <- which(res@xw[,2]>0);
  res@xw <- res@xw[ok,,drop=FALSE];
  # putting things in order
  oo <- order(res@xw[,1]);
  res@xw <- res@xw[oo,,drop=FALSE];
  # merging identical x values
  dd <- diff(res@xw[,1]);
  for (gg in rev(bf(dd))) {
    if (dd[gg] == 0) {
      # gg and gg+1 are identical
      ## 1: merging
      res@xw[gg,2] <- res@xw[gg,2] + res@xw[gg+1,2];
      ## 2: removing
      res@xw <- res@xw[-(gg+1),,drop=FALSE];
    }
  }
  res@xw[,2] <- res@xw[,2]/sum(res@xw[,2]) * res@repa[1];
}
#
## the continuous part
if (res@repa[2] == 0) {
  res@xy <- matrix(0,0,2);
} else {
  # computing first the sum
  my <- (res@xy[-1,2] + res@xy[-nrow(res@xy),2]) / 2;
  susu <- sum(diff(res@xy[,1]) * my);
  res@xy[,2] <- res@xy[,2]/susu * res@repa[2];
}
#
## removing values outside the support
ra <- range(res@xw[,1],res@xy[,1]);
eli <- (ra[1] < res@sup[1]) | (res@sup[2] < ra[2]);
if (eli) {
  if (res@repa[1]>0) {
    # discrete part
    ok <- which((res@sup[1] <= res@xw[,1]) &
                (res@xw[,1] <= res@sup[2]));
    res@xw <- res@xw[ok,,drop=FALSE];
  }
  if (res@repa[2]>0) {
    # continuous part
    rem <- numeric(0);
    # continuous part on the lower side
    cb <- sum(res@xy[,1] < res@sup[1]);
    if (cb > 0) {
      # something has to be removed
      if (cb == nrow(res@xy)) {
        # everything
        res@xy <- matrix(0,0,2);
      } else {
        # creating the new lower bound
        if (res@xy[cb,2] > 0) {
          # a new value has to be given
          res@xy[cb,2] <- res@xy[cb,2] + diff(res@xy[cb+(0:1),2]) *
                                         (res@sup[1]-res@xy[cb,1])/
                                         diff(res@xy[cb+(0:1),1]);
        }
        res@xy[cb,1] <- res@sup[1];
        res@xy[ 1,1] <- res@sup[1];
        # something to leave out
        if (cb > 2) {
          rem <- c(rem,2:(cb-1));
        }
      }
    }
    # continuous part on the upper side
    cb <- sum(res@xy[,1] > res@sup[2]);
    if (cb > 0) {
      # something has to be removed
      if (cb == nrow(res@xy)) {
        # everything
        res@xy <- matrix(0,0,2);
      } else {
        nr <- nrow(res@xy)-cb+1;
        # creating the new upper bound
        if (res@xy[nr,2] > 0) {
          # a new value has to be given
          res@xy[nr,2] <- res@xy[nr-1,2] + diff(res@xy[nr-1+(0:1),2]) *
                                         (res@xy[nr+1,1]-res@sup[2])/
                                         diff(res@xy[nr-1+(0:1),1]);
        }
        res@xy[nr     ,1] <- res@sup[2];
        res@xy[nr+cb-1,1] <- res@sup[2];
        # something to leave out
        if (cb > 2) {
          rem <- c(rem,nr+(1:(cb-2)));
        }
      }
    }
    # removing the useless row if any
    if (length(rem)>0) {
      res@xy <- res@xy[-rem,];
    }
  }
}
if (res@repa[2]>0) {
  #
  ## removing useless cordinates for the continuous part
  dd <- function(x) {
    x <- diff(x)==0;
    x <- (1-x)*bf(x);
    x <- diff(x)==0;
    1+which(x);
  }
  uu <- dd(res@xy[,2]);
  while (length(uu) > 0) {
    res@xy <- res@xy[-uu,,drop=FALSE];
    uu <- dd(res@xy[,2]);
  }
}
#
# in case of elimination a new normalization is necessary
if (eli) {
  res <- Recall(res);
}
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
range5empirical <- function(empi)
#TITLE  returns the ranges of an /empirical/
#DESCRIPTION
#   Returns the different supports of an /empirical/
# that is the minimum segment where stand
# the variation of the random variable,
# taking into account (or not) the \code{@sup} slot.
# See the section 'value' for details.
#DETAILS
# \code{empi} is not normalized, so the range
# can exceed the \code{@sup}.
#KEYWORDS 
#INPUTS
#{empi} <<The \code{empirical} to consider.>>
#[INPUTS]
#VALUE
# A list with four components:
# \code{$rcon}: the range of the continuous part
# (code{c(NA,NA)} if it does not exist).
# \code{$rdis}: the set of possible values for the
# discrete part (code{numeric(0)} if it does not exist).
# \code{$remp}: the range of both parts not taking into
# account the slot \code{@sup}.
# \code{$rang}: the range of both parts taking into
# account the slot \code{@sup}, identical to \code{$remp}
# after normalisation.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# range5empirical(new("empirical"));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_18
#REVISED 11_01_18
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
    erreur(che,"/empirical/ is not valid");
  }
}
#
# continuous part
if (empi@repa[2]>0) {
  rcon <- range(empi@xy[,1]);
} else {
  rcon <- c(NA,NA);
}
#
# discrete part
if (empi@repa[1]>0) {
  rdis <- unique(empi@xw[,1]);
} else {
  rdis <- numeric(0);
}
#
# global ranges
remp <- range(c(rcon,rdis),na.rm=TRUE);
rang <- range(c(remp,empi@sup));
#
# returning
list(rcon=rcon,rdis=rdis,remp=remp,rang=rang);
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
continuous2density <- function(xyco,x,dupli="m")
#TITLE  density of the continuous part
#DESCRIPTION
#   This function returns the density
# of a continous part of an /empirical/ for the
# values of \code{x}.
#DETAILS
# Values are not normalized, to obtain them
# just use \code{normalize8empirical} before
# extracting the continuous part.
#KEYWORDS 
#INPUTS
#{xyco} <<The \code{[email protected]} to consider.>>
#{x} << Values for which the values must be computed.>>
#[INPUTS]
#{dupli} << When a \code{x} value is exactly in a
#           vertical part of the density. Several choices
#           are possible monitored by \code{dupli}.
#           If "d", the point will be duplicated,
#              "1", the first value will be chosen,
#              "2", the second value will be chosen,
#              "m", the mean value will be chosen.>>
#VALUE
# The resulting values. When \code{dupli=="d"}, they
# are not returned as a vector but as a matrix the two columns
# of which are the possibly duplicated \code{x} and
# the resulting values.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# xx <- seq(0,6,0.3);
# xy <- matrix(c(1:5,0,1,3,3,0),5);
# cbind(xx,continuous2density(xy,xx));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_13
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  veri <- new("empirical",repa=0:1,xy=xyco,sup=range(xyco));
  che <- valid8empirical(veri);
  if (!identical(che,TRUE)) {
      erreur(che,"This 'xyco' is not valid");
  }
  #
  check4tyle(dupli,"character",1,c("d","1","2","m"),"'dupli' not accepted");
}
#
# degenerate cases
if (nrow(xyco) == 0) {
  return(rep(0,length(x)));
}
if (length(x) == 0) { return(numeric(0));}
#
ix <- findInterval(x,xyco[,1]);
#
res <- rep(0,length(x));
#
dup <- matrix(NA,0,3);
for (ii in bc(nrow(xyco)-1)) {
  uu <- which(ix == ii);
  if (length(uu)>0) {
    ma <- xyco[ii+(0:1),,drop=FALSE];
    if (ma[1,1]==ma[2,1]) {
      # the jumping case (vertical segment)
      res[uu] <- ma[2,2];
    } else {
      if (ma[1,2]==ma[2,2]) {
        # the horizontal case (uniform part)
        res[uu] <- mean(ma[,2]);
        if (dupli=="1") {res[uu] <- ma[1,2];}
        if (dupli=="2") {res[uu] <- ma[2,2];}
        if (dupli=="d") {
          res[uu] <- NA;
          dup <- rbind(dup,cbind(uu,matrix(ma[,2],length(uu),2,byrow=TRUE)));
        }
      } else {
        # the linearly [in|de]creasing case
        res[uu] <- ma[1,2] + diff(ma[,2]) *
                             (x[uu]-ma[1,1])/diff(ma[,1]);
      }
    }
  }
}
#
# the delicate case of duplication
if (dupli=="d") {
  res <- cbind(x,res);
  if (nrow(dup) > 0) {
    dup <- dup[order(dup[,1]),];
    res <- rbind(res,matrix(NA,nrow(dup),2));
    nn <- nrow(res);
    for (dd in rev(bc(nrow(dup)))) {
      ii <- dup[dd,1];
      i1 <- bd(ii,nn-1);
      res[i1+1,] <- res[i1,] ;
      res[ii,] <- dup[dd,2:3];
    }
  }
}
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
discrete2density <- function(xwdi,x)
#TITLE  density of the discrete part
#DESCRIPTION
#   This function returns the density
# of a discrete part of an /empirical/ for the
# values of \code{x}.
#DETAILS
# Values are not normalized, to obtain them
# just use \code{normalize8empirical} before
# extracting the continuous part.
#KEYWORDS 
#INPUTS
#{xwdi} <<The \code{[email protected]} to consider.>>
#{x} << Values for which the values must be computed.>>
#[INPUTS]
#VALUE
# The resulting values. If the \code{xwdi}
# was not normalized, possible multiple defined
# weight are cumulated.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# xx <- c(4,6,3);
# xw <- matrix(c(1:3,3:4,0,1,3,3,0),5);
# cbind(xx,ww=discrete2density(xw,xx));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_13
#REVISED 11_01_13
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  veri <- new("empirical",repa=1:0,xw=xwdi,sup=range(xwdi));
  che <- valid8empirical(veri);
  if (!identical(che,TRUE)) {
      erreur(che,"This 'xwdi' is not valid");
  }
}
#
# degenerate cases
if (nrow(xwdi) == 0) {
  return(rep(0,length(x)));
}
if (length(x) == 0) { return(numeric(0));}
#
# summing first
xx <- unique(xwdi[,1]);
if (nrow(xwdi) > length(xx)) {
  # some cumulation has to be done
  xxw <- cbind(xx,0);
  for (ix in bf(xx)) {
    xxw[ix,2] <- sum(xwdi[xwdi[,1]==xx[ix],2]);
  }
  xwdi <- xxw;
}
# computing second
res <- rep(0,length(x));
for (ix in bc(nrow(xwdi))) {
  res[x==xwdi[ix,1]] <- xwdi[ix,2];
}
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
continuous2cumulative <- function(xyco,x,norma=FALSE)
#TITLE  cumulative of the continuous part
#DESCRIPTION
#   This function returns the values of the cumulative
# of a continous part of an /empirical/ for the
# values of \code{x}.
#DETAILS
# Values are not normalized, to obtain them
# just use \code{normalize8empirical} before
# extracting the continuous part.
#KEYWORDS 
#INPUTS
#{xyco} <<The \code{[email protected]} to consider.>>
#{x} << Values for which the values must be computed.>>
#[INPUTS]
#{norma} <<Indicates if a normalization to one must be done.>>
#VALUE
# The resulting values.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# xx <- seq(0,6,0.3);
# xy <- matrix(c(1:5,0,1,3,3,0),5);
# cbind(xx,continuous2cumulative(xy,xx));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_14
#REVISED 11_01_21
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  veri <- new("empirical",repa=0:1,xy=xyco,sup=range(xyco[,1]));
  che <- valid8empirical(veri);
  if (!identical(che,TRUE)) {
      erreur(che,"This 'xyco' is not valid");
  }
}
#
# degenerate cases
if (nrow(xyco) == 0) {
  return(rep(0,length(x)));
}
if (length(x) == 0) { return(numeric(0));}
## computing some value for each interval
# number of intevals
nbi <- nrow(xyco)-1;
# slope of each interval (corrected for null increases of x)
ra <- diff(xyco[,2]) / diff(xyco[,1]);
ra[is.infinite(ra)] <- 0; 
# cumulative at the beginning of each interval
cu <- cumsum(diff(xyco[,1])*(xyco[-(nbi+1),2]+xyco[-1,2])/2);
cu <- c(0,cu[-nbi]);
## adding an x-value for the right outside elements
x <- c(x,max(xyco[,1]));
## affectation to the intervals
ix <- findInterval(x,xyco[,1]);
# upper bounds must be reaffected
ix[ix==nbi+1] <- nbi
## initialization of the result
res <- rep(0,length(x));
## dealing for each interval in turn
for (ii in bc(nbi)) {
  jj <- which(ix==ii);
  res[jj] <- cu[ii] + (x[jj]-xyco[ii,1])*
                      (xyco[ii,2]+(x[jj]-xyco[ii,1])*ra[ii]/2);
}
## putting right outside elements to the maximum
res[x>=max(xyco[,1])] <- res[length(res)];
if (norma) {
  ## normalizing
  res <- res / res[length(res)];
}
# removing the added element
res <- res[-length(res)];
#
## returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
discrete2cumulative <- function(xwdi,x,norma=FALSE)
#TITLE  cumulative of the discrete part
#DESCRIPTION
#   This function returns the values of the cumulative
# of a discrete part of an /empirical/ for the
# values of \code{x}.
#DETAILS
# Values are not normalized, to obtain them
# just use \code{normalize8empirical} before
# extracting the continuous part.
#KEYWORDS 
#INPUTS
#{xwdi} <<The \code{[email protected]} to consider.>>
#{x} << Values for which the values must be computed.>>
#[INPUTS]
#{norma} <<Indicates if a normalization to one must be done.>>
#VALUE
# The resulting values. If the \code{xwdi}
# was not normalized, possible multiple defined
# weight are cumulated.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# xx <- c(4,6,3);
# xw <- matrix(c(1:3,3:4,0,1,3,3,0),5);
# cbind(xx,ww=discrete2cumulative(xw,xx));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_14
#REVISED 11_01_14
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  veri <- new("empirical",repa=1:0,xw=xwdi,sup=range(xwdi));
  che <- valid8empirical(veri);
  if (!identical(che,TRUE)) {
      erreur(che,"This 'xwdi' is not valid");
  }
}
#
# degenerate cases
if (nrow(xwdi) == 0) {
  return(rep(0,length(x)));
}
if (length(x) == 0) { return(numeric(0));}
#
# summing first
xx <- unique(xwdi[,1]);
if (nrow(xwdi) > length(xx)) {
  # some cumulation has to be done
  xxw <- cbind(xx,0);
  for (ix in bf(xx)) {
    xxw[ix,2] <- sum(xwdi[xwdi[,1]==xx[ix],2]);
  }
  xwdi <- xxw;
}
# sorting second
oo <- order(xwdi[,1]);
xwdi <- xwdi[oo,,drop=FALSE];
# computing third
nbi <- nrow(xwdi)-1;
# adding an element for the right outside elements
x <- c(x,xwdi[nbi+1,1]);
# finding the affectations
ix <- findInterval(x,xwdi[,1]);
# initializing
res <- rep(0,length(x));
for (ii in bc(nbi+1)) {
  jj <- (ix==ii);
  res[jj] <- sum(xwdi[bc(ii),2]);
}
## putting right outside elements to the maximum
res[x>=max(xwdi[,1])] <- res[length(res)];
if (norma) {
  ## normalizing 
  res <- res[-length(res)] / res[length(res)];
}
# removing the added element
res <- res[-length(x)];
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
density2cumulative <- function(empi,N=1000)
#TITLE  returns the cumulative of an /empirical/
#DESCRIPTION
#   This function returns the cumulative probability
# distribution function of an /empirical/ object.
# That is returns the coordinates of points
# defining this function. Be aware that the calculus is
# made numerically, which is not the case for \code{pempirical},
# so better use it.
#DETAILS
# \code{N} must be understood
# as a simple way to give the level of precision the user
# wants to see apply for the returned curve.
#KEYWORDS
#INPUTS
#{empi} <<The \code{empirical} object to consider.>>
#[INPUTS]
#{N} <<The number of steps for the curve all
#      around the range of the random variable.>>
#VALUE
# A matrix with two columns \code{x,y}.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# density2cumulative(new("empirical"));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_13
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
    erreur(che,"/empirical/ is not valid");
  }
}
# degenerate case
if (sum(empi@repa)==0) { return();}
# computing the range
ra <- numeric(0);
if (empi@repa[1] > 0) {ra <- range(ra,empi@xw);}
if (empi@repa[2] > 0) {ra <- range(ra,empi@xy);}
# computing the points
N <- max(N,10);
xx <- seq(ra[1],ra[2],length=N);
if (empi@repa[1]>0) {
  xx <- sort(unique(c(xx,empi@xw[,1])));
}
# easying by the normalisation
empi <- normalize8empirical(empi);
# computing the values due to the continuous part
if (empi@repa[2] > 0) {
  yy <- cumsum(continuous2density(empi@xy,xx));
  yy <- yy / max(yy) * empi@repa[2];
} else {
  yy <- rep(0,length(xx));
}
# adding the jump for the discrete part
if (empi@repa[1]>0) {for (ii in bc(nrow(empi@xw))) {
  ou <- (xx >= empi@xw[ii,1]);
  yy[ou] <- yy[ou] + empi@xw[ii,2];
}}
# returning
matrix(cbind(xx,yy),ncol=2,dimnames=list(NULL,c("X","Y")));
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
solve2degree <- function(y,ky,kx2,kx,kk,dx=NULL,x0=NULL)
#TITLE  solves a 'generalized' degree two polynomial
#DESCRIPTION
#   This function returns the root (or two roots) of 
# the equation \code{ky*y + kx2*x^2 + kx*x + kk = 0}.
# When \code{dx} is not null, it is supposed to give
# the interval where the root lies, in that case only
# one root is returned.\cr
# The first parameter can be a vector of any
# length and all computations are vectorized.\cr
# Only real roots are considered.
#DETAILS
# When \code{dx} is defined  only one root is returned,
# belonging to the interval; if it is not possible (root(s)
# exist(s) and do(es) not comply, then a fatal error
# is issued.\cr
# When every real number complies with the equation, according
# to available arguments, the returning value is \code{x0},
# \code{mean(dx)} or \code{0}.
# When \code{is.null(dx)} either one or two roots is 
# returned with \code{NA} when the solution is complex.
#KEYWORDS
#INPUTS
#{y} <<Vector of values for which the equation must be satisfied.>>
#{ky} <<Coefficient for \code{y}.>>
#{kx2} <<Coefficient for \code{x^2}.>>
#{kx} <<Coefficient for \code{x}.>>
#{kk} <<Constant coefficient.>>
#[INPUTS]
#{dx} <<\code{NULL} or the interval (\code{numeric(2)}) for the roots.>>
#{x0} <<\code{NULL} or a proposal in case of indetermination.>>
#VALUE
# A matrix having one or two columns according to the values of
# \code{ky,kx2,kx,kk,dx}.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# solve2degree(1:10, 1,1,0,-20);
# solve2degree(   3,-1,1,1,  1);
# solve2degree(   3,-1,1,1,  1,c(0.5,1.5));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_18
#REVISED 11_01_21
#--------------------------------------------
{
# some checking
if (rbsb.mck) {
  check4tyle(y,"numeric",-1,message="solve2degree: non accepted 'y'");
  check4tyle(ky,"numeric",1,message="solve2degree: non accepted 'ky'");
  check4tyle(kx2,"numeric",1,message="solve2degree: non accepted 'kx2'");
  check4tyle(kx,"numeric",1,message="solve2degree: non accepted 'kx'");
  check4tyle(kk,"numeric",1,message="solve2degree: non accepted 'kk'");
  if (!is.null(dx)) {
    check4tyle(dx,"numeric",2,message="solve2degree: non accepted 'dx'");
  }
  if (!is.null(x0)) {
    check4tyle(x0,"numeric",1,message="solve2degree: non accepted 'x0'");
  }
}
# number of equations
ne <- length(y);
# degenerate case
if (ne==0) { return(numeric(0));}
# modifying the constant
kk <- ky*y + kk;
#
# exploring the case
if (kx2==0) {
  # 1rst degree at most
  if (kx==0) {
    # 0 degree
    if (all(kk==0)) {
      # any real is root
      if (!is.null(x0)) {
        res <- matrix(x0,ne,1);
      } else {
        if (!is.null(dx)) {
          res <- matrix(mean(dx),ne,1);
        } else {
          res <- matrix(0,ne,1);
        }
      }
    } else {
      # no root
      res <- matrix(NA,ne,1);
      erreur(list(ky,kx2,kx,kk),"solve2degree: no solution for the proposed equation");
    }
  } else {
    # 1rst degree
    res <- matrix(-kk/kx,ne,1);
  }
} else {
  # 2d degree
  disc <- kx^2 - 4*kx2*kk;
  rro <- (disc >= 0);
  res <- matrix(NA,ne,2);
  res[rro,] <- (-kx + outer(sqrt(disc[rro]),c(-1,1),"*")) / (2*kx2);
  if (!is.null(dx)) {
    ou1 <- (res[rro,1]-dx[1])*(res[rro,1]-dx[2]);
    ou2 <- (res[rro,2]-dx[2])*(res[rro,2]-dx[1]);
    ou <- (ou1 <= ou2);
    res[rro[!ou],1] <- res[rro[!ou],2];
    res <- res[,1,drop=FALSE];
  }
}
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
dempirical <- function(x,empi)
#TITLE  returns the density of an /empirical/
#DESCRIPTION
#   This function returns the probability
# density of an /empirical/ object as a list
# with two components: \code{w} for the discrete part
# and \code{y} for the continuous part associated
# to different values of the random variable.
#DETAILS
# the proposed \code{empi} is normalized before
# the computation. So if you don't what that, use
# by yourself the two functions \code{continuous2density}
# and \code{discrete2density}.
#KEYWORDS
#INPUTS
#{x} <<Vector of values for which the density is desired.>>
#{empi} <<The \code{empirical} object to consider.>>
#[INPUTS]
#VALUE
# A list with two components \code{$w,$y} having each the
# same length as \code{x};
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# dempirical(seq(0,4,0.5),new("empirical"));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_13
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
      erreur(che,"/empirical/ is not valid");
  }
}
# degenerate case
if (length(x)==0) { return(list(w=numeric(0),y=numeric(0)));}
# normalizing
empi <- normalize8empirical(empi);
# the continuous part
if (empi@repa[2]>0) {
  yy <- continuous2density(empi@xy,x);
} else {
  yy <- rep(0,length(x));
}
# the discrete part
if (empi@repa[1]>0) {
  ww <- discrete2density(empi@xw,x);
} else {
  ww <- rep(0,length(x));
}
#
# returning
list(w=ww,y=yy);
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
pempirical <- function(x,empi)
#TITLE  returns the distribution function of an /empirical/
#DESCRIPTION
#   This function returns the cumulative distribution
# function of an /empirical/ object for different values of
# the random variable.
#DETAILS 
#KEYWORDS
#INPUTS
#{x} <<Vector of values for which the distribution is wanted.>>
#{empi} <<The \code{empirical} object to consider.>>
#[INPUTS]
#VALUE
# A vector with same length as \code{x};
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# pempirical(seq(0,4,0.5),new("empirical"));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_10
#REVISED 11_01_18
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
      erreur(che,"/empirical/ is not valid");
  }
}
# degenerate case
if (length(x)==0) { return(numeric(0));}
# normalizing
empi <- normalize8empirical(empi);
#
# the continuous part
if (empi@repa[2]>0) {
  res <- continuous2cumulative(empi@xy,x,FALSE);
} else {
  res <- rep(0,length(x));
}
# the discrete part
if (empi@repa[1]>0) {
  res <- res + discrete2cumulative(empi@xw,x,FALSE);
}
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
qempirical <- function(p,empi)
#TITLE  returns the quantiles of an /empirical/
#DESCRIPTION
#   This function returns the quantiles of
# an /empirical/ object for different probabilities.
#DETAILS 
#KEYWORDS
#INPUTS
#{p} <<Vector of probabilities for which the quantile are wanted.>>
#{empi} <<The \code{empirical} object to consider.>>
#[INPUTS]
#VALUE
# A vector with same length as \code{p} containing the quantiles.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# pempirical(seq(0,1,0.2),new("empirical"));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_19
#REVISED 11_01_24
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  check4tyle(p,"numeric",-1,c(0,1),"qempirical: some 'p' are not probabilities");
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
    erreur(che,"/empirical/ is not valid");
  }
}
# degenerate case
if (length(p)==0) { return(numeric(0));}
# normalization
empi <- normalize8empirical(empi);
#
# describing the distribution function with pieces of parabolae
#
# getting the abscissae
xx <- sort(unique(c(empi@xw[,1],empi@xy[,1])));
#
# getting the ordinates
if (empi@repa[2] == 0) {
  # only the discrete part
  yy <- discrete2cumulative(empi@xw,xx);
} else {
  # for the continuous part
  yy <- continuous2cumulative(empi@xy,xx);
  if (empi@repa[1] > 0) {
    # also must be added the discrete part
    for (ii in bc(nrow(empi@xw))) {
      # looking for the place into xx where to insert
      # the corresponding probability
      xi <- empi@xw[ii,1];
      jj <- which(xx==xi)[1];
      # adding it on both xx and yy
      xx <- c(xx[1:jj],xi,xx[bd(jj+1,length(xx))]);
      yy <- c(yy[1:jj],yy[bd(jj,length(yy))]+empi@xw[ii,2]);
    }
  }
}
#
# adding a starting point when the discrete part
# has got the lowest value
if (yy[1] > 0) {
  xx <- c(xx[1],xx);
  yy <- c(    0,yy);
}
dx <- diff(xx);
dy <- diff(yy);
#
# Classifying the different types of intervals
nbi <- length(xx) - 1;
# horizontal cases
# (= no proba at all)
cah <- which(dy==0);
# vertical cases
# (= discrete proba)
cav <- which(dx==0);
# increasing case cases
if (empi@repa[2]>0) {
  # (= continuous proba)
  cas <- setdiff(1:nbi,union(cah,cav));
  d2x <- cbind(xx[cas],xx[cas])+t(outer((1:2)/3,dx[cas],"*"));
  vco <- matrix(continuous2density(empi@xy,d2x),ncol=2);
  d2y <- vco[,2]-vco[,1];
  d2x <- d2x[,2]-d2x[,1];
  d2x[d2x==0] <- 1; 
  cde <- cas[d2y <0]; # decreasing ones
  cco <- cas[d2y==0]; #  constant ones
  cin <- cas[d2y >0]; #  inreasing ones
  # introducing the derivative in the same way
  ax <- rep(0,nbi);
  ax[cas] <- d2y/d2x;
  d2y <- ax;
} else {
  cde <- cco <- cin <- numeric(0);
  # and discrete weight must be duplicated
  xx <- rep(xx,each=2);
  yy <- rep(yy,each=2);
  xx <- xx[-c(1,3)];
  yy <- yy[-c(1,length(yy))];
  # updating in consequence the classification
  dx <- diff(xx);
  dy <- diff(yy);
  cah <- which(dy==0);
  cav <- which(dx==0);
}
# coding the type of curve of each different interval
# from starting points 'cbind(xx,yy,c(indica,10));'
indica <- rep(0,nbi);
indica[cin] <- indica[cin ] + 1;
indica[cco] <- indica[cco ] + 2;
indica[cde] <- indica[cde ] + 3;
indica[cav] <- indica[cav ] + 4;
## getting the corresponding coefficient for the curve pieces
# initializing with the horizontal case
nbir <- length(yy)-1;
ky <- rep(-1,nbir);
kx2 <- kx <- rep(0,nbir);
kk <- - yy[1:nbir];
if (empi@repa[2]>0) {
  # straight case
  kx[cco] <-   dy[cco] / dx[cco];
  kk[cco] <-   yy[cco] - dy[cco] / dx[cco] * xx[cco];
  # increasing and decreasing case
  ca2 <- c(cde,cin);
  kx2[ca2] <- d2y[ca2]/2;
  kx[ca2] <- dy[ca2]/dx[ca2] - kx2[ca2]*(xx[ca2]+xx[ca2+1]);
  kk[ca2] <- (xx[ca2]*yy[ca2+1]-xx[ca2+1]*yy[ca2]) / (xx[ca2]-xx[ca2+1]) +
             kx2[ca2]*xx[ca2]*xx[ca2+1];
}
# vertical case
ky[cav] <-  0;
kx[cav] <- -1;
kk[cav] <- xx[cav];
#print(cbind(xx[-(nbir+1)],yy[-(nbir+1)],indica,d2y,dy,dx,ky,kx2,kx,kk));
#
# starting the inversion
#
# getting the partition over the nbir intervals
ou <- findInterval(p,yy);
ou[ou>nbir] <- nbir;
ou[ou==0] <- 1;
#
# looping on the nbir interval
res <- rep(NA,length(p));
for (ii in bc(nbir)) { if (ii %in% ou) {
  oui <- which(ou==ii);
  res[oui] <-  solve2degree(p[oui],
                              ky[ii],kx2[ii],kx[ii],kk[ii],
                              xx[ii+(0:1)],xx[ii]
                             );
}}
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
rempirical <- function(n,empi)
#TITLE  returns pseudo-random draws from an /empirical/
#DESCRIPTION
#   This function returns \code{n} independent pseudo-random
# draw from the /empirical/ object \code{empi}. The user is advised
# to manage the pseudo-random seed before calling this function.
# The object is normalized before being used (see \code{normalize8empirical}).
#DETAILS
# Just a call to \code{qempirical} with draw from the basic
# \code{runif}. But the bad convention of using \code{length(n)}
# when it is greater than 1 is not followed in case of \code{rbsb.mck}.
#KEYWORDS
#INPUTS
#{n} <<Wanted number of draws.>>
#{empi} <<The \code{empirical} object to consider.>>
#[INPUTS]
#VALUE
# A vector with length \code{n} containing the drawn values.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# dempirical(100,new("empirical"));
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_23
#REVISED 11_01_23
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  check4tyle(n,"numeric",1,c(0,Inf),message="qempirical: some 'p' are not probabilities");
  che <- valid8empirical(empi);
  if (!identical(che,TRUE)) {
    erreur(che,"/empirical/ is not valid");
  }
}
# degenerate case
if (n==0) { return(numeric(0));}
# normalization
empi <- normalize8empirical(empi);
#
# drawing uniform values
yy <- runif(n);
#
# deducing the values
res <- qempirical(yy,empi);
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
triangular2empirical <- function(a=-1,b=0,c=1,des=NULL)
#TITLE  returns the /empirical/ associated to a triangular distribution
#DESCRIPTION
#   This function returns  the /empirical/ object corresponding
# to the triangular distribution where \code{a} is the minimum value,
# \code{b} is the mode value and \code{c} is the maximum value.
#DETAILS
# The resulting object is normalized.
#KEYWORDS
#INPUTS
#[INPUTS]
#{a} <<Minimal value>>
#{b} <<Most probable value>>
#{c} <<Maximal value>>
#{des} << /des/ to give to the object, when \code{NULL}
# a standard description will be provided.>>
#VALUE
# An /empirical/ object.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# plot(triangular2empirical(1,5,6),what="d");
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_24
#REVISED 11_01_24
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  check4tyle(a,"numeric",1,message="triangular2empirical: bad 'a'");
  check4tyle(b,"numeric",1,message="triangular2empirical: bad 'b'");
  check4tyle(c,"numeric",1,message="triangular2empirical: bad 'c'");
  if(!all(diff(c(a,b,c))>=0)) {
    erreur(c(a,b,c),"The three arguments must be non decreasing");
  }
}
# creation
if (is.null(des)) {
  des <- char2des("Triangular");
}
res <- new("empirical",des=des,
                       sup=c(a,c),
                       repa=c(0,1),
                       xw=matrix(NA,0,2),
                       xy=matrix(c(a,b,c,0,1,0),3,2)
          );
# normalization
res <- normalize8empirical(res);
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
trapezoid2empirical <- function(a=-1,b=-0.5,c=0.5,d=1,des=NULL)
#TITLE  returns the /empirical/ associated to a trapezoid distribution
#DESCRIPTION
#   This function returns  the /empirical/ object corresponding
# to the trapezoid distribution where \code{a} is the minimum value,
# \code{b} and \code{c} are the mode values and \code{d} is the maximum value.
#DETAILS
# The resulting object is normalized.
#KEYWORDS
#INPUTS
#[INPUTS]
#{a} <<Minimal value>>
#{b} <<Beginning of the mode interval>>
#{c} <<Ending of the mode interval>>
#{d} <<Maximal value>>
#{des} << /des/ to give to the object, when \code{NULL}
# a standard description will be provided.>>
#VALUE
# An /empirical/ object.
#EXAMPLE
# rbsb3k("RESET"); # for R checking compliance (useless)
# plot(trapezoid2empirical(1,5,6,6),what="d");
#REFERENCE
#SEE ALSO
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_01_24
#REVISED 11_01_24
#--------------------------------------------
{
if (rbsb.mck) {
  # some checking
  check4tyle(a,"numeric",1,message="trapezoid2empirical: bad 'a'");
  check4tyle(b,"numeric",1,message="trapezoid2empirical: bad 'b'");
  check4tyle(c,"numeric",1,message="trapezoid2empirical: bad 'c'");
  check4tyle(d,"numeric",1,message="trapezoid2empirical: bad 'd'");
  if(!all(diff(c(a,b,c,d))>=0)) {
    erreur(c(a,b,c),"The four arguments must be non decreasing");
  }
}
# creation
if (is.null(des)) {
  des <- char2des("Trapezoid");
}
res <- new("empirical",des=des,
                       sup=c(a,d),
                       repa=c(0,1),
                       xw=matrix(NA,0,2),
                       xy=matrix(c(a,b,c,d,0,1,1,0),4,2)
          );
# normalization
res <- normalize8empirical(res);
#
# returning
res;
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Try the rbsb package in your browser

Any scripts or data that you put into this service are public.

rbsb documentation built on May 2, 2019, 4:41 p.m.