Nothing
#' Joint MoTBF density learning
#'
#' Two functions for learning joint MoTBFs. The first one, \code{parametersJointMoTBF()},
#' gets the parameters by solving a quadratic optimization problem, minimizing
#' the mean squared error between the empirical joint CDF and the estimated CDF.
#' The density is obtained as the derivative od the estimated CDF.
#' The second one, \code{jointMoTBF()}, fixes the equation of the joint function using
#' the previously learned parameters and converting this \code{"character"} string into an
#' object of class \code{"jointmotbf"}.
#'
#' @name jointmotbf.learning
#' @rdname jointmotbf.learning
#' @param X A dataset of class \code{"data.frame"}.
#' @param ranges A \code{"numeric"} matrix containing the range of the varibles used to fit the function,
#' where each column corresponds to a variable. If not specified, the range of each variable is computed from the data.
#' @param dimensions A \code{"numeric"} vector containing the number of parameters of each varible.
#' @param object A list with the output of the function \code{parametersJointMoTBF()}.
#' @return
#' \code{parametersJointMoTBF()} returns a list with the following elements:
#' \bold{Parameters}, which contains the computed coefficients of the resulting function;
#' \bold{Dimension}, which is a \code{"numeric"} vector containing the number
#' of coefficients used for each variable;
#' \bold{Range} contains a \code{"numeric"} matrix with the domain of each variable, by columns;
#' \bold{Iterations} contains the number of iterations needed to solve the problem;
#' \bold{Time} contains the execution time.
#'
#' \code{jointMoTBF()} returns an object of class \code{"jointmotbf"}, which is a list whose only visible element
#' is the analytical expression of the learned density. It also contains the other aforementioned elements,
#' which can be retrieved using \code{attributes()}
#' @importFrom Matrix nearPD
#' @examples
#
#' ## 1. EXAMPLE
#' ## Generate a multinormal dataset
#' data <- data.frame(X1 = rnorm(100), X2 = rnorm(100))
#'
#' ## Joint learnings
#' dim <- c(2,3)
#' param <- parametersJointMoTBF(X = data, dimensions = dim)
#'
#' param$Parameters
#' length(param$Parameters)
#' param$Dimension
#' param$Range
#'
#' P <- jointMoTBF(param)
#' P
#' attributes(P)
#' class(P)
#'
#' ###############################################################################
#' ## MORE EXAMPLES ##############################################################
#' ###############################################################################
#' \donttest{
#' ## Generate a dataset
#' data <- data.frame(X1 = rnorm(100), X2 = rnorm(100), X3 = rnorm(100))
#'
#' ## Joint learnings
#' dim <- c(3,2,3)
#' param <- parametersJointMoTBF(X = data, dimensions = dim)
#'
#' param$Parameters
#' length(param$Parameters) ## prod(dim)
#' param$Dimension
#' param$Range
#' param$Time
#'
#' P <- jointMoTBF(param)
#' P
#' attributes(P)
#' class(P)
#' }
#' @rdname jointmotbf.learning
#' @export
# parametersJointMoTBF=function(X, ranges=NULL, dimensions=NULL, fitPoints = 10, constraints = 10)
parametersJointMoTBF=function(X, ranges=NULL, dimensions=NULL)
{
tm <- Sys.time()
if(is.null(ranges)) ranges <- sapply(X, range)
if(is.null(dimensions)){
Fx <- lapply(X, univMoTBF, "MOP")
ll <- list(); for(i in 1:length(Fx)) ll[[length(ll)+1]] <- length(coef(Fx[[i]]))
} else {
ll <- list(); for(i in 1:length(dimensions)) ll[[length(ll)+1]] <- dimensions[i]
}
dim <- c(); pos <- which(c(length(ll[[1]]), length(ll[[2]]))==min(length(ll[[1]]), length(ll[[2]])))
if(pos[1]==1) {pos1 <- 1; pos2 <- 2}
if(pos[1]!=1) {pos1 <- 2; pos2 <- 1}
for(i in 1:length(ll[[pos1]])) for(j in 1:length(ll[[pos2]])) dim <- c(dim,list(c(ll[[pos1]][i],ll[[pos2]][j])))
if(pos1==2) for(i in 1:length(dim)) dim[[i]] <- dim[[i]][length(dim[[i]]):1]
size <- prod(sapply(1:length(ll), function(i) length(ll[[i]])))
si <- lapply(1:ncol(X),function(i) rep(ll[[i]][1], size))
dim <- c(); for(i in 1:length(si)) dim <- cbind(dim,unlist(si[[i]]))
dim <- lapply(1:nrow(dim), function(i) dim[i,])
## Create the grid: depending on the n.record and the n.variables
#npointsgrid <- 100
# npointsgrid <- max(ceiling((2/(ncol(X)))^2*100), 10)
# npointsgrid <- 10
fitPoints <- max(ceiling((2/(ncol(X)))^2*100), 10)
eg <- list()
for(i in 1: ncol(X)){
eg[[i]] <- seq(ranges[1,i], ranges[2,i], length.out = fitPoints)
}
x <- expand.grid(eg)
## Cumulative densities
y <- jointCDF(X,x)
P <- c()
for(s in 1:length(dim)){
Xt <- ranges; pos <- coefExpJointCDF(dim[[s]])
n <- nrow(x); nterms <- length(pos)
# x1 <- c()
# for(j in 1:ncol(X)){
# xx=rep(1, n)
# for(i in 1:nterms){
# xi <- cbind((x[,j]^pos[[i]][j]))
# xx <- cbind(xx,xi)
# }
# x1[[length(x1)+1]] <- xx
# }
x1 <- list()
for(j in 1:ncol(X)){
pos2 <- c(0,unlist(lapply(pos, function(x){x[j]})))
x1[[j]] <- outer(x[,j], pos2, FUN = "^")
}
xx <-Reduce("*", x1)
## Dmat
XX <- t(xx)%*%xx
## dvec
Xy <- t(xx)%*%y
## Constrains
xNew <- list()
constraints = round(60000^(1/ncol(Xt)))
for( j in 1:ncol(Xt)){
# xnew <- seq(min(Xt[,j]), max(Xt[,j]), length=round(60000^(1/ncol(Xt))))
xnew <- seq(min(Xt[,j]), max(Xt[,j]), length=constraints)
if(xnew[length(xnew)]!=max(Xt[,j])) xnew <- c(xnew,max(Xt[,j]))
xNew[[length(xNew)+1]] <- xnew
}
xN <- expand.grid(xNew)
ma22 <- c()
for(i in 1:ncol(X)){
d <- sapply(1:nterms, function(j) pos[[j]][i]*(xN[,i]^(pos[[j]][i]-1)))
tf <- sapply(1:length(pos), function(j) pos[[j]][i]==0)
if(any(tf)) d[,which(tf)] <- 0
d <- cbind(rep(0, nrow(xN)),d)
ma22[[length(ma22)+1]] <- t(d)
}
ma2 <- Reduce("*", ma22)
ma3 <- c()
for(j in 1:ncol(Xt)){
ma <- 0
for(i in 1:nterms) ma <- rbind(ma,(max(Xt[,j])^pos[[i]][j]) - (min(Xt[,j])^pos[[i]][j]))
ma3 <- cbind(ma3,ma)
}
for(i in 1:(ncol(X)-1)) ma3[,i+1] <- ma3[,i] * ma3[,i+1]
ma33 <- cbind(ma3[,i+1])
## Amat
AA <- cbind(ma33, ma2)
## bvec
B <- c(1,rep(1.0E-5,ncol(ma2)))
tr <- tryCatch(solve.QP(XX, Xy, AA, B, meq=1), error = function(e) NULL)
if(is.null(tr)==T){
message("matrix D in quadratic function has been approximated to the nearest positive definite!")
tr <- solve.QP(nearPD(XX)$mat, Xy, AA, B, meq=1)
# return(NULL)
}
finaltm <- Sys.time() - tm
soluc <- tr
parameters <- soluc$solution
# if(is.null(tr)==T){
# return(NULL)
# }else{
# soluc <- tr
# parameters <- soluc$solution
# }
## Parameters PDF
a <- sapply(pos, prod)
param <- parameters[-1]*a
parameters <- param[which(param!=0)]
P <- list(Parameters = parameters,
Dimension = dim[[s]], Range = Xt,
Iterations = tr$iterations[1],
Time=finaltm)
}
return(P)
}
#' @rdname jointmotbf.learning
#' @export
jointMoTBF <- function(object){
parameters <- object$Parameters
dimensions <- object$Dimension
v <- letters[24:(24+length(dimensions)-1)]
v[is.na(v)] <- letters[1:length(v[is.na(v)])]
if(length(parameters)==1) {
s <- paste(parameters[1], "+0", sep="")
for(i in 1:length(v)) s <- paste(s, "*", v[i], sep="")
P <- list(Function = s, Domain = object$Range)
P <- jointmotbf(P)
return(P)
}
pos <- coefExpJointCDF(dimensions)
t <- list()
for(i in 1:length(pos)) if(prod(pos[[i]])!=0)
t[[length(t)+1]] <- pos[[i]] else next
pos <- t
if(any(dimensions==1)) posvardim1 <- which(dimensions==1)
if(parameters[1]==0) s <- c() else s <- parameters[1]
for(i in 2:length(parameters)){
if(is.null(s)){
s <- paste(s, parameters[i],sep="")
} else{
s <- paste(s, ifelse(parameters[i]>=0, "+", ""), parameters[i],sep="")
if(i != length(parameters)) for(p in 1:length(v)) s <- paste(s, ifelse((pos[[i]][p]-1)==0, "", paste("*",v[p], ifelse((pos[[i]][p]-1)==1, "", paste("^", pos[[i]][p]-1, sep="")), sep="")),sep="")
if(i == length(parameters)){
for(p in 1:length(v)) s <- paste(s, ifelse((pos[[i]][p]-1)==0, ifelse(any(p==posvardim1),
paste("*",v[p],"^",pos[[i]][p]-1, sep=""),""),
paste("*",v[p], ifelse((pos[[i]][p]-1)==1, "",
paste("^", pos[[i]][p]-1, sep="")), sep="")),sep="")
}
}
}
P <- noquote(s)
P <- list(Function = P, Domain = object$Range,
Iterations = object$Iterations, Time = object$Time)
P <- jointmotbf(P)
return(P)
}
#' Degree Function
#'
#' Compute the degree for each term of a joint CDF.
#'
#' @param dimensions A \code{"numeric"} vector including the number of parameters of each variable.
#' @return A list with n element. Each element contains a \code{numeric} vector with the degree for
#' each variable and each term of the joint CDF.
#' @export
#' @examples
#'
#' ## Dimension of the joint PDF of 2 variables
#' dim <- c(4,5)
#' ## Potentials of each term of the CDF
#' c <- coefExpJointCDF(dim)
#' length(c) + 1 ## plus 1 because of the constant coefficient
#'
#' ## Dimension of the joint density function of 2 variables
#' dim <- c(5,5,3)
#' ## Potentials of the cumulative function
#' coefExpJointCDF(dim)
#'
coefExpJointCDF <- function(dimensions)
{
t <- lapply(1:length(dimensions), function(i) 0:dimensions[i])
l <- sapply(1:length(t), function(i) length(t[[i]]))
mm <- c()
for(i in 1:length(t)){
m <- c()
for(j in 1:length(t[[i]])) m <- c(m,rep(t[[i]][j], prod(l[-c(1:i)])))
if(i!=1) m <- rep(m, l[i-1])
mm <- cbind(mm, m)
}
mm <- mm[-1,]
colnames(mm) <- NULL
if(is.matrix(mm)) pos <- lapply(1:nrow(mm), function(i) mm[i,])
else pos <-list(as.numeric(mm))
return(pos)
}
#' Joint MoTBFs CDFs
#'
#' Function to compute multivariate CDFs.
#'
#' @name jointCDF
#' @rdname jointCDF
#' @param df The dataset as an object of class \code{data.frame}.
#' @param grid a \code{data.frame} with the selected data points where the objective function
#' will be evaluated when optimizing the parameters.
#' @return \code{jointCDF()} returns a vector.
#' @examples
#'
#' ## Create dataset with 2 variables
#' n = 2
#' size = 50
#' df <- as.data.frame(matrix(round(rnorm(size*n),2), ncol = n))
#'
#' ## Create grid dataset
#' npointsgrid <- 10
#' ranges <- sapply(df, range)
#' eg <- list()
#' for(i in 1: ncol(df)){
#' eg[[i]] <- seq(ranges[1,i], ranges[2,i], length.out = npointsgrid)
#' }
#'
#' x <- expand.grid(eg)
#'
#' ## Joint cumulative values
#' jointCDF(df = df, grid = x)
#'
#' @rdname jointCDF
#' @export
#'
jointCDF <- function(df, grid){
n <- ncol(df)
if(ncol(df)>10){
apply(grid, MARGIN=1,
FUN = function(grid,df){
b = which(colSums(t(df)<=grid)==n)
out = length(b)/nrow(df)
},
df=df)
}else{
apply(grid, MARGIN=1,
FUN = function(grid,df){
n <- ncol(df)
b = switch(n-1,
which(df[,1]<=grid[1] & df[,2]<=grid[2]),
which(df[,1]<=grid[1] & df[,2]<=grid[2]&df[,3]<=grid[3]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4] & df[,5]<=grid[5]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4] & df[,5]<=grid[5] & df[,6]<=grid[6]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4] & df[,5]<=grid[5] & df[,6]<=grid[6] & df[,7]<=grid[7]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4] & df[,5]<=grid[5] & df[,6]<=grid[6] & df[,7]<=grid[7] & df[,8]<=grid[8]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4] & df[,5]<=grid[5] & df[,6]<=grid[6] & df[,7]<=grid[7] & df[,8]<=grid[8] & df[,9]<=grid[9]),
which(df[,1]<=grid[1] & df[,2]<=grid[2] & df[,3]<=grid[3] & df[,4]<=grid[4] & df[,5]<=grid[5] & df[,6]<=grid[6] & df[,7]<=grid[7] & df[,8]<=grid[8] & df[,9]<=grid[9] & df[,10]<=grid[10])
)
out = length(b)/nrow(df)
},
df=df)
}
}
# Alternative to jointCDF, slower.
# jointCDF <- function(data, grid) unlist(lapply(1:nrow(grid), posGrid, data, grid))
#
# posGrid <- function(nrows, data, grid)
# {
# p <- lapply(1:ncol(grid), function(i) which(data[,i]<=grid[nrows,i]))
# pp <- c()
# for(i in 1:length(p)){
# pp <- c(pp,p[[i]])
# if(i!=1) pp <- pp[duplicated(pp, last=T)]
# }
# return(length(pp)/nrow(data))
# }
#' Coerce a \code{"jointmotbf"} Object to a Function
#'
#' Takes a \code{"jointmotbf"} object and contructs an \R function to evaluate it at multidimensional points.
#'
#' @param x An object of class \code{"joinmotbf"}.
#' @param \dots Further arguments to be passed to or from the method. Not necessary for this method.
#' @details This is an \code{S3} method for the generic function \link{as.function}.
#' @return It returns a function to evaluate an object of class \code{"jointmotbf"}.
#' @seealso \link{parametersJointMoTBF} and \link{jointMoTBF}
#' @export
#' @examples
#' ## 1.EXAMPLE
#' ## Dataset
#' data <- data.frame(X = rnorm(100), Y = rexp(100))
#'
#' ## Joint function
#' dim <- c(3,2)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#' density <- as.function(P)(data[,1], data[,2])
#' density
#'
#' ## Log-likelihood
#' sum(log(density))
#'
#' #############################################################################
#' ## MORE EXAMPLES ############################################################
#' #############################################################################
#' \donttest{
#' ## Dataset
#' data <- data.frame(X = rnorm(100), Y = rexp(100), Z = rnorm(100))
#'
#' ## Joint function
#' dim <- c(2,3,4)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#' density <- as.function(P)(data[,1], data[,2], data[,3])
#' density
#'
#' ## Log-likelihood
#' sum(log(density))
#' }
as.function.jointmotbf <- function(x, ...)
{
P <- x[[1]]
v <- nVariables(P)
f <- function(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,
s,t,u,v,w,x,y,z) eval(parse(text=P))
formals(f) <- formals(f)[1:length(v)]
names(formals(f)) <- v
return(f)
}
#' Coefficients of a \code{"jointmotbf"} object
#'
#' Extracts the parameters of a joint MoTBF density.
#'
#' @param object An MoTBF function.
#' @param \dots Other arguments, unnecessary for this function.
#' @return A \code{"numeric"} vector with the parameters of the function.
#' @seealso \link{parametersJointMoTBF} and \link{jointMoTBF}
#' @export
#' @examples
#' ## Generate a dataset
#' data <- data.frame(X1 = rnorm(100), X2 = rnorm(100))
#'
#' ## Joint function
#' dim <-c(2,4)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#' P$Time
#'
#' ## Coefficients
#' coef(P)
#'
#' #############################################################################
#' ## MORE EXAMPLES ############################################################
#' #############################################################################
#' \donttest{
#' ## Generate a dataset
#' data <- data.frame(X1 = rnorm(100), X2 = rnorm(100), X3 = rnorm(100))
#'
#' ## Joint function
#' dim <-c(2,4,3)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#' P$Time
#'
#' ## Coefficients
#' coef(P)
#' }
coef.jointmotbf <- function(object, ...) coeffMOP(object)
#' Number of Variables in a Joint Function
#'
#' Compute the number of variables which are in a \code{jointmotbf} object.
#'
#' @param P An \code{"motbf"} object or a \code{"jointmotbf"} object.
#' @return A \code{"character"} vector with the names of the variables in the function.
#' @export
#' @examples
#'
#' # 1. EXAMPLE
#' ## Generate a dataset
#' data <- data.frame(X1 = rnorm(100), X2 = rnorm(100))
#'
#' ## Joint function
#' dim <-c(3,2)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#' P
#'
#' ## Variables
#' nVariables(P)
#'
#' ##############################################################################
#' ## MORE EXAMPLES #############################################################
#' ##############################################################################
#' \donttest{
#' ## Generate a dataset
#' data <- data.frame(X1 = rnorm(100), X2 = rnorm(100), X3 = rnorm(100))
#'
#' ## Joint function
#' dim <- c(2,1,3)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#'
#' ## Variables
#' nVariables(P)
#' }
nVariables=function(P)
{
suppressWarnings({
Letters <- letters[c(24:26,1:23)]
string <- as.character(P)
t <- strsplit(string, split="-", fixed = T)[[1]]
t2 <- c()
for(i in 1:length(t)){
t1 <- strsplit(t[i], split="+", fixed = T, perl = FALSE, useBytes = FALSE)[[1]]
t2 <- c(t2,t1)
}
t <- unlist(sapply(1:length(t2), function(i) strsplit(t2[i], split="*", fixed = T, perl = FALSE, useBytes = FALSE)[[1]]))
variables <- t[is.na(as.numeric(t))]
t <- unlist(strsplit(variables, split="^", fixed = T, perl = FALSE, useBytes = FALSE))
t <- unlist(strsplit(t, split="(", fixed = T, perl = FALSE, useBytes = FALSE))
variables <- t[is.na(as.numeric(t))]
variables <- unique(variables)
variables <- Letters[Letters%in%variables]
})
return(variables)
}
#' Dimension of MoTBFs
#'
#' Get the dimension of \code{"motbf"} and \code{"jointmotbf"} densities.
#'
#' @param P An object of class \code{"motbf"} and subclass 'mop' or \code{"jointmotbf"}.
#' @return Dimension of the function.
#' @seealso \link{univMoTBF} and \link{jointMoTBF}
#' @export
#' @examples
#' ## 1. EXAMPLE
#' ## Data
#' X <- rnorm(2000)
#'
#' ## Univariate function
#' subclass <- "MOP"
#' f <- univMoTBF(X, POTENTIAL_TYPE = subclass)
#' dimensionFunction(f)
#'
#' ## 2. EXAMPLE
#' ## Dataset with 2 variables
#' X <- data.frame(rnorm(100), rnorm(100))
#'
#' ## Joint function
#' dim <- c(2,3)
#' param <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param)
#'
#' ## Dimension of the joint function
#' dimensionFunction(P)
#'
dimensionFunction <- function(P){
suppressWarnings({
string <- as.character(P)
nVar <- nVariables(P)
if(length(nVar)==1){
str <- strsplit(string, split="^", fixed=T)[[1]]
return(as.numeric(str[length(str)])+1)
}
dimensions <- c()
for(i in 1:length(nVar)){
elements <- strsplit(string, split=NULL)[[1]]
pos <- grep(nVar[i], elements)
nele <- elements[(pos[length(pos)]+1):(pos[length(pos)]+3)]
if(is.na(as.numeric(nele[3]))) nele <- nele[-3]
if(length(which(nele%in%"^"))!=0){
#nele[1]=="^")
pos <- as.numeric(nele[c(2,3)])
if(all(!is.na(pos))){
num <- c()
for(h in 2:length(nele)) num <- paste(num,nele[h], sep="")
dimensions=c(dimensions, as.numeric(num))
}else dimensions=c(dimensions, as.numeric(nele[2]))
} else dimensions=c(dimensions, 1)
}
return(dimensions+1)
})
}
#' Integration with MoTBFs
#'
#' Integrate a \code{"jointmotbf"} object over an non defined domain. It is able to
#' get the integral of a joint function over a set of variables or over all
#' the variables in the function.
#'
#' @param P A \code{"jointmotbf"} object.
#' @param var A \code{"character"} vector containing the name of the variables that will be integrated out.
#' Instead of the names, the position of the variables can be given.
#' By default it's \code{NULL} then all the variables are integrated out.
#' @return A multiintegral of a joint function of class \code{"jointmotbf"}.
#' @export
#' @examples
#' ## 1. EXAMPLE
#' ## Dataset with 2 variables
#' X <- data.frame(rnorm(100), rnorm(100))
#'
#' ## Joint function
#' dim <- c(2,3)
#' param <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param)
#'
#' ## Integral
#' integralJointMoTBF(P)
#' integralJointMoTBF(P, var="x")
#' integralJointMoTBF(P, var="y")
#'
#' ##############################################################################
#' ## MORE EXAMPLES #############################################################
#' ##############################################################################
#' \donttest{
#' ## Dataset with 3 variables
#' X <- data.frame(rnorm(50), rnorm(50), rnorm(50))
#'
#' ## Joint function
#' dim <- c(2,1,3)
#' param <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param)
#'
#' ## Integral
#' integralJointMoTBF(P)
#' integralJointMoTBF(P, var="x")
#' integralJointMoTBF(P, var="y")
#' integralJointMoTBF(P, var="z")
#' integralJointMoTBF(P, var=c("x","z"))
#' }
integralJointMoTBF <- function(P, var=NULL)
{
suppressWarnings({
Letters <- c(letters[24:26], letters[1:23])
if(is.numeric(var)) var <- Letters[var]
if(is.null(var)) var <- nVariables(P)
nVar <- nVariables(P)
for(h in 1:length(var)){
string <- as.character(P)
parameters <- coef(P)
f1 <- substr(string, 1, 1)
t <- strsplit(string, split="-", fixed = T)[[1]]
for(i in 1:length(t)) t[i] <- paste("-", t[i], sep="")##le volvemos a añadir el simbolo negativo
if(f1!=substr(t[1], 1, 1)) t[1] <- substr(t[1], 2, nchar(t[1]))
if(f1 == "-"){
t <- t[-1]
}
t2 <- c()
for(i in 1:(length(t))){
t1 <- strsplit(t[i], split="+", fixed = T, perl = FALSE, useBytes = FALSE)[[1]]
t2 <- c(t2,t1)
}
pos <- grep("e", t2)
if(length(pos)!=0){
for(i in pos){
t2[i] <- paste(t2[i],t2[i+1], sep="")
t2[i+1] <- NA
}
}
str <- t2[!is.na(t2)]
s <- strsplit(str, split=as.character(parameters), fixed=T)
for (i in 1:length(s)){
if(is.na(s[[i]][2])) s[[i]][2]= paste("*", var[h], sep="")
else{
t <- strsplit(s[[i]][2], split=NULL)[[1]]
if(!any(t%in%var[h])){
pVar <- which(nVar%in%var[h])
pVars <- which(nVar%in%t)
if(any(pVar<pVars)){
d <- which(nVar[pVars[which(pVar<pVars)[1]]]==t)
t <- append(t, paste("*", var[h], sep=""), d-2)
f <- c()
for(j in 1:length(t)) f <- paste(f, t[j], sep="")
s[[i]][2] <- f
}else{
s[[i]][2] <- paste(s[[i]][2], "*", var[h], sep="")
}
}else{
p <- which(t%in%var[h])
if(is.na(as.numeric(t[p+2]))){
t <- append(t, paste("^" , 2, sep=""),after=p)
parameters[i] <- parameters[i]/2
}else{
t[p+2]=as.numeric(t[p+2])+1
parameters[i] <- parameters[i]/as.numeric(t[p+2])
}
f <- c()
for(j in 1:length(t)) f <- paste(f, t[j], sep="")
s[[i]][2] <- f
}
}
}
if(any(parameters[2:length(parameters)]>=0)) parameters[parameters[2:length(parameters)]>=0]
parameters[2:length(parameters)][parameters[2:length(parameters)]>=0] =
paste("+",parameters[2:length(parameters)][parameters[2:length(parameters)]>=0], sep="")
s <- unlist(s)
s[which(s%in%"")]=parameters
str <- c()
for(i in 1:length(s)) str <- paste(str, s[i], sep="")
str <- noquote(str)
if(length(nVariables(str))>1) P <- jointmotbf(str)
else P <- motbf(str)
}
return(P)
})
}
#' Evaluation of joint MoTBFs
#'
#' Evaluates a \code{"jointmotbf"} object at a specific point.
#'
#' @param P A \code{"jointmotbf"} object.
#' @param values A list with the name of the variables equal to the values to be evaluated.
#' @return If all the variables in the equation are evaluated then a \code{"numeric"} value
#' is returned. Otherwise, an \code{"motbf"} object or a \code{"jointmotbf"} object is returned.
#' @export
#' @examples
#' #' ## 1. EXAMPLE
#' ## Dataset with 2 variables
#' X <- data.frame(rnorm(100), rexp(100))
#'
#' ## Joint function
#' dim <- c(3,3) # dim <- c(5,4)
#' param <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param)
#' P
#'
#' ## Evaluation
#' nVariables(P)
#' val <- list(x = -1.5, y = 3)
#' evalJointFunction(P, values = val)
#' val <- list(x = -1.5)
#' evalJointFunction(P, values = val)
#' val <- list(y = 3)
#' evalJointFunction(P, values = val)
#'
#' ##############################################################################
#' ## MORE EXAMPLES #############################################################
#' ##############################################################################
#' \donttest{
#' ## Dataset with 3 variables
#' X <- data.frame(rnorm(100), rexp(100), rnorm(100, 1))
#'
#' ## Joint function
#' dim <- c(2,1,3)
#' param <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param)
#' P
#'
#' ## Evaluation
#' nVariables(P)
#' val <- list(x = 0.8, y = -2.1, z = 1.2)
#' evalJointFunction(P, values = val)
#' val <- list(x = 0.8, z = 1.2)
#' evalJointFunction(P, values = val)
#' val <- list(y = -2.1)
#' evalJointFunction(P, values = val)
#' val <- list(y = -2.1)
#' evalJointFunction(P, values = val)
#' }
evalJointFunction <- function(P, values)
{
nVar <- nVariables(P)
if(is.motbf(P)) return(as.numeric(as.function(P)(values[[1]])))
string <- as.character(P)
parameters <- coef(P)
f1 <- substr(string, 1, 1)
t <- strsplit(string, split="-", fixed = T)[[1]]
for(i in 1:length(t)) t[i] <- paste("-", t[i], sep="")##le volvemos a añadir el simbolo negativo
if(f1!=substr(t[1], 1, 1)) t[1] <- substr(t[1], 2, nchar(t[1]))
if(f1 == "-"){
t <- t[-1]
}
t2 <- c()
for(i in 1:(length(t))){
t1 <- strsplit(t[i], split="+", fixed = T, perl = FALSE, useBytes = FALSE)[[1]]
t2 <- c(t2,t1)
}
pos <- grep("e", t2)
if(length(pos)!=0){
for(i in pos){
t2[i] <- paste(t2[i],t2[i+1], sep="")
t2[i+1] <- NA
}
}
str <- t2[!is.na(t2)]
s <- strsplit(str, split=as.character(parameters), fixed=T)
var <- names(values)
for(i in 1:length(s)){
for(j in 1:length(var)){
if(length(grep(var[j],s[[i]]))==0){
if(j==1){s[[i]][1]=parameters[i]; if(s[[i]][1]>0) s[[i]][1]=paste("+", s[[i]][1], sep="")}
}else{
if(j==1){s[[i]][1]=parameters[i]; if(s[[i]][1]>0) s[[i]][1]=paste("+", s[[i]][1], sep="")}
if(!is.na(s[[i]][2])){
s1=strsplit(s[[i]][2], split=NULL)[[1]]
if(any(s1%in%var[j])){
if(!is.na(s1[which(s1==var[j])+1])){
if(s1[which(s1==var[j])+1]=="^"){
exp=as.numeric(s1[which(s1==var[j])+2])
s[[i]][1]=as.numeric(s[[i]][1])*(values[[j]]^exp)
if(s[[i]][1]>0) s[[i]][1]=paste("+", s[[i]][1], sep="")
s1[which(s1==var[j])+1]=""
s1[which(s1==var[j])-1]=""
s1[which(s1==var[j])+2]=""
s1[s1%in%var[j]]=""
s2=c(); for(h in 1:length(s1)) s2=paste(s2, s1[h], sep="")
s[[i]][2]=s2
}else{
s[[i]][1]=as.numeric(s[[i]][1])*values[[j]]
if(s[[i]][1]>0) s[[i]][1]=paste("+", s[[i]][1], sep="")
s1[which(s1%in%var[j])-1]=""
s1[which(s1%in%var[j])]=""
s2=c(); for(h in 1:length(s1)) s2=paste(s2, s1[h], sep="")
s[[i]][2]=s2
}
}else{
s[[i]][1]=as.numeric(s[[i]][1])*values[[j]]
if(s[[i]][1]>0) s[[i]][1]=paste("+", s[[i]][1], sep="")
s1[which(s1%in%var[j])-1]=""
s1[which(s1%in%var[j])]=""
s2=c(); for(h in 1:length(s1)) s2=paste(s2, s1[h], sep="")
s[[i]][2]=s2
}
}
}
}
}
}
if(is.na(s[[1]][2])) s[[1]][2]=""
exp <- sapply(1:length(s), function(i) s[[i]][2])
exp <- unique(exp)
s <- unlist(s)
param <- c()
for(i in 1:length(exp)) param <- c(param,sum(as.numeric(s[which(s==exp[i])-1])))
sign <- ""
if(length(param)>1) for(i in 2:length(param)) if(param[i]>0) sign <- c(sign,"+") else sign <- c(sign,"")
str <- c()
for(i in 1:length(exp)) str <- paste(str, sign[i], param[i], exp[i], sep="")
if(all(nVar%in%names(values))){
message("The method 'as.function()' can be used. \n")
return(eval(parse(text=str)))
} else{
f <- noquote(str)
l <- length(nVariables(f))
if(l==1) f <- motbf(f)
if(l>1) f <- jointmotbf(f)
return(f)
}
}
#' Marginalization of MoTBFs
#'
#' Computes the marginal densities from a \code{"jointmotbf"}
#' object.
#'
#'@param P An object of class \code{"jointmotbf"}, i.e., the joint density function.
#'@param var The \code{"numeric"} position or the \code{"character"} name of the marginal variable.
#'@return The marginal of a \code{"jointmotbf"} function. The result is an object of class \code{"motbf"}.
#'@seealso \link{jointMoTBF} and \link{evalJointFunction}
#'@export
#'@examples
#' ## 1. EXAMPLE
#' ## Dataset with 2 variables
#' X <- data.frame(rnorm(100), rnorm(100))
#'
#' ## Joint function
#' dim <- c(4,3)
#' param <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param)
#' P
#'
#' ## Marginal
#' marginalJointMoTBF(P, var = "x")
#' marginalJointMoTBF(P, var = 2)
#'
#' ##############################################################################
#' ## MORE EXAMPLES #############################################################
#' ##############################################################################
#' \donttest{
#' ## Generate a dataset with 3 variables
#' data <- data.frame(rnorm(100), rnorm(100), rnorm(100))
#'
#' ## Joint function
#' dim <- c(2,1,3)
#' param <- parametersJointMoTBF(data, dimensions = dim)
#' P <- jointMoTBF(param)
#' nVariables(P)
#'
#' ## Marginal
#' marginalJointMoTBF(P, var="x")
#' marginalJointMoTBF(P, var="y")
#' marginalJointMoTBF(P, var="z")
#' }
marginalJointMoTBF <- function(P, var)
{
suppressWarnings({
Letters <- c(letters[24:26], letters[1:23])
if(is.numeric(var)) var <- Letters[var] else var <- var
varN <- nVariables(P); noVar=varN[!varN%in%var]
l <- length(noVar)
for(j in 1:l){
varN <- nVariables(P)
posnoVar <- which(varN%in%noVar[j])
posvar <- which(varN%in%var)
min <- lapply(posnoVar, function(i) P$Domain[1,i])
names(min) <- noVar[j]
max <- lapply(posnoVar, function(i) P$Domain[2,i])
names(max) <- noVar[j]
iP <- integralJointMoTBF(P, noVar[j])
minF <- evalJointFunction(iP, min)
maxF <- evalJointFunction(iP, max)
if(is.numeric(minF)&&is.numeric(maxF)){
P <- list(Function=(maxF-minF), Subclass="mop", Domain= P$Domain[,which(!varN%in%noVar[j])])
P <- motbf(P)
return(P)
}
## Min part
nVar <- nVariables(minF)
string <- as.character(minF)
f1 <- substr(string, 1, 1)
t <- strsplit(string, split="-", fixed = T)[[1]]
for(i in 1:length(t)) t[i] <- paste("-", t[i], sep="")
if(f1!=substr(t[1], 1, 1)) t[1] <- substr(t[1], 2, nchar(t[1]))
if(t[1]=="-") t=t[-1]
t2 <- c()
for(i in 1:(length(t))){
t1 <- strsplit(t[i], split="+", fixed = T, perl = FALSE, useBytes = FALSE)[[1]]
t2 <- c(t2,t1)
}
pos <- grep("e", t2)
if(length(pos)!=0){
for(i in pos){
t2[i] <- paste(t2[i],t2[i+1], sep="")
t2[i+1] <- NA
}
}
str <- t2[!is.na(t2)]
coefMIN <- coef(minF)
sMIN <- unlist(strsplit(str, split=coefMIN, fixed=T))
varseq=(unique(sMIN))[-1]
sMIN[which(sMIN=="")]=coefMIN
paramMIN <-c()
for(i in 1:length(varseq)){
pos <- which(sMIN%in%varseq[i])
paramMIN <- c(paramMIN, sum(as.numeric(sMIN[pos-1])))
sMIN <- sMIN[-c(pos, pos-1)]
}
paramMIN <- c(sum(as.numeric(sMIN)), paramMIN)
## Max part
nVar <- nVariables(maxF)
string <- as.character(maxF)
f1 <- substr(string, 1, 1)
t <- strsplit(string, split="-", fixed = T)[[1]]
for(i in 1:length(t)) t[i] <- paste("-", t[i], sep="")
if(f1!=substr(t[1], 1, 1)) t[1] <- substr(t[1], 2, nchar(t[1]))
if(t[1]=="-") t=t[-1]
t2 <- c()
for(i in 1:(length(t))){
t1 <- strsplit(t[i], split="+", fixed = T, perl = FALSE, useBytes = FALSE)[[1]]
t2 <- c(t2,t1)
}
pos <- grep("e", t2)
if(length(pos)!=0){
for(i in pos){
t2[i] <- paste(t2[i],t2[i+1], sep="")
t2[i+1] <- NA
}
}
str <- t2[!is.na(t2)]
coefMAX <- coef(maxF)
sMAX <- unlist(strsplit(str, split=coefMAX, fixed=T))
varseq=(unique(sMAX))[-1]
sMAX[which(sMAX=="")]=coefMAX
paramMAX <-c()
for(i in 1:length(varseq)){
pos <- which(sMAX%in%varseq[i])
paramMAX <- c(paramMAX, sum(as.numeric(sMAX[pos-1])))
sMAX <- sMAX[-c(pos, pos-1)]
}
paramMAX <- c(sum(as.numeric(sMAX)), paramMAX)
## New parameters
newparam <- paramMAX - paramMIN
sign=""
for(v in 2:length(newparam)) if(newparam[v]>0) sign=c(sign,"+") else sign=c(sign,"")
st <- newparam[1]
for(v in 2:length(newparam)) st <- paste(st, sign[v],newparam[v], varseq[v-1], sep="")
if(length(nVariables(st))>1){
P <- list(Function=noquote(st), Domain= P$Domain[,which(!varN%in%noVar[j])])
P <- jointmotbf(P)
}else{
P <- list(Function=noquote(st), Subclass="mop", Domain= P$Domain[,which(!varN%in%noVar[j])])
P <- motbf(P)
if(length(unique(coeffPol(P)))==1){
st <- paste(sum(coef(P)), ifelse(unique(coeffPol(P))==0, "",paste("*", var, ifelse(unique(coeffPol(P))==1, "", paste("^", unique(coeffPol(P)), sep="")), sep="")), sep="")
P <- list(Function=noquote(st), Subclass="mop", Domain= P$Domain)
P <- motbf(P)
}
}
}
return(P)
})
}
#' Bidimensional plots for \code{'jointmotbf'} objects
#'
#' PLot the perpective and the contour plots for joint MoTBF functions.
#'
#' @param x An object of class \code{'jointmotbf'}.
#' @param type A \code{"character"} string, either \emph{contour} or \emph{perspective}. It is set to \code{"contour"} by default.
#' @param ranges A \code{"numeric"} matrix containing the domain of the variables, by columns, which is used to specify the plotting range.
#' @param orientation A \code{"numeric"} vector indicating the perpective of the plot in degrees.
#' By default, it is set to \code{(5,-30)}.
#' @param data An object of class \code{"data.frame"} containing two columns only.
#' This argument is used to draw the points over the main plot. By default, it is set to \code{NULL}.
#' @param filled A logical argument; it is only used if \code{type = "contour"}.
#' is active. By default, it is \code{TRUE}, so filled contours are plotted.
#' @param ticktype A \code{"character"} string, either \emph{simple} or \emph{detailed}. By default, it is set to \code{"simple"},
#' which draws just an arrow parallel to the axis to indicate direction of increase.
#' In contrast, \code{"detailed"} draws normal ticks. This argument is only used in the \code{"perspective"} plot.
#' @param \dots Further arguments to be passed to \link{plot}.
#' @method plot jointmotbf
#' @return A plot of the joint MoTBF.
#' @seealso \link{jointMoTBF}
#' @export
#' @examples
#' ## 1 .EXAMPLE
#' ## Dataset
#' X <- data.frame(rnorm(500), rnorm(500))
#'
#' ## Joint function
#' dim <- c(3,3)
#' param1 <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param1)
#' P
#'
#' ## Plots
#' plot(P)
#' plot(P, type = "perspective", orientation = c(90,0))
#'
#' #############################################################################
#' ## MORE EXAMPLES ############################################################
#' #############################################################################
#' \donttest{
#' ## Dataset
#' X <- data.frame(rnorm(200,2), rexp(200, 1))
#'
#' ## Joint function
#' dim <- c(4,5)
#' param2 <- parametersJointMoTBF(X, dimensions = dim)
#' P <- jointMoTBF(param2)
#' P
#'
#' ## Plots
#' plot(P)
#' plot(P, filled = FALSE, data = X)
#' plot(P, type = "perspective", orientation = c(10,180))
#' }
plot.jointmotbf <- function(x, type="contour", ranges=NULL, orientation=c(5,-30),
data=NULL, filled=TRUE, ticktype="simple", ...)
{
opar <- par(no.readonly =TRUE)
on.exit(par(opar))
varname <- attr(x$Domain, "dimnames")[[2]]
if(length(nVariables(x))!=2) stop("It is not possible plotting a joint function with more than 2 variables.")
if(is.null(ranges)) ranges <- x$Domain
#{ranges <- x$Domain; ranges[1,] <- ranges[1,]+2E-1; ranges[2,] <- ranges[2,]-2E-1}
X <- seq(min(ranges[,1]), max(ranges[,1]), length = 20)
Y <- seq(min(ranges[,2]), max(ranges[,2]), length = 25)
Z <- outer(X, Y, as.function(x))
Z[Z<=0]=1.0E-10
if(type == "perspective"){
nrz <- nrow(Z)
ncz <- ncol(Z)
## Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("blue", "green") )
## Generate the desired number of colors from this palette
nbcol <- 100; color <- jet.colors(nbcol)
## Compute the z-value at the facet centres
zfacet <- Z[-1, -1] + Z[-1, -ncz] + Z[-nrz, -1] + Z[-nrz, -ncz]
## Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
persp(X, Y, Z, col = color[facetcol], phi = orientation[1], theta = orientation[2], ticktype=ticktype,
xlab = varname[1], ylab = varname[2], zlab = '',...)
}
if(type == "contour"){
if(!filled){
if(!is.null(data)) plot(data, xlab = varname[1], ylab = varname[2], xlim=ranges[,1], ylim=ranges[,2],...)
if(is.null(data)) plot(NULL, xlim=ranges[,1], ylim=ranges[,2], xlab="", ylab="")
contour(X,Y,Z, col=terrain.colors(15), xlab = varname[1], ylab = varname[2], add=T)
abline(h =round(ranges[,1])[1]:round(ranges[,1])[2] , v =round(ranges[,2])[1]:round(ranges[,2])[2],
col = "gray", lty = 2, lwd = 0.1)
}
if(filled){
zlim <- range(Z, finite=TRUE)
zlim[1] <- 0
nlevels <- 20
levels <- pretty(zlim, nlevels)
nlevels <- length(levels)
color <- colorRampPalette(c("#FFFFD9", "#EDF8B1", "#C7E9B4", "#7FCDBB", "#41B6C4", "#1D91C0", "#225EA8", "#253494", "#081D58"))
filled.contour(X,Y,Z,nlevels=nlevels,levels=levels,
las=1,
col=color(nlevels),
#col=rainbow(nlevels, alpha=0.8),
plot.title = {title(xlab = varname[1], ylab = varname[2], cex.lab = 1.5)},...
)
if(!is.null(data)){
mar.orig <- par("mar")
w <- (3 + mar.orig[2]) * par("csi") * 2.54
layout(matrix(c(2, 1), ncol = 2), widths = c(1, lcm(w)))
points(data)
par(mfrow=c(1,1))
}
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.