Nothing
## Fit density model
ssden1 <- function(formula,type=NULL,data=list(),alpha=1.4,
weights=NULL,subset,na.action=na.omit,
id.basis=NULL,nbasis=NULL,seed=NULL,
domain=as.list(NULL),quad=NULL,
prec=1e-7,maxiter=30)
{
## Obtain model frame and model terms
mf <- match.call()
mf$type <- mf$alpha <- NULL
mf$id.basis <- mf$nbasis <- mf$seed <- NULL
mf$domain <- mf$quad <- mf$quad <- NULL
mf$prec <- mf$maxiter <- NULL
mf[[1]] <- as.name("model.frame")
mf <- eval(mf,sys.frame(sys.parent()))
cnt <- model.weights(mf)
mf$"(weights)" <- NULL
## Use ssden for 1-D estimation
if (dim(mf)[2]==1) stop("use ssden to estimate 1-D density")
## Generate sub-basis
nobs <- dim(mf)[1]
if (is.null(id.basis)) {
if (is.null(nbasis)) nbasis <- max(30,ceiling(10*nobs^(2/9)))
if (nbasis>=nobs) nbasis <- nobs
if (!is.null(seed)) set.seed(seed)
id.basis <- sample(nobs,nbasis,prob=cnt)
}
else {
if (max(id.basis)>nobs|min(id.basis)<1)
stop("gss error in ssden1: id.basis out of range")
nbasis <- length(id.basis)
}
## Set domain and type, generate rho and quadrature
if (is.null(quad)) quad <- as.list(NULL)
rho <- rho.log <- as.list(NULL)
rho.int <- rho.int2 <- NULL
for (xlab in names(mf)) {
x <- mf[[xlab]]
if (is.factor(x)) {
## factor variable
domain[[xlab]] <- NULL
wt <- as.numeric(table(x))
rho[[xlab]] <- wt/sum(wt)
quad[[xlab]] <- list(pt=unique(x),wt=rho[[xlab]])
rho.log[[xlab]] <- log(rho[[xlab]])
rho.int <- c(rho.int,sum(rho[[xlab]]*log(rho[[xlab]])))
rho.int2 <- c(rho.int2,sum(rho[[xlab]]*(log(rho[[xlab]]))^2))
}
if (is.vector(x)&!is.factor(x)) {
## numerical vector
if (is.null(domain[[xlab]])) {
mn <- min(x)
mx <- max(x)
domain[[xlab]] <- c(mn,mx)+c(-1,1)*(mx-mn)*.05
}
else domain[[xlab]] <- c(min(domain[[xlab]]),max(domain[[xlab]]))
if (is.null(type[[xlab]]))
type[[xlab]] <- list("cubic",domain[[xlab]])
else {
if (length(type[[xlab]])==1)
type[[xlab]] <- list(type[[xlab]][[1]],domain[[xlab]])
}
form <- as.formula(paste("~",xlab))
rho[[xlab]] <- ssden(form,data=mf,type=type[xlab],
domain=data.frame(domain[xlab]),
alpha=2,id.basis=id.basis)
qd.wk <- rho[[xlab]]$quad
rho.wk <- dssden(rho[[xlab]],qd.wk$pt)
qd.wk$pt <- qd.wk$pt[[1]]
qd.wk$wt <- rho.wk*qd.wk$wt
quad[[xlab]] <- qd.wk
rho.log[[xlab]] <- log(rho.wk)
rho.int <- c(rho.int,sum(log(rho.wk)*qd.wk$wt))
rho.int2 <- c(rho.int2,sum((log(rho.wk))^2*qd.wk$wt))
}
if (is.matrix(x)) {
## numerical matrix
if (is.null(quad[[xlab]])|is.null(quad))
stop("gss error in ssden1: no default quadrature")
else {
qd.wk <- quad[[xlab]]
qd.wk$pt <- data.frame(I(qd.wk$pt))
colnames(qd.wk$pt) <- xlab
form <- as.formula(paste("~",xlab))
rho[[xlab]] <- ssden(form,data=mf,type=type[xlab],quad=qd.wk,
alpha=2,id.basis=id.basis)
rho.wk <- dssden(rho[[xlab]],qd.wk$pt)
quad[[xlab]]$wt <- rho.wk*quad[[xlab]]$wt
rho.log[[xlab]] <- log(rho.wk)
rho.int <- c(rho.int,sum(log(rho.wk)*quad[[xlab]]$wt))
rho.int2 <- c(rho.int2,sum((log(rho.wk))^2*quad[[xlab]]$wt))
}
}
}
## Generate terms
term <- mkterm(mf,type)
term$labels <- term$labels[term$labels!="1"]
int <- mkint(mf,type,id.basis,quad,term,rho.log,rho.int)
## Generate s, r, and q
s <- r <- NULL
nq <- 0
for (label in term$labels) {
x <- mf[,term[[label]]$vlist]
x.basis <- mf[id.basis,term[[label]]$vlist]
nphi <- term[[label]]$nphi
nrk <- term[[label]]$nrk
if (nphi) {
phi <- term[[label]]$phi
for (i in 1:nphi) s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
}
if (nrk) {
rk <- term[[label]]$rk
for (i in 1:nrk) {
nq <- nq+1
r <- array(c(r,rk$fun(x,x.basis,nu=i,env=rk$env,out=TRUE)),c(nobs,nbasis,nq))
}
}
}
if (!is.null(s)) {
nnull <- dim(s)[2]
## Check s rank
if (qr(s)$rank<nnull)
stop("gss error in ssden1: unpenalized terms are linearly dependent")
}
## Use ssden for 1-D estimation
if (nq==1) stop("use ssden to estimate density on 1-D continuous domain")
## Fit the model
z <- mspdsty1(s,r,id.basis,cnt,int,prec,maxiter,alpha)
## Brief description of model terms
desc <- NULL
for (label in term$labels)
desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
desc <- rbind(desc,apply(desc,2,sum))
rownames(desc) <- c(term$labels,"total")
colnames(desc) <- c("Unpenalized","Penalized")
## Return the results
obj <- c(list(call=match.call(),mf=mf,cnt=cnt,terms=term,desc=desc,alpha=alpha,
domain=domain,rho=rho,rho.int=rho.int,rho.int2=rho.int2,quad=quad,
id.basis=id.basis,int=int),z)
class(obj) <- c("ssden1","ssden")
obj
}
## Fit multiple smoothing parameter density
mspdsty1 <- function(s,r,id.basis,cnt,int,prec,maxiter,alpha)
{
nq <- dim(r)[3]
## initialization
theta <- -log10(apply(r[id.basis,,],3,function(x)sum(diag(x))))
r.wk <- int.r.wk <- 0
for (i in 1:nq) {
r.wk <- r.wk + 10^theta[i]*r[,,i]
int.r.wk <- int.r.wk + 10^theta[i]*int$r[,i]
}
int.wk <- list(r=int.r.wk,s=int$s)
## theta adjustment
z <- sspdsty1(s,r.wk,r.wk[id.basis,],cnt,int.wk,prec,maxiter,alpha)
theta <- theta + z$theta
r.wk <- int.r.wk <- 0
for (i in 1:nq) theta[i] <- 2*theta[i] + log10(t(z$c)%*%r[id.basis,,i]%*%z$c)
for (i in 1:nq) {
r.wk <- r.wk + 10^theta[i]*r[,,i]
int.r.wk <- int.r.wk + 10^theta[i]*int$r[,i]
}
int.wk <- list(r=int.r.wk,s=int$s)
## lambda search
z <- sspdsty1(s,r.wk,r.wk[id.basis,],cnt,int.wk,prec,maxiter,alpha)
lambda <- z$lambda
theta <- theta + z$theta
cd <- c(z$c,z$d)
scal <- NULL
## return
z$theta <- theta
return(z)
}
## Fit single smoothing parameter density
sspdsty1 <- function(s,r,q,cnt,int,prec,maxiter,alpha)
{
nxi <- dim(r)[2]
nobs <- dim(r)[1]
if (!is.null(s)) nnull <- dim(s)[2]
else nnull <- 0
nxis <- nxi+nnull
if (is.null(cnt)) cnt <- 0
if (sum(cnt)) wt <- cnt/sum(cnt)
else wt <- 1/nobs
## cv function
cv <- function(lambda) {
fit <- .Fortran("dnewton10",
cd=as.double(cd), as.integer(nxis),
as.double(10^(lambda+theta)*q), as.integer(nxi),
as.double(cbind(10^theta*r,s)), as.integer(nobs),
as.double(sum(cnt)), as.double(cnt),
as.double(c(10^theta*int$r,int$s)),
as.double(prec), as.integer(maxiter),
as.double(.Machine$double.eps), integer(nxis),
wk=double(2*nobs+nxis*(nxis+3)),
info=integer(1),PACKAGE="gss")
if (fit$info==1) stop("gss error in ssden1: Newton iteration diverges")
if (fit$info==2) warning("gss warning in ssden1: Newton iteration fails to converge")
aa <- fit$wk[1:nobs]
assign("cd",fit$cd,inherits=TRUE)
eta0 <- cbind(10^theta*r,s)%*%cd
wwt <- wt*exp(-eta0)
wwt <- wwt/sum(wwt)
assign("scal",sum(wt*exp(-eta0)),inherits=TRUE)
trc <- sum(wwt*exp(aa/(1-aa)))-1
cv <- sum(c(10^theta*int$r,int$s)*cd) + log(scal) + alpha*trc
alpha.wk <- max(0,log.la0-lambda-5)*(3-alpha) + alpha
alpha.wk <- min(alpha.wk,3)
adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*trc,0)
cv+adj
}
## initialization
mu.r <- apply(wt*r,2,sum)
v.r <- apply(wt*r^2,2,sum)
mu.s <- apply(wt*s,2,sum)
v.s <- apply(wt*s^2,2,sum)
if (is.null(s)) theta <- 0
else theta <- log10(sum(v.s-mu.s^2)/nnull/sum(v.r-mu.r^2)*nxi) / 2
log.la0 <- log10(sum(v.r-mu.r^2)/sum(diag(q))) + theta
## lambda search
cd <- rep(0,nxi+nnull)
la <- log.la0
tol <- 0
scal <- NULL
mn0 <- log.la0-6
mx0 <- log.la0+6
repeat {
mn <- max(la-1,mn0)
mx <- min(la+1,mx0)
zz <- nlm0(cv,c(mn,mx))
if ((min(zz$est-mn,mx-zz$est)>=1e-1)||
(min(zz$est-mn0,mx0-zz$est)<1e-1)) break
else la <- zz$est
}
## return
jk1 <- cv(zz$est)
c <- cd[1:nxi]
if (nnull) d <- cd[nxi+(1:nnull)]
else d <- NULL
list(lambda=zz$est,theta=theta,c=c,d=d,scal=scal,cv=zz$min)
}
## Calculate integrals of phi and rk for ssden1
mkint <- function(mf,type,id.basis,quad,term,rho,rho.int)
{
## Obtain model terms
mt <- attr(mf,"terms")
xvars <- as.character(attr(mt,"variables"))[-1]
xfacs <- attr(mt,"factors")
term.labels <- labels(mt)
vlist <- xvars[as.logical(apply(xfacs,1,sum))]
## Set types for marginals
var.type <- NULL
for (xlab in vlist) {
x <- mf[,xlab]
if (!is.null(type[[xlab]])) {
## Check consistency and set default parameters
type.wk <- type[[xlab]][[1]]
if
(!(type.wk%in%c("ordinal","nominal","cubic","linear","per",
"cubic.per","linear.per","tp","sphere","custom")))
stop("gss error in mkint: unknown type")
if (type.wk%in%c("ordinal","nominal")) {
par.wk <- NULL
if (!is.factor(x))
stop("gss error in mkint: wrong type")
}
if (type.wk%in%c("cubic","linear")) {
if (length(type[[xlab]])==1) {
mn <- min(x)
mx <- max(x)
par.wk <- c(mn,mx)+c(-1,1)*.05*(mx-mn)
}
else par.wk <- type[[xlab]][[2]]
if (is.factor(x)|!is.vector(x))
stop("gss error in mkint: wrong type")
}
if (type.wk%in%c("per","cubic.per","linear.per")) {
if (type.wk=="per") type.wk <- "cubic.per"
if (length(type[[xlab]])==1)
stop("gss error in mkint: missing domain of periodicity")
else par.wk <- type[[xlab]][[2]]
if (is.factor(x)|!is.vector(x))
stop("gss error in mkint: wrong type")
}
if (type.wk=="tp") {
if (length(type[[xlab]])==1)
par.wk <- list(order=2,mesh=x,weight=1)
else {
par.wk <- par.wk1 <- type[[xlab]][[2]]
if (length(par.wk1)==1)
par.wk <- list(order=par.wk1,mesh=x,weight=1)
if (is.null(par.wk$mesh)) par.wk$mesh <- x
if (is.null(par.wk$weight)) par.wk$weight <- 1
}
if (dim(as.matrix(x))[2]!=dim(as.matrix(par.wk$mesh))[2])
stop("gss error in mkint: wrong dimension in normalizing mesh")
}
if (type.wk=="sphere") {
if (length(type[[xlab]])==1)
par.wk <- 2
else par.wk <- type[[xlab]][[2]]
if (!(par.wk%in%(2:4)))
stop("gss error in mkint: spherical order not implemented")
}
if (type.wk=="custom") par.wk <- type[[xlab]][[2]]
}
else {
## Set default types
if (is.factor(x)) {
## categorical variable
if (is.ordered(x)) type.wk <- "ordinal"
else type.wk <- "nominal"
par.wk <- NULL
}
else {
## numerical variable
if (is.vector(x)) {
type.wk <- "cubic"
mn <- min(x)
mx <- max(x)
par.wk <- c(mn,mx)+c(-1,1)*.05*(mx-mn)
}
else {
type.wk <- "tp"
par.wk <- list(order=2,mesh=x,weight=1)
}
}
}
var.type[[xlab]] <- list(type.wk,par.wk)
}
## Create phi and rk
nbasis <- length(id.basis)
nvar <- length(names(mf))
s <- r <- s.rho <- r.rho <- NULL
ns <- nq <- 0
for (label in term.labels) {
ns <- ns+term[[label]]$nphi
nq <- nq+term[[label]]$nrk
vlist <- xvars[as.logical(xfacs[,label])]
x <- mf[,vlist]
dm <- length(vlist)
phi <- rk <- NULL
if (dm==1) {
type.wk <- var.type[[vlist]][[1]]
xx <- mf[id.basis,vlist]
xmesh <- quad[[vlist]]$pt
if (type.wk%in%c("nominal","ordinal")) {
## factor variable
if (type.wk=="nominal") fun <- mkrk.nominal(levels(x))
else fun <- mkrk.ordinal(levels(x))
if (nlevels(x)>2) {
## rk
rk <- fun$fun(xmesh,xx,fun$env,TRUE)
}
else {
## phi
wk <- as.factor(names(fun$env$code)[1])
phi <- fun$fun(xmesh,wk,fun$env)
}
}
if (type.wk=="cubic") {
## cubic splines
range <- var.type[[vlist]][[2]]
## phi
phi.fun <- mkphi.cubic(range)
phi <- phi.fun$fun(xmesh,1,phi.fun$env)
## rk
rk.fun <- mkrk.cubic(range)
rk <- rk.fun$fun(xmesh,xx,rk.fun$env,TRUE)
}
if (type.wk%in%c("cubic.per","linear","linear.per","sphere")) {
## cubic periodic, linear, and linear periodic splines
range <- var.type[[vlist]][[2]]
## rk
if (type.wk=="cubic.per") rk.fun <- mkrk.cubic.per(range)
if (type.wk=="linear") rk.fun <- mkrk.linear(range)
if (type.wk=="linear.per") rk.fun <- mkrk.linear.per(range)
if (type.wk=="sphere") rk.fun <- mkrk.sphere(range)
rk <- rk.fun$fun(xmesh,xx,rk.fun$env,TRUE)
}
if (type.wk=="tp") {
## thin-plate splines
par <- var.type[[vlist]][[2]]
order <- par$order
mesh <- par$mesh
weight <- par$weight
if (is.vector(x)) xdim <- 1
else xdim <- dim(x)[2]
## phi
phi.fun <- mkphi.tp(xdim,order,mesh,weight)
nphi <- choose(xdim+order-1,xdim)-1
if (nphi>0) {
for (nu in 1:nphi) {
phi <- cbind(phi,phi.fun$fun(xmesh,nu,phi.fun$env))
}
}
## rk
rk.fun <- mkrk.tp(xdim,order,mesh,weight)
rk <- rk.fun$fun(xmesh,xx,rk.fun$env,TRUE)
}
if (type.wk=="custom") {
## user-defined
par <- var.type[[vlist]][[2]]
nphi <- par$nphi
if (nphi>0) {
phi.fun <- par$mkphi(par$env)
for (nu in 1:nphi) {
phi <- cbind(phi,phi.fun$fun(xmesh,nu,phi.fun$env))
}
}
rk.fun <- par$mkrk(par$env)
rk <- rk.fun$fun(xmesh,xx,rk.fun$env,TRUE)
}
wmesh <- quad[[vlist]]$wt
if (!is.null(phi)) {
s.rho.wk <- rho.int*sum(wmesh*phi)
s.rho.wk[names(mf)==vlist] <- sum(wmesh*rho[[vlist]]*phi)
s <- c(s,sum(wmesh*phi))
s.rho <- c(s.rho,sum(s.rho.wk))
}
if (!is.null(rk)) {
r.rho.wk <- outer(apply(wmesh*rk,2,sum),rho.int)
r.rho.wk[,names(mf)==vlist] <- apply(wmesh*rho[[vlist]]*rk,2,sum)
r <- cbind(r,apply(wmesh*rk,2,sum))
r.rho <- cbind(r.rho,apply(r.rho.wk,1,sum))
}
}
else {
bin.fac <- n.phi <- phi.list <- rk.list <- NULL
for (i in 1:dm) {
type.wk <- var.type[[vlist[i]]][[1]]
if (type.wk%in%c("nominal","ordinal")) {
## factor variable
if (type.wk=="nominal")
rk.wk <- mkrk.nominal(levels(x[[i]]))
else rk.wk <- mkrk.ordinal(levels(x[[i]]))
phi.wk <- rk.wk
n.phi <- c(n.phi,0)
bin.fac <- c(bin.fac,!(nlevels(x[[i]])>2))
}
if (type.wk=="cubic") {
## cubic or linear splines
range <- var.type[[vlist[i]]][[2]]
## phi
phi.wk <- mkphi.cubic(range)
n.phi <- c(n.phi,1)
## rk
rk.wk <- mkrk.cubic(range)
bin.fac <- c(bin.fac,0)
}
if (type.wk%in%c("cubic.per","linear","linear.per","sphere")) {
## cubic periodic, linear, or linear periodic splines
range <- var.type[[vlist[i]]][[2]]
n.phi <- c(n.phi,0)
phi.wk <- NULL
if (type.wk=="cubic.per") rk.wk <- mkrk.cubic.per(range)
if (type.wk=="linear") rk.wk <- mkrk.linear(range)
if (type.wk=="linear.per") rk.wk <- mkrk.linear.per(range)
if (type.wk=="sphere") rk.wk <- mkrk.sphere(range)
bin.fac <- c(bin.fac,0)
}
if (type.wk=="tp") {
## thin-plate splines
par <- var.type[[vlist[i]]][[2]]
order <- par$order
mesh <- par$mesh
weight <- par$weight
if (is.vector(x[[i]])) xdim <- 1
else xdim <- dim(x[[i]])[2]
phi.wk <- mkphi.tp(xdim,order,mesh,weight)
n.phi <- c(n.phi,choose(xdim+order-1,xdim)-1)
rk.wk <- mkrk.tp(xdim,order,mesh,weight)
bin.fac <- c(bin.fac,0)
}
if (type.wk=="custom") {
## user-defined
par <- var.type[[vlist[i]]][[2]]
n.phi <- c(n.phi,par$nphi)
if (par$nphi>0) phi.wk <- par$mkphi(par$env)
else phi.wk <- NULL
rk.wk <- par$mkrk(par$env)
bin.fac <- c(bin.fac,0)
}
phi.list <- c(phi.list,list(phi.wk))
rk.list <- c(rk.list,list(rk.wk))
}
## phi
id0 <- names(mf)%in%vlist
nphi <- term[[label]]$nphi
iphi <- term[[label]]$iphi
if (nphi>0) {
for (nu in 1:nphi) {
ind <- nu - 1
s.wk <- 1
s.rho.wk <- rho.int
s.rho.wk[id0] <- 1
for (i in 1:dm) {
phi.wk <- phi.list[[i]]
xmesh <- quad[[vlist[i]]]$pt
if (bin.fac[i]) {
wk <- as.factor(names(phi.wk$env$code)[1])
phi <- phi.wk$fun(xmesh,wk,phi.wk$env)
}
else {
code <- ind%%n.phi[i] + 1
ind <- ind%/%n.phi[i]
phi <- phi.wk$fun(xmesh,code,phi.wk$env)
}
wmesh <- quad[[vlist[i]]]$wt
s.wk <- s.wk*sum(wmesh*phi)
id1 <- names(mf)==vlist[i]
s.rho.wk[id1] <- s.rho.wk[id1]*sum(wmesh*rho[[vlist[i]]]*phi)
s.rho.wk[!id1] <- s.rho.wk[!id1]*sum(wmesh*phi)
}
s <- c(s,s.wk)
s.rho <- c(s.rho,sum(s.rho.wk))
}
}
## rk
n.rk <- ifelse(n.phi,2,1)
nrk <- prod(n.rk) - as.logical(nphi)
if (nrk>0) {
for (nu in 1:nrk) {
ind <- nu - !nphi
r.wk <- 1
r.rho.wk <- outer(rep(1,nbasis),rho.int)
r.rho.wk[,id0] <- 1
for (i in 1:dm) {
code <- ind%%n.rk[i] + 1
ind <- ind%/%n.rk[i]
xx <- mf[id.basis,vlist[[i]]]
xmesh <- quad[[vlist[i]]]$pt
if (code==n.rk[i]) {
rk.wk <- rk.list[[i]]
rk <- rk.wk$fun(xmesh,xx,rk.wk$env,TRUE)
}
else {
rk <- 0
phi.wk <- phi.list[[i]]
for (j in 1:n.phi[i]) {
phix <- phi.wk$fun(xmesh,j,phi.wk$env)
phiy <- phi.wk$fun(xx,j,phi.wk$env)
rk <- rk + outer(phix,phiy)
}
}
wmesh <- quad[[vlist[i]]]$wt
r.wk <- r.wk*apply(wmesh*rk,2,sum)
id1 <- names(mf)==vlist[i]
r.rho.wk[,id1] <- r.rho.wk[,id1]*apply(wmesh*rho[[vlist[i]]]*rk,2,sum)
r.rho.wk[,!id1] <- r.rho.wk[,!id1]*apply(wmesh*rk,2,sum)
}
r <- cbind(r,r.wk)
r.rho <- cbind(r.rho,apply(r.rho.wk,1,sum))
}
}
}
}
list(s=s,r=r,s.rho=s.rho,r.rho=r.rho,var.type=var.type)
}
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.