R/normalizeGaussian_severalstations.R
In ecor/RMRAINGEN: RMRAINGEN (R Multi-site RAINfall GENeretor): a package to generate daily time series of rainfall from monthly mean values
# TODO: Add comment
#
# Author: ecor
###############################################################################
NULL
#' Wilsks Correctio for Precipitation Marginal Gausianization
#'
#' Wrapper Gaussianization of \code{\link{normalizeGaussian_severalstations}}
#'
#'
#'
#' @param x see \code{\link{normalizeGaussian_severalstations}}, e. g. precipitation depth values
#' @param data see \code{\link{normalizeGaussian_severalstations}}
#'
#' @param prec_tolerance tolerance used for precipitation value
#' @param iterations number of iteration proposed for 'Wilks Gaussianization'
#' @param force.precipitation.value logical value. If it is \code{TRUE} (Default) gaussianized values corresponding to precipitation days are forced to fit the observed precipitation values.
#' @param seed seed used for random generation.
#'
#' @param valmin,tolerance see \code{\link{CCGamma}}
#' @param ... further arguments for \code{\link{normalizeGaussian_severalstations}}
#'
#'
#' @export
#'
#'
#'
#'
#'
#'
#' @seealso \code{\link{normalizeGaussian_severalstations}},\code{\link{CCGamma}}
#'
#' @examples
#'
#'
#'library(RMRAINGEN)
#'
#'
#'data(trentino)
#'
#'year_min <- 1961
#'year_max <- 1990
#'origin <- paste(year_min,1,1,sep="-")
#'
#'period <- PRECIPITATION$year>=year_min & PRECIPITATION$year<=year_max
#'station <- names(PRECIPITATION)[!(names(PRECIPITATION) %in% c("day","month","year"))]
#'prec_mes <- PRECIPITATION[period,station]
#'
#'
#'
#' ## removing nonworking stations (e.g. time series with NA)
#'accepted <- array(TRUE,length(names(prec_mes)))
#' names(accepted) <- names(prec_mes)
#'for (it in names(prec_mes)) {
#' accepted[it] <- (length(which(!is.na(prec_mes[,it])))==length(prec_mes[,it]))
#'}
#'
#'prec_mes <- prec_mes[,accepted]
#'valmin <- 0.5
#'prec_mes_gaussWilks <- WilksGaussianization(x=prec_mes, data=prec_mes,valmin=valmin,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin,force.precipitation.value="both")
#' prec_mes_gauss <- normalizeGaussian_severalstations(x=prec_mes, data=prec_mes,step=0,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin)
#### prec_mes_gauss_nullstep <- normalizeGaussian_severalstations(x=prec_mes, data=prec_mes,step=NULL,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin)
# CCGamma <- CCGamma(data=prec_mes, lag = 0,valmin=valmin,only.matrix=TRUE,tolerance=0.001)
#' prec_mes_ginv <- list()
#' prec_mes_ginv$unforced <- normalizeGaussian_severalstations(x=prec_mes_gaussWilks$unforced, data=prec_mes,step=0,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin,inverse=TRUE)
#' prec_mes_ginv$forced <- normalizeGaussian_severalstations(x=prec_mes_gaussWilks$forced, data=prec_mes,step=0,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin,inverse=TRUE)
#'
#' str(prec_mes_ginv)
#'
#' plot(prec_mes[,1],prec_mes_ginv$unforced[,1])
#' plot(prec_mes[,19],prec_mes_ginv$unforced[,19])
#'
#' CCGamma <- CCGamma(data=prec_mes, lag = 0,valmin=valmin,only.matrix=TRUE,tolerance=0.001)
#'
#'
# plot(prec_mes[,5],prec_mes_ginv[,5])
# plot(prec_mes_gaussWilks[,1],prec_mes_gauss[,1])
# plot(prec_mes_gaussWilks[,2],prec_mes_gauss[,2])
# plot(prec_mes_gaussWilks[,19],prec_mes_gauss[,19])
# plot(prec_mes_gaussWilks[,10],prec_mes_gauss[,10])
#' plot(cor(prec_mes_gaussWilks$unforced),cor(prec_mes_gaussWilks$forced))
#' abline(0,1)
#'
#'
#' VARselect(prec_mes_gaussWilks$forced)
#'
#'
#' # u <- apply(prec_mes_ginv$forced,2,rank)/(nrow(prec_mes_ginv$forced)+1)
#' #copula <- normalCopula(dim=ncol(u), disp = "un",param=P2p(CCGamma))
#' #out <- fitCopula(copula, data=u, method = "ml" )
#' # start = NULL, lower = NULL, upper = NULL,
#' # optim.method = "BFGS", optim.control = list(maxit=1000),
#' # estimate.variance = TRUE, hideWarnings = TRUE)
# #### c("mpl", "ml", "itau", "irho"),
#
#'
# ### ,step = valmin, prec = 10^-4, type = 3,
# extremes = TRUE, sample = "monthly", origin_x = origin, origin_data = origin)
#
#
#
#'
#
#
WilksGaussianization <- function (x,data=x,gauss=NULL,valmin=0.5,tolerance=0.001,prec_tolerance=c(0.05,2.5),iterations=20,force.precipitation.value=TRUE,seed=1234,shuffle=list(e1=1:ncol(x)),p=1,args_var=NULL,...) {
if (!is.list(shuffle) | is.data.frame(shuffle)) shuffle <- list(v1=shuffle)
if (length(shuffle)>1) {
cc <- length(shuffle)
for (i in 1:(cc-1)) {
gauss <- WilksGaussianization(x=x,data=data,gauss=gauss,valmin=valmin,tolerance=tolerance,iterations=iterations,force.precipitation.value=TRUE,seed=seed,shuffle=shuffle[i],...)
}
shuffle <- shuffle[cc]
### TO FINISH .....
}
if (!is.list(shuffle)) shuffle <- list(v1=shuffle)
set.seed(seed)
lag=0
names <- names(x)
x <- x[,shuffle[[1]]]
x[x<valmin] <- 0
data[data<valmin] <- 0
if (is.null(gauss)) {
gauss <- normalizeGaussian_severalstations(x=x,data=data,step=0,...)
} else if (!all.equal(names,names(gauss))) {
stop("Error in WilksGaussianization: different names between gauss and x!!!")
}
ccgamma <- CCGamma(data=x, lag = lag,valmin=valmin,only.matrix=TRUE,tolerance=tolerance)
### change ccgamma
ccgamma_block <- CCGammaToBlockmatrix(data=x,lag=0,valmin=valmin,tolerance=tolerance,p=p+1)
CCGamma0 <- as.blockmatrix(ccgamma_block[1:p,1:p],nrow=p,ncol=p)
CCGamma1 <- as.blockmatrix(ccgamma_block[(1:p),(1:p)+1],nrow=p,ncol=p)
CCGamma_1t <- t(CCGamma1)
CCGamma_0t <- t(CCGamma0)
A <- t(solve(CCGamma_0t,CCGamma_1t,symm.precond=TRUE)) ### A blockmatrix from Yule-Walker Equation
eigenvs <- abs(eigen(as.matrix(A))$values)
stable <- eigenvs<=1
stable <- length(which(stable))==length(stable)
ccgamma <- ccgamma_block[1,1]
### END change ccgamma
chol_ccg <- t(as.matrix(chol(ccgamma)))
out <- gauss
for (iter in 1:iterations) {
##print(iter)
cor <- nearPD(cor(gauss),cor=TRUE)$mat
chol_cor <- t(as.matrix(chol(cor)))
inv_chol_cor <- solve(chol_cor)
for (r in 1:nrow(out)) {
temp <- as.vector(gauss[r,])
temp <- t(as.matrix(temp))
xv <- inv_chol_cor %*% temp
out[r,] <- chol_ccg %*% as.matrix(xv)
}
message <- sprintf("Iteration: %d on %d",iter,iterations)
x_rec <- normalizeGaussian_severalstations(x=out, data=data,step=0,inverse=TRUE,...)
if (length(prec_tolerance)==1) prec_tolerance[2] <- 50*prec_tolerance[2]
cond <- abs(x-x_rec)<prec_tolerance[1]
condhigh <- abs(x-x_rec)<prec_tolerance[2]
points <- length(cond)
adjusted <- length(which(!cond))
adjustedhigh <- length(which((!condhigh) & (x>valmin)))
message2 <- sprintf("Points adjusted: %d (%d) on %d",adjusted,adjustedhigh,points)
print(message)
print(message2)
gauss[cond] <- out[cond]
#### INSERT AUTOREGRESSION WITH VARs
if (p>0) {
###
args_var <- c(list(y=out,p=p),args_var)
var <- do.call("VAR",args=args_var)
res <- as.data.frame(residuals(var))
str(res)
stop()
###
if (!stable) {
print(eigenves)
stop("Error A eigenvalues greater than 1 in modulus!!!")
}
xprev <-
print(A)
print(eigen(as.matrix(A))$values)
### fare un generate.blockmatrix ....
stop()
}
}
if (force.precipitation.value==TRUE) {
cond <- x>valmin
out[cond] <- gauss[cond]
out <- out[,names]
} else if (force.precipitation.value=="both") {
gauss_unforced <- out
out <- list()
out$forced <- gauss_unforced
out$unforced <- gauss_unforced
cond <- x>valmin
out$forced[cond] <- gauss[cond]
out$forced <- out$forced[,names]
out$unforced <- out$unforced[,names]
}
### }
return(out)
}
NULL
#'
#'
#'
#'
#'
#'
generate.blockmatrix <- function(x,xprev=NULL,noise=noise,n=10,x.noise.gen=NULL,...) {
p <- nrow(x)-1
if (p>1) count <- 1:p
if (!is.null(noise)) n <- nrow(noise)
if (is.null(noi))
if (p>1) {
noisep <- NULL
for (i in 0:(p-1)) {
count <- (p+1):n-l
noisep <- cbind(noisep,noise[count,])
}
print(noisep)
###### TO DOOO
### GENERARE DA P+1 IN POI ......
}
}
ecor/RMRAINGEN documentation built on May 15, 2019, 8:53 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.
# TODO: Add comment
#
# Author: ecor
###############################################################################
NULL
#' Wilsks Correctio for Precipitation Marginal Gausianization
#'
#' Wrapper Gaussianization of \code{\link{normalizeGaussian_severalstations}}
#'
#'
#'
#' @param x see \code{\link{normalizeGaussian_severalstations}}, e. g. precipitation depth values
#' @param data see \code{\link{normalizeGaussian_severalstations}}
#'
#' @param prec_tolerance tolerance used for precipitation value
#' @param iterations number of iteration proposed for 'Wilks Gaussianization'
#' @param force.precipitation.value logical value. If it is \code{TRUE} (Default) gaussianized values corresponding to precipitation days are forced to fit the observed precipitation values.
#' @param seed seed used for random generation.
#'
#' @param valmin,tolerance see \code{\link{CCGamma}}
#' @param ... further arguments for \code{\link{normalizeGaussian_severalstations}}
#'
#'
#' @export
#'
#'
#'
#'
#'
#'
#' @seealso \code{\link{normalizeGaussian_severalstations}},\code{\link{CCGamma}}
#'
#' @examples
#'
#'
#'library(RMRAINGEN)
#'
#'
#'data(trentino)
#'
#'year_min <- 1961
#'year_max <- 1990
#'origin <- paste(year_min,1,1,sep="-")
#'
#'period <- PRECIPITATION$year>=year_min & PRECIPITATION$year<=year_max
#'station <- names(PRECIPITATION)[!(names(PRECIPITATION) %in% c("day","month","year"))]
#'prec_mes <- PRECIPITATION[period,station]
#'
#'
#'
#' ## removing nonworking stations (e.g. time series with NA)
#'accepted <- array(TRUE,length(names(prec_mes)))
#' names(accepted) <- names(prec_mes)
#'for (it in names(prec_mes)) {
#' accepted[it] <- (length(which(!is.na(prec_mes[,it])))==length(prec_mes[,it]))
#'}
#'
#'prec_mes <- prec_mes[,accepted]
#'valmin <- 0.5
#'prec_mes_gaussWilks <- WilksGaussianization(x=prec_mes, data=prec_mes,valmin=valmin,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin,force.precipitation.value="both")
#' prec_mes_gauss <- normalizeGaussian_severalstations(x=prec_mes, data=prec_mes,step=0,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin)
#### prec_mes_gauss_nullstep <- normalizeGaussian_severalstations(x=prec_mes, data=prec_mes,step=NULL,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin)
# CCGamma <- CCGamma(data=prec_mes, lag = 0,valmin=valmin,only.matrix=TRUE,tolerance=0.001)
#' prec_mes_ginv <- list()
#' prec_mes_ginv$unforced <- normalizeGaussian_severalstations(x=prec_mes_gaussWilks$unforced, data=prec_mes,step=0,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin,inverse=TRUE)
#' prec_mes_ginv$forced <- normalizeGaussian_severalstations(x=prec_mes_gaussWilks$forced, data=prec_mes,step=0,sample="monthly",extremes=TRUE,origin_x = origin, origin_data = origin,inverse=TRUE)
#'
#' str(prec_mes_ginv)
#'
#' plot(prec_mes[,1],prec_mes_ginv$unforced[,1])
#' plot(prec_mes[,19],prec_mes_ginv$unforced[,19])
#'
#' CCGamma <- CCGamma(data=prec_mes, lag = 0,valmin=valmin,only.matrix=TRUE,tolerance=0.001)
#'
#'
# plot(prec_mes[,5],prec_mes_ginv[,5])
# plot(prec_mes_gaussWilks[,1],prec_mes_gauss[,1])
# plot(prec_mes_gaussWilks[,2],prec_mes_gauss[,2])
# plot(prec_mes_gaussWilks[,19],prec_mes_gauss[,19])
# plot(prec_mes_gaussWilks[,10],prec_mes_gauss[,10])
#' plot(cor(prec_mes_gaussWilks$unforced),cor(prec_mes_gaussWilks$forced))
#' abline(0,1)
#'
#'
#' VARselect(prec_mes_gaussWilks$forced)
#'
#'
#' # u <- apply(prec_mes_ginv$forced,2,rank)/(nrow(prec_mes_ginv$forced)+1)
#' #copula <- normalCopula(dim=ncol(u), disp = "un",param=P2p(CCGamma))
#' #out <- fitCopula(copula, data=u, method = "ml" )
#' # start = NULL, lower = NULL, upper = NULL,
#' # optim.method = "BFGS", optim.control = list(maxit=1000),
#' # estimate.variance = TRUE, hideWarnings = TRUE)
# #### c("mpl", "ml", "itau", "irho"),
#
#'
# ### ,step = valmin, prec = 10^-4, type = 3,
# extremes = TRUE, sample = "monthly", origin_x = origin, origin_data = origin)
#
#
#
#'
#
#
WilksGaussianization <- function (x,data=x,gauss=NULL,valmin=0.5,tolerance=0.001,prec_tolerance=c(0.05,2.5),iterations=20,force.precipitation.value=TRUE,seed=1234,shuffle=list(e1=1:ncol(x)),p=1,args_var=NULL,...) {
if (!is.list(shuffle) | is.data.frame(shuffle)) shuffle <- list(v1=shuffle)
if (length(shuffle)>1) {
cc <- length(shuffle)
for (i in 1:(cc-1)) {
gauss <- WilksGaussianization(x=x,data=data,gauss=gauss,valmin=valmin,tolerance=tolerance,iterations=iterations,force.precipitation.value=TRUE,seed=seed,shuffle=shuffle[i],...)
}
shuffle <- shuffle[cc]
### TO FINISH .....
}
if (!is.list(shuffle)) shuffle <- list(v1=shuffle)
set.seed(seed)
lag=0
names <- names(x)
x <- x[,shuffle[[1]]]
x[x<valmin] <- 0
data[data<valmin] <- 0
if (is.null(gauss)) {
gauss <- normalizeGaussian_severalstations(x=x,data=data,step=0,...)
} else if (!all.equal(names,names(gauss))) {
stop("Error in WilksGaussianization: different names between gauss and x!!!")
}
ccgamma <- CCGamma(data=x, lag = lag,valmin=valmin,only.matrix=TRUE,tolerance=tolerance)
### change ccgamma
ccgamma_block <- CCGammaToBlockmatrix(data=x,lag=0,valmin=valmin,tolerance=tolerance,p=p+1)
CCGamma0 <- as.blockmatrix(ccgamma_block[1:p,1:p],nrow=p,ncol=p)
CCGamma1 <- as.blockmatrix(ccgamma_block[(1:p),(1:p)+1],nrow=p,ncol=p)
CCGamma_1t <- t(CCGamma1)
CCGamma_0t <- t(CCGamma0)
A <- t(solve(CCGamma_0t,CCGamma_1t,symm.precond=TRUE)) ### A blockmatrix from Yule-Walker Equation
eigenvs <- abs(eigen(as.matrix(A))$values)
stable <- eigenvs<=1
stable <- length(which(stable))==length(stable)
ccgamma <- ccgamma_block[1,1]
### END change ccgamma
chol_ccg <- t(as.matrix(chol(ccgamma)))
out <- gauss
for (iter in 1:iterations) {
##print(iter)
cor <- nearPD(cor(gauss),cor=TRUE)$mat
chol_cor <- t(as.matrix(chol(cor)))
inv_chol_cor <- solve(chol_cor)
for (r in 1:nrow(out)) {
temp <- as.vector(gauss[r,])
temp <- t(as.matrix(temp))
xv <- inv_chol_cor %*% temp
out[r,] <- chol_ccg %*% as.matrix(xv)
}
message <- sprintf("Iteration: %d on %d",iter,iterations)
x_rec <- normalizeGaussian_severalstations(x=out, data=data,step=0,inverse=TRUE,...)
if (length(prec_tolerance)==1) prec_tolerance[2] <- 50*prec_tolerance[2]
cond <- abs(x-x_rec)<prec_tolerance[1]
condhigh <- abs(x-x_rec)<prec_tolerance[2]
points <- length(cond)
adjusted <- length(which(!cond))
adjustedhigh <- length(which((!condhigh) & (x>valmin)))
message2 <- sprintf("Points adjusted: %d (%d) on %d",adjusted,adjustedhigh,points)
print(message)
print(message2)
gauss[cond] <- out[cond]
#### INSERT AUTOREGRESSION WITH VARs
if (p>0) {
###
args_var <- c(list(y=out,p=p),args_var)
var <- do.call("VAR",args=args_var)
res <- as.data.frame(residuals(var))
str(res)
stop()
###
if (!stable) {
print(eigenves)
stop("Error A eigenvalues greater than 1 in modulus!!!")
}
xprev <-
print(A)
print(eigen(as.matrix(A))$values)
### fare un generate.blockmatrix ....
stop()
}
}
if (force.precipitation.value==TRUE) {
cond <- x>valmin
out[cond] <- gauss[cond]
out <- out[,names]
} else if (force.precipitation.value=="both") {
gauss_unforced <- out
out <- list()
out$forced <- gauss_unforced
out$unforced <- gauss_unforced
cond <- x>valmin
out$forced[cond] <- gauss[cond]
out$forced <- out$forced[,names]
out$unforced <- out$unforced[,names]
}
### }
return(out)
}
NULL
#'
#'
#'
#'
#'
#'
generate.blockmatrix <- function(x,xprev=NULL,noise=noise,n=10,x.noise.gen=NULL,...) {
p <- nrow(x)-1
if (p>1) count <- 1:p
if (!is.null(noise)) n <- nrow(noise)
if (is.null(noi))
if (p>1) {
noisep <- NULL
for (i in 0:(p-1)) {
count <- (p+1):n-l
noisep <- cbind(noisep,noise[count,])
}
print(noisep)
###### TO DOOO
### GENERARE DA P+1 IN POI ......
}
}
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.