Nothing
#(tres <- twUtestF(combinePops,"test.twDEMC.apply"))
twDEMCApply <- function(
### applies function to each chain in $parms of an twDEMC
vtwdemc,
fun,
### function(x,name) applied to each variable in twDEMC
otherMargin = 1,
### margin for that each index the an array slice is provided as x to fun, defaults to 1: Vars (2: generations, 3:chains)
...
){
##details<<
## the arguments of fun are \describe{
## \item{x}{a matrix corresponging to margin (default with rows correponding to generations and cols corresponding to) }
## \item{name}{the name of the remaining margin, or the number of the item for otherMargin=2,3}}
## the return value must be of same size as x
if( length(otherMargin)!= 1)
stop(paste("twDEMC.applyParms: margin must be of type c(i,j) with i!=j and i,j in 1:3"))
res <- vtwdemc
#for( iVarName in unlist(dimnames(vtwdemc$parms)[otherMargin]) ){
.ind <- 1:(dim(vtwdemc$parms)[otherMargin])
.names <- unlist(dimnames(vtwdemc$parms)[otherMargin])
if( is.null(.names) ) .names=.ind
for( iVar in .ind ){
#aperm(res$parms,c(otherMargin,margin)) <- fun( aperm(vtwdemc$parms,c(otherMargin,margin))[iVar,,],iVar,...)
iVarName = .names[iVar]
if( otherMargin==1) res$parms[iVar,,] <- fun( vtwdemc$parms[iVar,,],iVarName,...)
else if( otherMargin==2) res$parms[,iVar,] <- fun( vtwdemc$parms[,iVar,],iVarName,...)
else if( otherMargin==3) res$parms[,,iVar] <- fun( vtwdemc$parms[,,iVar],iVarName,...)
}
res
}
#mtrace(twDEMCApply)
test.twDEMCApply <- function(){
.varDistr <- c(a="lognorm",b=NULL)
.tmp <- twDEMCApply(twdemcEx1,function(x,varname){ transOrigPopt.default(x,.varDistr[varname])})
checkEquals(dim(twdemcEx1$parms),dim(.tmp$parms))
checkEquals(.tmp$parms["a",,], exp(twdemcEx1$parms["a",,]) )
checkEquals(.tmp$parms["b",,], (twdemcEx1$parms["b",,]) )
.tmp2 <- twDEMCApply(twdemcEx1,otherMargin=2, fun=function(x,varname){ x+1 })
checkEquals(.tmp2$parms[,,], (twdemcEx1$parms[,,]+1) )
}
twDEMC <- function(
### Differential Evolution Markov Chain
Zinit,...
){
##details<<
## calls correct initializer of \code{\link{twDEMCBlockInt}} depending on the class of Zinit
UseMethod("twDEMC"
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{twDEMC.default}}
## \code{\link{twDEMC.twDEMC}}
)}
#gamma <- popt.tmp; gamma0 <- tmp.m
rotateInd <- function(gamma, gamma0, A){
# tansInd
#
# transpose parameter vector gamma to independent parameters xi
# gamma parameters to transpose
# gamma0 approximate mean of the posterior
# A matrix square root, usually obtained by chol(cov), where cov is the covariance of the posterior
gd <- gamma - gamma0
xi <- backsolve( A, gd )
xi
}
rotateIndBack <- function(xi,gamma0,A){
# tansIndBack
#
# transpose independent parameters xi back to original parameters (inverse of rotateInd)
# xi parameters to transpose back
# gamma0 approximate mean of the posterior
# A upper diagonal matrix square root, usually obtained by chol(cov), where cov is the covariance of the posterior
gamma = A %*% xi + gamma0
gamma
}
#------------ rotating parameter vectors
tmp.f <- function(){
?chol
tmp.cov <- cov( pss[,-1] ) #sample covariance
tmp.m <- colMeans(pss[,-1]) #sample mean
A = chol(tmp.cov) #upper t
xi <- transInd(popt.tmp,tmp.m, A)
xi2 <- solve(A) %*% gd
xi - xi2
tmp <- rotateIndBack(xi,tmp.m, A)
popt.tmp - tmp
mtrace(transInd)
pssXi <- t(transInd( t(pss[sample(1:nrow(pss),20),-1]), tmp.m, A ))
cor(pssXi)
corOrdered(pssXi) #why still highly dependent?
}
checkConvergenceFalse <- function(res, addArgs=list() ){ FALSE }
checkConvergenceGelman <- function(res, addArgs=list() ){
# checkConvergenceGelman
#
# gelman RHat criterion applied to second half of each chain of a population
# res see return value of DEMCzsp ($parms (d x nStep x nChain) )
# value: TRUE/FALSE
# see Gelman04 (twutz:Gelman04_3#Inference_and_assessing_convergence)
# addArgs
# $burninFrac fraction of the chain to be discarded (default 0.5)
# $rHatMin rHat criterion, returns TRUE if smaller than this
if( is.null(addArgs$burninFrac) ) addArgs$burninFrac=0.5
if( is.null(addArgs$rHatMin) ) addArgs$rHatMin = 1.1
l <- dim(res$parms)[2]
res2 <- res$parms[ ,(l*addArgs$burninFrac):l, ] # second half of all the chains
n <- dim(res2)[2] # number of steps
m <- dim(res2)[3] # number of chains
rl <- sapply( 1:dim(res2)[1], function(vn){ #over all variables (estimands)
muyj <- sapply(1:m, function(j){ mean(res2[vn,,j]) }) # mean of each chain
muy <- mean(muyj) # mean across chains
B <- n/(m-1)*sum((muyj-muy)^2) # between chain variance
s2j <- sapply(1:m, function(j){sum((res2[vn,,j]-muyj[j])^2) })/(n-1) #squared sums of each chain
W <- sum(s2j)/m # within chain variance
VarPlus=(n-1)/n*W + 1/n*B
#c( B=B, W=W, VarPlus=VarPlus, Rhat=round(sqrt(VarPlus/W),3), effSize=min(n*m,n*m*VarPlus/B) )
c( Rhat2=VarPlus/W )
})
all(rl <= (addArgs$rHatMin)^2 )
}
#mtrace( checkConvergenceGelman )
.testVectorApply <- function(){
tmp.f <- function(x){ 1:4 }
x <- matrix(1:3, nrow=3, ncol=5)
(.tmp <- apply( x, 2, tmp.f )) #each vector becomes a column
tmp.f1 <- function(x){ 1 }
(.tmp <- apply( x, 2, tmp.f1 )) #each scaler becomes component of a vector
}
.tmp.f <- function(
doLocalAdaptation=FALSE
### if TRUE, then states for differential proposal are sampled with heigher weight from current chain and recent past
, doAdaptT=FALSE
### if TRUE, Temperature will be adapted to acceptance rate
){
}
.doDEMCSteps <- function( iGenTStart, ctrl,
X, xStep,
fDiscrProp, argsFDiscrProp,
logDenCompX, argsFLogDen, fLogDenScale, debugSequential, ...
){
for( iGenT in (1:ctrl$thin) ){
iGen = iGenTStart+iGenT
xProp = X + xStep[,,iGenT]
# if necessary, discretize proposal
if( is.function(fDiscrProp))
xProp = do.call(fDiscrProp,xProp,argsFDiscrProp, quote=TRUE)
##details<<
## in order to support a two-level Metropolis desition
#twCalcLogDenPar expects parameters in columns, need transpose
.resLogDenPar <- twCalcLogDenPar(fLogDen=fLogDen, xProp=t(xProp), logDenCompX=logDenCompX, argsFLogDen=argsFLogDen, fLogDenScale=fLogDenScale, debugSequential=debugSequential, ...)
logfitness_x_prop <- .resLogDenPar$logDen
logDenCompProp <- .resLogDenPar$logDenComp
TcurStep = TstepFixed[iGen,]
#logr = (logfitness_x_prop+rExtra - logfitness_X)/Tstep[iGen]
logr = (logfitness_x_prop+rExtra[,iGenT] - logfitness_X) / rep(TcurStep, each=nChainPop)
# print(c(logfitness_X[i], logfitness_x_prop ,logr,rExtra))
#Metropolis step for each chain
for (i in iseq){ #body must affect outer variables: cannot use apply instead of for
if ( is.finite(logr[i]) & (logr[i]) > log(runif(1)) ){
acceptN[i] = acceptN[i]+1
X[,i] = xProp[,i]
logfitness_X[i] = logfitness_x_prop[i]
if( 0<length(logDenCompX) )
logDenCompX[i,] = logDenCompProp[i,] #here chains are rows
acceptWindow[ acceptPos, i ] <- TRUE
}
} #Metropolis step for each chain
}# iGenT within thinning interval
resDo <- list( acceptN=acceptN, X, logfitness_X, logDenCompX, acceptWindow )
### list with components \describe{
### \item{acceptN}{vector number of accepted steps per chain in thinning interval}
### \item{X}{matrix current position for each chain (column?)}
### \item{logfitness_X}{vector current logDen of chains}
### \item{logDenCompX}{matrix: result of fLogDen for last accepted state per chain}
### \item{acceptWindow}{}
}
.generateXPropThin <- function(
nPops,
Z,mZ,
ctrl,
...
){
d <- as.list(structure(dim(Z),names=c("parms","gen","chains")))
d$gen <- mZ
d$steps <- ctrl$thin
nChainPop = d$chains %/% nPops
##details<<
## Random states for chains for difference vector generation are within subsets of populations.
## This allows simulating independent population of chains.
## The acceptance rate may differ amonng populations. Hence, the set of previous generations to
## randomly select from also differs between poplations.
# integer array (thinSteps*nChain*3) sample of chains, within populations
rrChains <- abind( lapply( 1:nPops, function(iPop){
sChains <- ((iPop-1)*nChainPop+1):(iPop*nChainPop)
rrChainsPop <- array( sample(sChains, d$steps*nChainPop*3, replace = TRUE), dim=list(gen=d$steps,chain=nChainPop,iPop=3) )
}), along=2 )
# integer array same dimension of sample of generations
rrGen <- abind( lapply( 1:nPops, function(iPop){
sGens <- (mZ-g[iPop]+1):mZ
rrGenPop <- array( sample(sGens, d$steps*nChainPop*3, replace = TRUE), dim=list(gen=d$steps,chain=nChainPop,i=3) )
}), along=2 )
}
calcCorrParsDf <- function(parms, normpopt){
# calcCorrParsDf
#
# caluclates the values of correlated parameters on transformed, i.e. normal, scale
# based on dataframe parms$corrParsDf, which has columns y,x (char), a, b
# normpcorr[y] = a*normpopt[x] +b
# which is the result of linApprox in R/sensitivity.R of this package
tmp5 <- parms$corrParsDf
if( is.data.frame(tmp5) && (nrow(tmp5) > 0)){
tmp <- sapply( 1:nrow(tmp5), function(i){
row = tmp5[i,]
row$a * normpopt[ row$x ] + row$b
})
names(tmp) <- tmp5$y
tmp
}else c()
}
calcCorrParsDfMat <- function(parms, normpopt){
# calcCorrParsDfMat
#
# same as calcCorrParsDf but with normpopt specified as a matrix with mutliple rows
# caluclates the values of correlated parameters on transformed, i.e. normal, scale
# based on dataframe parms$corrParsDf, which has columns y,x (char), a, b
# normpcorr[y] = a*normpopt[x] +b
# which is the result of linApprox in R/sensitivity.R of this package
tmp5 <- parms$corrParsDf
if( is.data.frame(tmp5) && (nrow(tmp5) > 0)){
tmp <- sapply( 1:nrow(tmp5), function(i){
row = tmp5[i,]
row$a * normpopt[, row$x ] + row$b
})
colnames(tmp) <- tmp5$y
tmp
}else c()
}
.test.ggstat <- function(){
str(diamonds)
p1 <- ggplot(diamonds, aes(x=color,y=carat))
p2 <- p1 + geom_point()
p1 + stat_summary(aes(colour="mean"), fun.y=mean, geom="point")
+ scale_colour_hue("Mean")
x <- rnorm(100)
mtrace(dnorm)
p1 <- ggplot(x)
p1 + geom_
qplot(x, geom="density") + stat_function(fun = dnorm, colour="red")
}
R.methodsS3::setMethodS3("transOrigPopt","twDEMC", function(
### Applies \code{\link{transOrigPopt.default}} to each column of parameters in \code{vtwdemc}.
vtwdemc,
### list of class twDEMC with $parms in transformed scale
poptDistr=eval(parse(text="parDistr$trans")),
### character vector of kind of transformation ("lognorm"/"logitnorm") for each column of normpopt
...
){
##seealso<<
## \code{\link{transOrigPopt.default}}
res <- vtwdemc
normpopt <- vtwdemc$parms
##details<<
## Either \code{poptDistr} has names for each column name,
## or \code{poptDistr} has the same length as \code{colnames(normpopt)}.
if( !is.null(names(poptDistr)) )
if( all(rownames(normpopt) %in% names(poptDistr)) )
poptDistr=poptDistr[ rownames(normpopt) ]
if( length(poptDistr) != nrow(normpopt) )
stop("poptDistr must have entries for each row or must have the same length as nrow(normpopt$parms).")
for( vi in seq(along.with=rownames(normpopt)) ){
res$parms[vi,,] <- transOrigPopt.default( normpopt[vi,,], poptDistr[vi] )
if( 0 < length(vtwdemc$Y) )
res$Y[ rownames(normpopt)[vi] ,,] <- transOrigPopt.default( res$Y[ rownames(normpopt)[vi] ,,], poptDistr[vi] )
}
res
})
# link on ?twDEMC, better delete resulting twDEMC.Rd
R.methodsS3::setMethodS3("twDEMC","default", function(
### Initialize and call \code{\link{twDEMCBlockInt}}.
Zinit, ##<< initial population: a numeric array (d x M0 x Npop) or an object of class twDEMC
... ##<< further arguments to \code{\link{twDEMCBlockInt}}
){
stop("unknown type of initialization argument Zinit")
##seealso<<
## \code{\link{twDEMCBlockInt}}
## \code{\link{twDEMC.array}}
## \code{\link{twDEMC.twDEMC}}
})
.makeTableDensityTest <- function(){
df <- c(1:10,15,20,30,50)
tmp <- -sapply( df, function(dfi){ signif(getRLogDenQuantile(NULL, maxLogDen=0, df=dfi),3) })
plot( tmp~df)
ds <- data.frame(df=df,diffLogDen=tmp)
cat(twDf2wikiTable(ds)) # package twMisc
lm1 <- lm(tmp~df)
abline(lm1)
lm2 <- lm(tmp[-(1:5)]~df[-(1:5)])
abline(lm2, col="red") #about 0.6 per parameter
}
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.