swh:1:snp:ffdd0a7d2c8ea15ad41d45b3b178f668bd942287
Tip revision: 4af1e2789bcea7df3c1775a53cd05b37ec3185d0 authored by Derek Young on 05 December 2022, 13:30:02 UTC
version 2.0.0
version 2.0.0
Tip revision: 4af1e27
spregmix.R
spregmix = function (lmformula, bw = NULL, constbw = FALSE,
bwmult = 0.9, # gives Silverman's rule; or use 1.06 for Scott's
z.hat = NULL, symm = TRUE, betamethod="LS",
m = ifelse(is.null(z.hat), 2, ncol(z.hat)),
epsilon = 1e-04, maxit = 1000, verbose = FALSE,
...)
{
# m, not k, is the number of components here
y <- model.frame(lmformula, ...)[,1]
x <- model.matrix(lmformula, ...)
n <- length(y)
p <- ncol(x)
tt0 <- proc.time() # for total time
# Now initialize z.hat, which is the only thing needed to start the iterations
# as long as we use a method such as least squares to find beta
if(is.null(z.hat)) {
z.hat <- matrix(runif(n*m), n, m)
z.hat <- sweep(z.hat, 1, rowSums(z.hat), "/")
}
L1norm <- function(beta, y, x, p) sum(p*abs(y - x %*% beta))
nploglik <- function(beta, y, x, p, bw, res, wts, symm) {
sum(p*log(wkde(x=res, u=as.vector(y-x%*%beta), w=wts, bw=bw, sym=symm)))
}
lambda <- matrix(0, maxit, m)
beta <- matrix(0, p, m)
for (j in 1:m) { # Initialize beta to weighted LS solution
wx <- sweep(x, 1, z.hat[,j], "*")
beta[,j] <- solve(t(wx) %*% x, t(wx) %*% y)
}
iter <- 0
finished <- FALSE
# j <- 0
#plot(tonedata)
SPweight <- ifelse(betamethod=="transition", 0, 1)
while (!finished) {
# j <- j + 1; if (j>m) {
#browser();
#plot(x[,2],y);
# j <- 1}
iter <- iter + 1
t0 <- proc.time()
## M-step
lambda[iter, ] <- colMeans(z.hat)
oldbeta <- beta
for (j in 1:m) {
#abline(beta[1,j],beta[2,j], col=1+j, lwd=3, lty=2)
if (betamethod=="LS") {
wx <- sweep(x, 1, z.hat[,j], "*")
beta[,j] <- solve(t(wx) %*% x, t(wx) %*% y) # Weighted least squares solution
} else if (betamethod=="L1") { # Weighted L1 solution
beta[,j] <- optim(par=beta[,j], fn=L1norm, y=y, x=x, p=z.hat[,j])$par
} else if (betamethod=="transition") { # transition from LS to NP
wx <- sweep(x, 1, z.hat[,j], "*")
beta[,j] <- SPweight *
optim(par=beta[,j], fn=nploglik, y=y, x=x, p=z.hat[,j],
bw=bw, res=as.vector(y-x%*%beta),
wts=as.vector(z.hat), symm=symm,
control=list(fnscale=-1))$par + # NP loglik
(1-SPweight) *
solve(t(wx) %*% x, t(wx) %*% y) # Weighted least squares
SPweight <- min(1, SPweight + (1e-4)*2^(iter-1))
} else { # Nonparametric loglikelihood
beta[,j] <- optim(par=beta[,j], fn=nploglik, y=y, x=x, p=z.hat[,j],
bw=bw, res=as.vector(y-x%*%beta),
wts=as.vector(z.hat), symm=symm,
control=list(fnscale=-1))$par
#res=as.vector(y-x%*%beta)
#plot(res,wkde(x=res, u=res, w=as.vector(z.hat), bw=bw, sym=symm))
#browser()
}
## Here, we might try other methods of estimating beta.
}
## update the bandwidth, if necssary
xbetavec <- as.double(x %*% beta)
zhatvec <- as.double(z.hat)
if (is.null(bw) || !constbw) {
res <- y-xbetavec
np.sigma <- sqrt(sum(res^2 * zhatvec)/(n-1))
bw <- bwmult / n^(1/5) * min(np.sigma, wiqr<-wIQR(wt=zhatvec, x=res)/1.34)
#print(c(np.sigma, wiqr))
}
## density estimation step
if (symm) {
ans <- .C(C_KDEsymloc2, n=as.integer(n), m=as.integer(m),
mu = xbetavec, y=as.double(y), bw=as.double(bw),
z=as.double(z.hat), f=double(n*m))
} else {
ans <- .C(C_KDEloc2, n=as.integer(n), m=as.integer(m),
mu = xbetavec, y=as.double(y), bw=as.double(bw),
z=as.double(z.hat), f=double(n*m))
}
fkernel <- matrix(ans$f, ncol=m)
lambda.f <- sweep(fkernel, 2, lambda[iter,], "*")
## E-step (for next iteration)
z.hat <- lambda.f/rowSums(lambda.f)
## Test convergence criteria
finished <- (iter > 1 && SPweight == 1 &&
(iter >= maxit ||
max(max(abs(beta-oldbeta)),
max(abs(lambda[iter,]-lambda[iter-1,]))) < epsilon
))
## Print message for each iteration
if (verbose) {
t1 <- proc.time()
cat("iteration ", iter, " lambda ", round(lambda[iter,], 4),
" bandwidth ", round(bw, 3))
cat(" time", (t1 - t0)[3], "\n")
}
}
## Print final message after convergence
if (verbose) {
tt1 <- proc.time()
cat("lambda ", round(lambda[iter,], 4))
cat(", total time", (tt1 - tt0)[3], "s\n")
}
res <- ans$y - ans$mu
xx <- res
yy <- ans$f
if (symm) {
xx <- c(xx, -res) # exploiting symmetry
yy <- c(yy, yy)
}
ox <- order(xx)
np.sigma <- sqrt(sum(res^2 * ans$z)/(n-1))
loglik <- rep(0, n)
for(j in 1:m) {
loglik <- loglik + z.hat[,j]*wkde(as.vector(y-x%*%beta[,j]),
w=as.vector(z.hat[,j]), bw=bw, sym=symm)
}
loglik <- sum(log(loglik))
a <- list(x=x, y=y, lambda = lambda[1:iter,], beta = beta,
posterior = z.hat, np.stdev = np.sigma, bandwidth=bw,
density.x = xx[ox], density.y = yy[ox], symmetric=symm,
loglik = loglik,
ft="regmixEM")
class(a) = "npEM"
a
}