Nothing
R/f3.rbsb4.code.r
In rbsb: set of ancillary functions
Defines functions
valid8empirical
print8empirical
plot8empirical
normalize8empirical
range5empirical
continuous2density
discrete2density
continuous2cumulative
discrete2cumulative
density2cumulative
solve2degree
dempirical
pempirical
qempirical
rempirical
triangular2empirical
trapezoid2empirical
Documented in
continuous2cumulative
continuous2density
dempirical
density2cumulative
discrete2cumulative
discrete2density
normalize8empirical
pempirical
plot8empirical
print8empirical
qempirical
range5empirical
rempirical
solve2degree
trapezoid2empirical
triangular2empirical
valid8empirical
###########################################
###########################################
########
#((((((( 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{empirical@xy} 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{empirical@xw} 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{empirical@xy} 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{empirical@xw} 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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions valid8empirical print8empirical plot8empirical normalize8empirical range5empirical continuous2density discrete2density continuous2cumulative discrete2cumulative density2cumulative solve2degree dempirical pempirical qempirical rempirical triangular2empirical trapezoid2empirical
Documented in continuous2cumulative continuous2density dempirical density2cumulative discrete2cumulative discrete2density normalize8empirical pempirical plot8empirical print8empirical qempirical range5empirical rempirical solve2degree trapezoid2empirical triangular2empirical valid8empirical
###########################################
###########################################
########
#((((((( 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{empirical@xy} 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{empirical@xw} 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{empirical@xy} 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{empirical@xw} 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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.