Nothing
### PTdensity.R
### Fit a Mixture of Polya trees for density estimation
###
### Copyright: Alejandro Jara and Tim Hanson, 2006-2012.
###
### Last modification: 25-01-2010.
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### This program is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
### General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
###
### The authors's contact information:
###
### Alejandro Jara
### Department of Statistics
### Facultad de Matematicas
### Pontificia Universidad Catolica de Chile
### Casilla 306, Correo 22
### Santiago
### Chile
### Voice: +56-2-3544506 URL : http://www.mat.puc.cl/~ajara
### Fax : +56-2-3547729 Email: atjara@uc.cl
###
### Tim Hanson
### Department of Statistics
### University of South Carolina
### 216 LeConte College,
### 1523 Greene St
### Columbia, SC 29208
### Voice: 803-777-3853 URL : http://www.stat.sc.edu/~hansont/
### Fax : 803-777-4048 Email: hansont@stat.sc.edu
###
PTdensity <- function(y,ngrid=1000,grid=NULL,prior,mcmc,state,status,data=sys.frame(sys.parent()),na.action=na.fail)
UseMethod("PTdensity")
PTdensity.default <- function(y,ngrid=1000,grid=NULL,prior,mcmc,state,status,data,na.action=na.fail)
{
#########################################################################################
# call parameters
#########################################################################################
cl <- match.call()
y <- na.action(as.matrix(y))
#########################################################################################
# data structure
#########################################################################################
nrec <- nrow(y)
nvar <- ncol(y)
if(nvar>1)
{
left <- rep(0,2)
right <- rep(0,2)
left[1] <- min(y[,1])-0.5*sqrt(var(y[,1]))
right[1] <- max(y[,1])+0.5*sqrt(var(y[,1]))
left[2] <- min(y[,2])-0.5*sqrt(var(y[,2]))
right[2] <- max(y[,2])+0.5*sqrt(var(y[,2]))
}
else
{
left <- min(y)-0.5*sqrt(var(y))
right <- max(y)+0.5*sqrt(var(y))
}
#########################################################################################
# prior information
#########################################################################################
jfr <- c(0,0)
if(nvar==1)
{
tau <- c(-1,1)
m0 <- 0
S0 <- 1
}
else
{
nu0 <- -1
tinv <- diag(1,nvar)
m0 <- rep(0,nvar)
S0 <- diag(1,nvar)
}
if(is.null(prior$mu))
{
murand <- 1
if(is.null(prior$m0))
{
jfr[1] <- 1
}
else
{
m0 <- prior$m0
S0 <- prior$S0
}
}
else
{
murand <- 0
mu <- prior$mu
}
if(is.null(prior$sigma))
{
sigmarand <- 1
if(nvar==1)
{
if(is.null(prior$tau1))
{
jfr[2] <- 1
}
else
{
tau <- c(prior$tau1,prior$tau2)
}
}
else
{
if(is.null(prior$nu0))
{
jfr[2] <- 1
}
else
{
nu0 <- prior$nu0
tinv <- prior$tinv
}
}
}
else
{
sigmarand <- 0
sigma <- prior$sigma
}
if(is.null(prior$a0))
{
ca <- -1
cb <- -1
cpar <- prior$alpha
crand <- 0
}
else
{
ca <- prior$a0
cb <- prior$b0
cpar <- 10
crand <- 1
}
ab <- c(ca,cb)
#########################################################################################
# mcmc specification
#########################################################################################
mcmcvec <- c(mcmc$nburn,mcmc$nskip,mcmc$ndisplay)
nsave <- mcmc$nsave
if(is.null(mcmc$tune1))
{
tune1 <- 1.1
}
else
{
tune1 <- mcmc$tune1
}
if(is.null(mcmc$tune2))
{
tune2 <- 1.1
}
else
{
tune2 <- mcmc$tune2
}
if(is.null(mcmc$tune3))
{
tune3 <- 1.1
}
else
{
tune3 <- mcmc$tune3
}
#########################################################################################
# output
#########################################################################################
acrate <- rep(0,3)
if(nvar==1)
{
f <- rep(0,ngrid)
}
else
{
ngrid <- as.integer(sqrt(ngrid))
f <- matrix(0,nrow=ngrid,ncol=ngrid)
fun1 <- rep(0,ngrid)
fun2 <- rep(0,ngrid)
}
thetasave <- matrix(0,nrow=nsave,ncol=nvar+nvar*(nvar+1)/2+1)
randsave <- matrix(0,nrow=nsave,ncol=nvar)
#########################################################################################
# parameters depending on status
#########################################################################################
if(status==TRUE)
{
if(nvar==1)
{
if(murand==1)
{
mu <- mean(y)
}
if(sigmarand==1)
{
sigma <- sqrt(var(y))
}
}
else
{
if(murand==1)
{
mu <- apply(y,2,mean)
}
if(sigmarand==1)
{
sigma <- matrix(var(y),nvar,nvar)
}
}
}
if(status==FALSE)
{
cpar <- state$alpha
mu <- state$mu
sigma <- state$sigma
}
#########################################################################################
# working space
#########################################################################################
acrate <- rep(0,3)
cpo <- rep(0,nrec)
if(nvar==1)
{
if(is.null(grid))
{
grid <- seq(left,right,length=ngrid)
}
else
{
grid <- as.matrix(grid)
ngrid <- nrow(grid)
grid <- as.vector(grid)
}
}
else
{
if(is.null(grid))
{
grid1<-seq(left[1],right[1],length=ngrid)
grid2<-seq(left[2],right[2],length=ngrid)
}
else
{
grid <- as.matrix(grid)
ngrid <- nrow(grid)
grid1 <- grid[,1]
grid2 <- grid[,2]
}
}
seed <- c(sample(1:29000,1),sample(1:29000,1))
#########################################################################################
# calling the fortran code
#########################################################################################
if(nvar==1)
{
if(is.null(prior$M))
{
whicho <- rep(0,nrec)
whichn <- rep(0,nrec)
foo <- .Fortran("ptdensityu",
ngrid =as.integer(ngrid),
nrec =as.integer(nrec),
y =as.double(y),
ab =as.double(ab),
murand =as.integer(murand),
sigmarand =as.integer(sigmarand),
jfr =as.integer(jfr),
m0 =as.double(m0),
s0 =as.double(S0),
tau =as.double(tau),
mcmcvec =as.integer(mcmcvec),
nsave =as.integer(nsave),
tune1 =as.double(tune1),
tune2 =as.double(tune2),
tune3 =as.double(tune3),
acrate =as.double(acrate),
f =as.double(f),
thetasave =as.double(thetasave),
cpo =as.double(cpo),
cpar =as.double(cpar),
mu =as.double(mu),
sigma =as.double(sigma),
grid =as.double(grid),
seed =as.integer(seed),
whicho =as.integer(whicho),
whichn =as.integer(whichn),
PACKAGE ="DPpackage")
}
else
{
nlevel <- prior$M
ninter <- 2**nlevel
assign <- matrix(0,nrow=nrec,ncol=nlevel)
accums <- matrix(0,nrow=nlevel,ncol=ninter)
counter <- matrix(0,nrow=nlevel,ncol=ninter)
endp <- rep(0,ninter-1)
intpn <- rep(0,nrec)
intpo <- rep(0,nrec)
prob <- rep(0,ninter)
rvecs <- matrix(0,nrow=nlevel,ncol=ninter)
foo <- .Fortran("ptdensityup",
ngrid =as.integer(ngrid),
nrec =as.integer(nrec),
y =as.double(y),
ab =as.double(ab),
murand =as.integer(murand),
sigmarand =as.integer(sigmarand),
jfr =as.integer(jfr),
m0 =as.double(m0),
s0 =as.double(S0),
tau =as.double(tau),
nlevel =as.integer(nlevel),
ninter =as.integer(ninter),
mcmcvec =as.integer(mcmcvec),
nsave =as.integer(nsave),
tune1 =as.double(tune1),
tune2 =as.double(tune2),
tune3 =as.double(tune3),
acrate =as.double(acrate),
f =as.double(f),
thetasave =as.double(thetasave),
cpo =as.double(cpo),
cpar =as.double(cpar),
mu =as.double(mu),
sigma =as.double(sigma),
grid =as.double(grid),
intpn =as.integer(intpn),
intpo =as.integer(intpo),
accums =as.double(accums),
assign =as.integer(assign),
counter =as.integer(counter),
endp =as.double(endp),
prob =as.double(prob),
rvecs =as.double(rvecs),
seed =as.integer(seed),
PACKAGE ="DPpackage")
}
}
else
{
iflag <- rep(0,nvar)
limw <- rep(0,nvar)
linf <- rep(0,nvar)
lsup <- rep(0,nvar)
muc <- rep(0,nvar)
narea <- 2**nvar
mass <- rep(0,narea)
massi <- rep(0,narea)
parti <- rep(0,nvar)
pattern <- rep(0,nvar)
patterns <- rep(0,nvar)
propv <- matrix(0,nrow=nvar,ncol=nvar)
propv1 <- matrix(0,nrow=nvar,ncol=nvar)
propv2 <- matrix(0,nrow=nvar,ncol=nvar)
s <- matrix(0,nrow=nvar,ncol=nvar)
sigmac <- matrix(0,nrow=nvar,ncol=nvar)
sigmainv <- matrix(0,nrow=nvar,ncol=nvar)
sigmainvc <- matrix(0,nrow=nvar,ncol=nvar)
vv <- rep(0,nvar)
whicho <- rep(0,nrec)
whichn <- rep(0,nrec)
workh1 <- rep(0,nvar*(nvar+1)/2)
workh2 <- rep(0,nvar*(nvar+1)/2)
workmh <- rep(0,nvar*(nvar+1)/2)
workm1 <- matrix(0,nrow=nvar,ncol=nvar)
workm2 <- matrix(0,nrow=nvar,ncol=nvar)
ybar <- rep(0,nvar)
z <- matrix(0,nrow=nrec,ncol=nvar)
zc <- matrix(0,nrow=nrec,ncol=nvar)
zwork <- rep(0,nvar)
s0 <- solve(S0)
if(is.null(prior$M))
{
foo <- .Fortran("ptmdensity",
ngrid =as.integer(ngrid),
nrec =as.integer(nrec),
nvar =as.integer(nvar),
y =as.double(y),
ab =as.double(ab),
murand =as.integer(murand),
sigmarand =as.integer(sigmarand),
jfr =as.integer(jfr),
m0 =as.double(m0),
s0 =as.double(s0),
nu0 =as.integer(nu0),
tinv =as.double(tinv),
mcmcvec =as.integer(mcmcvec),
nsave =as.integer(nsave),
tune1 =as.double(tune1),
tune2 =as.double(tune2),
tune3 =as.double(tune3),
acrate =as.double(acrate),
cpo =as.double(cpo),
f =as.double(f),
randsave =as.double(randsave),
thetasave =as.double(thetasave),
cpar =as.double(cpar),
mu =as.double(mu),
sigma =as.double(sigma),
grid1 =as.double(grid1),
grid2 =as.double(grid2),
iflag =as.integer(iflag),
whicho =as.integer(whicho),
whichn =as.integer(whichn),
limw =as.double(limw),
linf =as.double(linf),
lsup =as.double(lsup),
narea =as.integer(narea),
mass =as.double(mass),
massi =as.integer(massi),
parti =as.integer(parti),
pattern =as.integer(pattern),
patterns =as.integer(patterns),
s =as.double(s),
sigmainv =as.double(sigmainv),
sigmainvc =as.double(sigmainvc),
ybar =as.double(ybar),
z =as.double(z),
zc =as.double(zc),
zwork =as.double(zwork),
vv =as.double(vv),
workmh =as.double(workmh),
workh1 =as.double(workh1),
workh2 =as.double(workh2),
workm1 =as.double(workm1),
workm2 =as.double(workm2),
muc =as.double(muc),
sigmac =as.double(sigmac),
propv =as.double(propv),
propv1 =as.double(propv1),
propv2 =as.double(propv2),
seed =as.integer(seed),
PACKAGE ="DPpackage")
}
else
{
foo <- .Fortran("ptmdensityp",
ngrid =as.integer(ngrid),
nrec =as.integer(nrec),
nvar =as.integer(nvar),
y =as.double(y),
ab =as.double(ab),
nlevel =as.integer(prior$M),
murand =as.integer(murand),
sigmarand =as.integer(sigmarand),
jfr =as.integer(jfr),
m0 =as.double(m0),
s0 =as.double(s0),
nu0 =as.integer(nu0),
tinv =as.double(tinv),
mcmcvec =as.integer(mcmcvec),
nsave =as.integer(nsave),
tune1 =as.double(tune1),
tune2 =as.double(tune2),
tune3 =as.double(tune3),
acrate =as.double(acrate),
cpo =as.double(cpo),
f =as.double(f),
randsave =as.double(randsave),
thetasave =as.double(thetasave),
cpar =as.double(cpar),
mu =as.double(mu),
sigma =as.double(sigma),
grid1 =as.double(grid1),
grid2 =as.double(grid2),
iflag =as.integer(iflag),
whicho =as.integer(whicho),
whichn =as.integer(whichn),
limw =as.double(limw),
linf =as.double(linf),
lsup =as.double(lsup),
narea =as.integer(narea),
mass =as.double(mass),
massi =as.integer(massi),
parti =as.integer(parti),
pattern =as.integer(pattern),
patterns =as.integer(patterns),
s =as.double(s),
sigmainv =as.double(sigmainv),
sigmainvc =as.double(sigmainvc),
ybar =as.double(ybar),
z =as.double(z),
zc =as.double(zc),
zwork =as.double(zwork),
vv =as.double(vv),
workmh =as.double(workmh),
workh1 =as.double(workh1),
workh2 =as.double(workh2),
workm1 =as.double(workm1),
workm2 =as.double(workm2),
muc =as.double(muc),
sigmac =as.double(sigmac),
propv =as.double(propv),
propv1 =as.double(propv1),
propv2 =as.double(propv2),
seed =as.integer(seed),
PACKAGE ="DPpackage")
}
}
#########################################################################################
# save state
#########################################################################################
model.name <- "Bayesian Density Estimation Using MPT"
varnames<-colnames(y)
if(is.null(varnames))
{
varnames<-all.vars(cl)[1:nvar]
}
state <- list(alpha=foo$cpar,
mu=foo$mu,
sigma=matrix(foo$sigma,nrow=nvar,ncol=nvar)
)
thetasave <- matrix(foo$thetasave,nrow=nsave,ncol=(nvar+nvar*(nvar+1)/2+1))
if(nvar>1)
{
randsave <- matrix(foo$randsave,nrow=nsave,ncol=nvar)
colnames(randsave) <- varnames
}
coeff<-apply(thetasave,2,mean)
pnames1<-NULL
for(i in 1:nvar)
{
pnames1<-c(pnames1,paste("mu",varnames[i],sep=":"))
}
pnames2<-NULL
for(i in 1:nvar)
{
for(j in i:nvar)
{
if(i==j)
{
tmp <- varnames[i]
}
else
{
tmp <- paste(varnames[i],varnames[j],sep="-")
}
pnames2 <- c(pnames2,paste("sigma",tmp,sep=":"))
}
}
names(coeff) <- c(pnames1,pnames2,"alpha")
colnames(thetasave) <- c(pnames1,pnames2,"alpha")
save.state <- list(thetasave=thetasave,randsave=randsave)
if(crand==0)
{
acrate<-foo$acrate[1:2]
}
else
{
acrate<-foo$acrate
}
x1<-NULL
x2<-NULL
dens<-NULL
if(nvar==1)
{
x1<-foo$grid
dens<-foo$f
f<-foo$f
grid1<-foo$grid
grid2<-NULL
fun1<-foo$f
fun2<-NULL
}
else
{
x1 <- grid1
x2 <- grid2
dens <- matrix(foo$f,nrow=ngrid,ncol=ngrid)
f <- matrix(foo$f,nrow=ngrid,ncol=ngrid)
dist1 <- grid2[2]-grid2[1]
dist2 <- grid1[2]-grid1[1]
fun1 <- (dist1/2)*(dens[,1]+dens[,ngrid]+2*apply(dens[,2:(ngrid-1)],1,sum))
fun2 <- (dist2/2)*(dens[1,]+dens[ngrid,]+2*apply(dens[2:(ngrid-1),],2,sum))
}
z<-list(call=cl,y=y,varnames=varnames,modelname=model.name,cpo=foo$cpo,
prior=prior,mcmc=mcmc,state=state,save.state=save.state,nrec=foo$nrec,
nvar=nvar,crand=crand,coefficients=coeff,f=f,grid1=grid1,grid2=grid2,
fun1=fun1,fun2=fun2,x1=x1,x2=x2,dens=dens,acrate=acrate)
cat("\n\n")
class(z)<-"PTdensity"
return(z)
}
###
### Tools
###
### Copyright: Alejandro Jara Vallejos, 2006
### Last modification: 28-11-2006.
###
"print.PTdensity"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\n")
cat("Posterior Predictive Distributions (log):\n")
print.default(format(summary(log(x$cpo)), digits = digits), print.gap = 2,
quote = FALSE)
cat("\nPosterior Inference of Parameters:\n")
print.default(format(x$coefficients, digits = digits), print.gap = 2,
quote = FALSE)
cat("\nAcceptance Rate for Metropolis Step = ",x$acrate,"\n")
cat("\nNumber of Observations:",x$nrec)
cat("\nNumber of Variables:",x$nvar,"\n")
cat("\n\n")
invisible(x)
}
"summary.PTdensity"<-function(object, hpd=TRUE, ...)
{
stde<-function(x)
{
n<-length(x)
return(sd(x)/sqrt(n))
}
hpdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[1],a$aupp[1]))
}
pdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[2],a$aupp[2]))
}
thetasave<-object$save.state$thetasave
nvar<-object$nvar
if(object$crand==0)
{
dimen<-(nvar+nvar*(nvar+1)/2)
mat<-matrix(thetasave[,1:dimen],ncol=2)
}
else
{
dimen<-(nvar+nvar*(nvar+1)/2+1)
mat<-thetasave[,1:dimen]
}
coef.p<-object$coefficients[1:dimen]
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
names(coef.m)<-names(object$coefficients[1:dimen])
names(coef.sd)<-names(object$coefficients[1:dimen])
names(coef.se)<-names(object$coefficients[1:dimen])
names(coef.l)<-names(object$coefficients[1:dimen])
names(coef.u)<-names(object$coefficients[1:dimen])
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans <- c(object[c("call", "modelname")])
ans$coefficients<-coef.table
### CPO
ans$cpo<-object$cpo
ans$acrate<-object$acrate
ans$nrec<-object$nrec
ans$nvar<-object$nvar
class(ans) <- "summaryPTdensity"
return(ans)
}
"print.summaryPTdensity"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\n")
cat("Posterior Predictive Distributions (log):\n")
print.default(format(summary(log(x$cpo)), digits = digits), print.gap = 2,
quote = FALSE)
if (length(x$coefficients)) {
cat("\nBaseline parameters:\n")
print.default(format(x$coefficients, digits = digits), print.gap = 2,
quote = FALSE)
}
else cat("No coefficients\n")
cat("\nAcceptance Rate for Metropolis Step = ",x$acrate,"\n")
cat("\nNumber of Observations:",x$nrec)
cat("\nNumber of Variables:",x$nvar,"\n")
cat("\n\n")
invisible(x)
}
"plot.PTdensity"<-function(x, ask=TRUE, output="density", param=NULL, hpd=TRUE, nfigr=1, nfigc=1, col="#bdfcc9", ...)
{
fancydensplot1<-function(x, hpd=TRUE, npts=200, xlab="", ylab="", main="",col="#bdfcc9", ...)
# Author: AJV, 2006
#
{
dens <- density(x,n=npts)
densx <- dens$x
densy <- dens$y
meanvar <- mean(x)
densx1 <- max(densx[densx<=meanvar])
densx2 <- min(densx[densx>=meanvar])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ymean <- densy1 + ((densy2-densy1)/(densx2-densx1))*(meanvar-densx1)
if(hpd==TRUE)
{
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
n<-length(x)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(x),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
xlinf<-a$alow[1]
xlsup<-a$aupp[1]
}
else
{
xlinf <- quantile(x,0.025)
xlsup <- quantile(x,0.975)
}
densx1 <- max(densx[densx<=xlinf])
densx2 <- min(densx[densx>=xlinf])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylinf <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlinf-densx1)
densx1 <- max(densx[densx<=xlsup])
densx2 <- min(densx[densx>=xlsup])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylsup <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlsup-densx1)
plot(0.,0.,xlim = c(min(densx), max(densx)), ylim = c(min(densy), max(densy)),
axes = F,type = "n" , xlab=xlab, ylab=ylab, main=main, cex=1.2)
xpol<-c(xlinf,xlinf,densx[densx>=xlinf & densx <=xlsup],xlsup,xlsup)
ypol<-c(0,ylinf,densy[densx>=xlinf & densx <=xlsup] ,ylsup,0)
polygon(xpol, ypol, border = FALSE,col=col)
lines(c(min(densx), max(densx)),c(0,0),lwd=1.2)
segments(min(densx),0, min(densx),max(densy),lwd=1.2)
lines(densx,densy,lwd=1.2)
segments(meanvar, 0, meanvar, ymean,lwd=1.2)
segments(xlinf, 0, xlinf, ylinf,lwd=1.2)
segments(xlsup, 0, xlsup, ylsup,lwd=1.2)
axis(1., at = round(c(xlinf, meanvar,xlsup), 2.), labels = T,pos = 0.)
axis(1., at = round(seq(min(densx),max(densx),length=15), 2.), labels = F,pos = 0.)
axis(2., at = round(seq(0,max(densy),length=5), 2.), labels = T,pos =min(densx))
}
"bivk"<-function(x, y, h, n = 25, lims = c(range(x), range(y)))
{
nx <- length(x)
if (length(y) != nx)
stop("Data vectors must be the same length")
gx <- seq(lims[1], lims[2], length = n)
gy <- seq(lims[3], lims[4], length = n)
if (missing(h))
h <- c(band(x), band(y))
h <- h/4
ax <- outer(gx, x, "-")/h[1]
ay <- outer(gy, y, "-")/h[2]
z <- matrix(dnorm(ax), n, nx) %*% t(matrix(dnorm(ay), n,
nx))/(nx * h[1] * h[2])
return(list(x = gx, y = gy, z = z))
}
"band"<-function(x)
{
r <- quantile(x, c(0.25, 0.75))
h <- (r[2] - r[1])/1.34
4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
}
if(is(x, "PTdensity"))
{
if(output=="density")
{
# Density estimation
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1),nrow=nfigr,ncol=nfigc,byrow=TRUE))
if(x$nvar==1)
{
title1<-paste("Density of",x$varnames[1],sep=' ')
aa<-hist(x$y[,1],plot=F,)
maxx<-max(aa$intensities+aa$density)+0.1*max(aa$intensities+aa$density)
miny<-min(x$y[,1])
maxy<-max(x$y[,1])
deltay<-(maxy-miny)*0.2
miny<-miny-deltay
maxy<-maxy+deltay
hist(x$y[,1],probability=T,xlim=c(min(x$grid1),max(x$grid1)),ylim=c(0,maxx),nclas=25,main=title1,xlab="values", ylab="density")
lines(x$x1,x$dens,lwd=2)
}
if(x$nvar==2)
{
title1<-paste("Density of",x$varnames[1],sep=' ')
aa<-hist(x$y[,1],plot=F,)
maxx<-max(aa$intensities+aa$density)+0.1*max(aa$intensities+aa$density)
miny<-min(x$y[,1])
maxy<-max(x$y[,1])
deltay<-(maxy-miny)*0.2
miny<-miny-deltay
maxy<-maxy+deltay
hist(x$y[,1],probability=T,xlim=c(min(x$grid1),max(x$grid1)),ylim=c(0,maxx),nclas=25,main=title1,xlab="values", ylab="density")
lines(x$grid1,x$fun1,lwd=2)
title1<-paste("Density of",x$varnames[2],sep=' ')
aa<-hist(x$y[,2],plot=F,)
maxx<-max(aa$intensities+aa$density)+0.1*max(aa$intensities+aa$density)
miny<-min(x$y[,2])
maxy<-max(x$y[,2])
deltay<-(maxy-miny)*0.2
miny<-miny-deltay
maxy<-maxy+deltay
hist(x$y[,2],probability=T,xlim=c(min(x$grid2),max(x$grid2)),ylim=c(0,maxx),nclas=25,main=title1,xlab="values", ylab="density")
lines(x$grid2,x$fun2,lwd=2)
varsn<-paste(x$varnames[1],x$varnames[2],sep="-")
title1<-paste("Predictive Density of ",varsn,sep='')
xx<-matrix(x$grid1,ncol=1)
yy<-matrix(x$grid2,ncol=1)
z<-x$f
colnames(xx)<-x$varnames[1]
colnames(yy)<-x$varnames[2]
contour(xx,yy,z,main=title1,xlab=x$varnames[1],ylab=x$varnames[2])
persp(xx,yy,z,xlab=x$varnames[1],ylab=x$varnames[2],zlab="density",theta=-30,phi=15,expand = 0.9, ltheta = 120,main=title1)
}
if(x$nvar>2)
{
for(i in 1:x$nvar)
{
title1<-paste("Density of",x$varnames[i],sep=' ')
aa<-hist(x$y[,i],plot=F)
maxx<-max(aa$intensities+aa$density)+0.1*max(aa$intensities+aa$density)
miny<-min(x$y[,i])-0.5*sqrt(var(x$y[,i]))
maxy<-max(x$y[,i])+0.5*sqrt(var(x$y[,i]))
deltay<-(maxy-miny)*0.2
miny<-miny-deltay
maxy<-maxy+deltay
hist(x$y[,i],probability=T,xlim=c(miny,maxy),ylim=c(0,maxx),nclas=25,main=title1,xlab="values", ylab="density")
lines(density(x$save.state$randsave[,i]),lwd=2)
}
for(i in 1:(x$nvar-1))
{
vectmp<-x$y[,i]
xlim<-c(min(vectmp)-0.5*sqrt(var(vectmp)),max(vectmp)+0.125*sqrt(var(vectmp)))
for(j in (i+1):x$nvar)
{
vectmp<-x$y[,j]
ylim<-c(min(vectmp)-0.5*sqrt(var(vectmp)),max(vectmp)+0.125*sqrt(var(vectmp)))
varsn<-paste(x$varnames[i],x$varnames[j],sep="-")
title1<-paste("Predictive Density of ",varsn,sep='')
xx<-x$save.state$randsave[,i]
yy<-x$save.state$randsave[,j]
est<-bivk(xx,yy,n=200)
contour(est,xlim=xlim,ylim=ylim,main=title1,xlab=x$varnames[i],ylab=x$varnames[j])
persp(est,theta=-30,phi=15,expand = 0.9, ltheta = 120,main=title1,
xlab=x$varnames[i],ylab=x$varnames[j],zlab="density")
}
}
}
}
else
{
if(is.null(param))
{
pnames<-colnames(x$save.state$thetasave)
n<-dim(x$save.state$thetasave)[2]
cnames<-names(x$coefficients)
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
for(i in 1:(n-1))
{
title1<-paste("Trace of",pnames[i],sep=" ")
title2<-paste("Density of",pnames[i],sep=" ")
plot(x$save.state$thetasave[,i],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot1(x$save.state$thetasave[,i],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
if(is.null(x$prior$a0))
{
cat("")
}
else
{
title1<-paste("Trace of",pnames[n],sep=" ")
title2<-paste("Density of",pnames[n],sep=" ")
plot(x$save.state$thetasave[,n],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot1(x$save.state$thetasave[,n],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
}
else
{
pnames<-colnames(x$save.state$thetasave)
n<-dim(x$save.state$thetasave)[2]
poss<-0
for(i in 1:n)
{
if(pnames[i]==param)poss=i
}
if (poss==0)
{
stop("This parameter is not present in the original model.\n")
}
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
title1<-paste("Trace of",pnames[poss],sep=" ")
title2<-paste("Density of",pnames[poss],sep=" ")
plot(x$save.state$thetasave[,poss],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot1(x$save.state$thetasave[,poss],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
}
}
}
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.