Nothing
logdbf=c("logd","b","f")
# dil.r is dilution ratio, previously named k
# verbose=1: only one line per iteration
# verbose=2: print after step 2, step 3a, step 3b
# verbose=4: print optim output
#reltol=1e-3; max.iter=20; verbose=T; model="sp"; grid.density=200; method.init="TLS"; init=NULL
prcsp=function(xvar, dil.x, yvar, dil.y,
model=c("sp","struct"), #sp: semiparametric, struct: structural
stop.when.dropping=FALSE, grid.density=200,
init=NULL, method.init=c("TLS","naive"), reltol=1e-3, max.iter=20,
try.additiona.support.sets=FALSE, # sp only
keep.history=FALSE, verbose=FALSE)
{
method.init=match.arg(method.init)
model=match.arg(model)
if(verbose) myprint(model, method.init)
# removing this line requires adding package namespace to Matrix and mosek function calls
if(model=="sp") {
have.Rmosek=requireNamespace("Rmosek")
if (!have.Rmosek) {
print("Rmosek does not load successfully")
return (NULL)
}
}
dil.r=dil.x/dil.y
stopifnot(length(xvar)==length(yvar))
n=length(xvar)
xvar=as.vector(xvar); yvar=as.vector(yvar) # strip off potential attributes; otherwise it may cause trouble with gnls
dat=data.frame(readout.x=unname(xvar), readout.y=unname(yvar))
res=list(dilution.ratio=dil.r, dilution.x=dil.x, dilution.y=dil.y, xvar=xvar, yvar=yvar)
class(res)=c("prc",class(res))
#########################################################
# Step 1: run prc to the end to initialize
if (is.null(init)) {
if (verbose) cat("\n========== run prc to initialize ================\n")
fit.init= prc (xvar, dil.x, yvar, dil.y, method=method.init, verbose=verbose)
if (!is.null(fit.init$error)) {
res$error="During initialization, prc fails to run."
cat(res$error, "\n", file=stderr())
return(res)
}
theta=c(coef(fit.init), sigma.sq=fit.init$sigma.sq)
theta=c(logc=unname(log(theta["c"])), logd=unname(log(theta["d"])), theta[-match(c("c","d"), names(theta))])
} else {
theta=init
if (length(theta)!=5) stop("init need to include sigma.sq")
}
nu=c(alpha=1,beta=1) # initial value for nu
if (verbose) {
cat("\n")
if (model=="struct") cat("initial nu:", nu, "\n")
cat("initial theta:", theta, "\n\n")
}
# define the grid and compute matrix A for the first time
u.0=seq(theta["logc"]*0.9, theta["logd"]*1.1, length=2+grid.density*1.2)
# loop through step 2 and 3.
# only theta has new.theta because we need to keep both update and previous values to decide convergence
# we only keep one iteration of support/p/nu
iterations=0
histories=list()
mix.liks.after.step.3=NULL
while(TRUE) {
iterations=iterations+1
if (verbose>=2) cat("=========== Iter "%.%iterations%.%" ==============\n") else cat(formatC(iterations,digits=2,), ") ", sep="")
#########################################################
# Step 2: update distribution of latent variable either nonparametrically or parametrically
u=u.0[ u.0<theta["logd"] & u.0>theta["logc"] ] # subset to the grid points in bound
A=compute.A (theta["logc"], theta["logd"], theta["b"], theta["f"], dil.r, theta["sigma.sq"], u, xvar, yvar)
#### Nonparametric estimate of distribution
if (model=="sp") {
sco1 <- list(sense = "min")
sco1$c <- rep(0,n)
sco1$A <- Matrix::Matrix (t(A), sparse = TRUE )
sco1$bc <- rbind(blc = rep(-Inf,length(u)), buc = rep(n,length(u)))
sco1$bx <- rbind(blx = rep(0,n), bux = rep(Inf,n))
opro <- matrix(list(), nrow=5, ncol=n, dimnames=list(c("type","j","f","g","h"), NULL))
for (i in 1:n) opro[,i] <- list("LOG", i, -1.0, 1.0, 0.0)
sco1$scopt <- list(opro=opro)
fit.mosek <- Rmosek::mosek(sco1, opts=list(verbose=1))
if (fit.mosek$sol$itr$solsta!="OPTIMAL") warning("fit.mosek$sol$itr$solsta!=OPTIMAL")
if (fit.mosek$sol$itr$prosta!="PRIMAL_AND_DUAL_FEASIBLE") warning("fit.mosek$sol$itr$prosta!=PRIMAL_AND_DUAL_FEASIBLE")
#mosek_write (sco1, "mosekfiles/yf_ex1_run_1.task", list(scofile="yf_ex1.sco"))
mix.lik.dual = sum(-log(fit.mosek$sol$itr$xx))
if (verbose>=2) {
myprint(mix.lik.dual)
plot(fit.mosek$sol$itr$xc, sco1$A %*% fit.mosek$sol$itr$xx); abline(0,1) # sanity check
title(main="Iteration "%.%iterations, outer=T, line=-1)
empty.plot()
#print(fit.mosek)
#print(table(fit.mosek$sol$itr$skc))
# print(table(fit.mosek$sol$itr$skx))
# print(summary(fit.mosek$sol$itr$xx))
# print(summary(fit.mosek$sol$itr$xc))
# print(cbind(fit.mosek$sol$itr$skx, fit.mosek$sol$itr$xx, fit.mosek$sol$itr$slx, fit.mosek$sol$itr$sux))
if(verbose>=4) print(cbind(fit.mosek$sol$itr$skc, fit.mosek$sol$itr$xc, fit.mosek$sol$itr$slc, fit.mosek$sol$itr$suc)[fit.mosek$sol$itr$skc=="UL",])
}
#print(sort(fit.mosek$sol$itr$xc, decreasing=TRUE))
# recover primal soln
support.set.1 = suppressWarnings(which(fit.mosek$sol$itr$skc=="UL")) # UL seems to choose less than ideal set
if (length(support.set.1)<=1) {
support.set.1 =order(fit.mosek$sol$itr$suc,decreasing=TRUE) [1:floor(n*.1)] # 10% sample size
support.set.2 =order(fit.mosek$sol$itr$suc,decreasing=TRUE) [1:floor(n*.2)] # 20% sample size
support.set.3 =order(fit.mosek$sol$itr$suc,decreasing=TRUE) [1:floor(n*.4)] # 40% sample size
} else {
support.set.2 =order(fit.mosek$sol$itr$suc,decreasing=TRUE) [1:min(length(support.set.1)*1.5,n*.8)] # twice as many as support.set.1, but not over 80% sample size
support.set.3 =order(fit.mosek$sol$itr$suc,decreasing=TRUE) [1:min(length(support.set.1)*2.0,n*.8)] # four times as many as support.set.1, but not over 80% sample size
}
# the following line may throw an error when mosek did not converge
support.set.5 =try(order(fit.mosek$sol$itr$suc,decreasing=TRUE) [1:which(cumsum(sort(fit.mosek$sol$itr$suc,decreasing=TRUE))>=1-1e-6)[1]]) # cumsum of suc just above 1
if (verbose>=2) cat("support sizes UL|1|2|3: ", suppressWarnings(sum(fit.mosek$sol$itr$skc=="UL")), length(support.set.1), length(support.set.2), length(support.set.3), "\n")
# Others tried:
# fit.mosek$sol$itr$xc==n sometimes return empty set
# a modification of set.1 that includes all with suc greater than any UL suc, does not make a big difference
# support.set.2 = which(fit.mosek$sol$itr$suc>0.01)
# support.set.3 = which(fit.mosek$sol$itr$suc>0.001)
# try three choices of support set and use the one with the best mix.lik.primal
support.set.list=list()
if (length(support.set.2)<length(support.set.1)) num.sets=1 else if (length(support.set.3)<length(support.set.2)) num.sets=2 else num.sets=3
if (!try.additiona.support.sets) num.sets=1
mix.lik.primals=numeric(num.sets)
for (support.idx in 1:num.sets) {
if (verbose>=2) myprint(support.idx)
support.set = get("support.set."%.%support.idx)
while (TRUE) {
p.fit=try(lm.fit(x = A[,support.set,drop=FALSE], y = 1/fit.mosek$sol$itr$xx))
if (inherits(p.fit, "try-error")) {mix.lik.primal=-Inf; break}
p.new=coef(p.fit)
if (any(is.na(p.new))) {
#stop("some p.new are NA")
p.new[is.na(p.new)]=-1
}
# p.new=fit.mosek$sol$itr$suc[support.set] # performance not good
# support.1=unique(r.x); p.1=rep(1/5,5) # mix.lik not as good as support
# support.1=support[-c(2,4)]; p.1=rep(1/5,5) # mix.lik not as good as support
# sum(log(compute.A (theta["logc"],theta["logd"],theta["b"],theta["f"], dil.r, theta["sigma.sq"], support.1, xvar, yvar) %*%p.1))
if (verbose>=2) cat("sum(p.new)=", round(sum(p.new),2), ", ", sep="")
p.new=p.new/sum(p.new)
mix.lik.primal=.Call("compute_mixlik", exp(theta["logc"]),exp(theta["logd"]),theta["b"],theta["f"], dil.r, theta["sigma.sq"], u[support.set], xvar, yvar, p.new)
if (verbose>=2) {
myprint(mix.lik.primal, digits=6)
plot(p.fit$fitted.values, 1/fit.mosek$sol$itr$xx) # making sure we are doing a good job at recovering primal soln
plot(0,0, type="n", ylim=range(c(0,p.new)), xlim=range(u[support.set]), xlab="r", ylab="p")#, main="Iteration "%.%iterations)
points(u[support.set], p.new, pch=19)
for (i in 1:length(p.new)) lines(rep(u[support.set][i],2), c(0,p.new[i]))
if (verbose>=3) {
print(cbind(fit.mosek$sol$itr$skc, fit.mosek$sol$itr$xc, fit.mosek$sol$itr$slc, fit.mosek$sol$itr$suc)[order(fit.mosek$sol$itr$suc,decreasing=TRUE)[1:n],])
myprint(sum(fit.mosek$sol$itr$suc[support.set]))
#myprint(sum(fit.mosek$sol$itr$suc[support.set.5]))
print(rbind(u[support.set], p.new))
}
}
support.set=support.set[p.new>0]
if (all(p.new>0)) break
} # end while (TRUE)
mix.lik.primals[support.idx]=mix.lik.primal
support.set.list[[support.idx]]=support.set # save this
if(verbose>=2) {old.par=par(no.readonly = TRUE); par(old.par)} # to start a new page with the existing graphics settings, but the order of mfrow vs mfcol is not saved
} # end support.idx loop
mix.lik.primal=max(mix.lik.primals)
if (verbose>=2) myprint(mix.lik.primal, digits=6)
if (mix.lik.primal<mix.lik.dual/2) {
res$error="unable to extract primal soln"
cat("\n", res$error, "\n", file=stderr())
break;
} else {
# updated support.set has been saved, p needs to be recomputed from that
support.set = support.set.list[[which.max(mix.lik.primals)]]
p.fit=lm.fit(x = A[,support.set,drop=FALSE], y = 1/fit.mosek$sol$itr$xx)
p=coef(p.fit)
p=p/sum(p)
support=u[support.set]
}
#### Structural relationship estimate
} else if (model=="struct"){
support=u
support.transf = (support-theta["logc"])/(theta["logd"]-theta["logc"])
# assume gamma distribution
optim.out = try(optim(
nu,
function(nu.fn,...) {
p=suppressWarnings(dbeta(support.transf, nu.fn["alpha"], nu.fn["beta"]))
p=p/sum(p)
-sum(log(A%*%p))
},
gr=NULL,
method="BFGS", control = list(trace=0), hessian = F
), silent=TRUE)
if(inherits(optim.out, "try-error")) {
# try again with Neld-Mead
optim.out = try(optim(
nu,
function(nu.fn,...) {
p=suppressWarnings(dbeta(support.transf, nu.fn["alpha"], nu.fn["beta"]))
p=p/sum(p)
-sum(log(A%*%p))
},
gr=NULL,
method="Nelder-Mead", control = list(trace=0), hessian = F
))
if(inherits(optim.out, "try-error")) {
res$error="fails to find nu"
cat("\n", res$error, "\n", file=stderr())
break
}
}
if (verbose>=4) print(optim.out)
nu = optim.out$par
p=dbeta(support.transf, nu["alpha"], nu["beta"])
p=p/sum(p)
if (verbose>=2) {
cat("opt nu : ", nu, ";", -optim.out$value, "\n")
}
}
# To summarize, step 2 produces two variables: support and p, and nu if structrual.
# For sp model, support is the set of grid points that have non-zero probablity estimate
# For struct model, support is the grid
new.theta=theta
max.inner.iter=25
inner.iter=0
debug.1=function(){ # trying to figure out why c_hat from struct has a problem compared to tls estimator
support.fn=u.0[u.0>new.theta["logc"] & u.0<new.theta["logd"]]
p.fn=dbeta((support.fn-new.theta["logc"])/(new.theta["logd"]-new.theta["logc"]), nu["alpha"], nu["beta"]); p.fn=p.fn/sum(p.fn)
cat(" support: ", length(support.fn), "from", support.fn[1], "to", last(support.fn), ", p[1:2 and last]:", p.fn[1:2], last(p.fn),"\n")
saved=(log(compute.A (new.theta["logc"], new.theta["logd"], new.theta["b"], new.theta["f"], dil.r, new.theta["sigma.sq"], support.fn, xvar, yvar) %*% p.fn ))
cat(saved, "\n"); saved
}
#saved.old=debug.1()
while (TRUE) {
inner.iter=inner.iter+1
new.theta.0=new.theta # new.theta.0 is used to decide convergence
#########################################################
# Step 3a: update curve parameters estimate
# optim (Nelder-Mead) performs ok, but is too slow
# optim (BFGS) has trouble with numerical differentiation
# nlminb performs well, but will return log(c) that is too close to min(support), causing problems for sigmasq if continue to step 3a
nlminb.out = nlminb(
new.theta[logdbf],
function(theta.fn,...) {
if (model=="sp") {
subset.=support>new.theta["logc"] & support<theta.fn["logd"]
support.fn=support[subset.]
p.fn=p[subset.]
p.fn=p.fn/sum(p.fn)
} else if (model=="struct") {
subset.=u.0>new.theta["logc"] & u.0<theta.fn["logd"]
support.fn=u.0[subset.]
p.fn=dbeta((support.fn-new.theta["logc"])/(theta.fn["logd"]-new.theta["logc"]), nu["alpha"], nu["beta"])
p.fn=p.fn/sum(p.fn)
}
# if divide by n, results are slightly different
-.Call("compute_mixlik", exp(new.theta["logc"]), exp(theta.fn["logd"]), theta.fn["b"], theta.fn["f"], dil.r, new.theta["sigma.sq"], support.fn, xvar, yvar, p.fn)
},
gradient=NULL, hessian = F,
control = list(trace=0)
)
new.theta[logdbf] = nlminb.out$par
if (verbose>=2) cat("opt theta : ", new.theta, ";", -nlminb.out$objective, "\n")
if (verbose>=4) print(nlminb.out)
#saved.2=debug.1(); plot(saved.2-saved.old); saved.old=saved.2
#########################################################
# Step 3b: update sigma.sq estimate
if (model=="sp") {
subset.=support>new.theta["logc"] & support<new.theta["logd"]
support.3=support[subset.]
p.3=p[subset.]
p.3=p.3/sum(p.3)
} else if (model=="struct") {
subset.=u.0>new.theta["logc"] & u.0<new.theta["logd"]
support.3 = u.0[subset.]
p.3=dbeta((support.3-new.theta["logc"])/(new.theta["logd"]-new.theta["logc"]), nu["alpha"], nu["beta"])
p.3=p.3/sum(p.3)
}
optimize.out = optimize(
f = function(sigma.sq.f,...) {
-.Call("compute_mixlik", exp(new.theta["logc"]), exp(new.theta["logd"]), new.theta["b"], new.theta["f"], dil.r, sigma.sq.f, support.3, xvar, yvar, p.3)
},
interval=c(new.theta["sigma.sq"]/5, new.theta["sigma.sq"]*5), # when the factor is 10, likelihood often drops
)
new.theta["sigma.sq"] = optimize.out$minimum
if (verbose>=2) cat("opt sigmasq: ", new.theta, ";", -optimize.out$objective, "\n")
if (verbose>=4) print(optimize.out)
# stopping rule for step 3
if (max(abs(1 - new.theta/new.theta.0)) < reltol) {
break;
} else if (inner.iter==max.inner.iter) {
break
}
} # end step 3 loop
mix.lik.after.step.3 = -optimize.out$objective
mix.liks.after.step.3=c(mix.liks.after.step.3, mix.lik.after.step.3)
if (model=="sp") histories[[iterations]]=list(support=support, p=p, theta=new.theta, mix.lik.after.step.3=mix.lik.after.step.3)
if (model=="struct") histories[[iterations]]=list(support=support, p=p, nu=nu, theta=new.theta, mix.lik.after.step.3=mix.lik.after.step.3)
if (verbose>=2) {
if (model=="sp") cat("support size:", length(support), "\n")
cat("\n")
} else if (verbose) {
cat("ll", formatC(mix.lik.after.step.3,format="g",width=7))
if (model=="sp") cat(", #support", formatC(length(support),format="g",digits=ceiling(log10(n))), sep="")
if (model=="struct") cat(", nu", nu, sep=" ")
cat(", theta", new.theta, "\n")
}
# stopping rules
if (max(abs(1 - new.theta/theta)) < reltol) {if (verbose) cat("converged\n"); break; }
theta = new.theta # update theta for the new iteration
if (iterations>=max.iter) {
cat("Stop: max iter reached\n"); break;
} else if (stop.when.dropping & mix.lik.after.step.3<mix.liks.after.step.3[ifelse(length(mix.liks.after.step.3)-1<1,1,length(mix.liks.after.step.3)-1)]) {
cat("Stop: likelihood starts decreasing.\n"); break;
} else if (mix.lik.after.step.3 < 0.5*max(mix.liks.after.step.3)) {
cat("Stop: likelihood decreases to half of the maximum.\n"); break;
}
} # end while loop
if (length(mix.liks.after.step.3)==0) {
return (res)
} else {
res$error=NULL # if there are already iterations, we need to remove res$error because this is used to determine overall fit status
}
# all parameter taken from best.iter in the history
best.iter=histories[[which.max(mix.liks.after.step.3)]]
res$support=best.iter$support
res$p=best.iter$p
res$coefficients=best.iter$theta[c("logc",logdbf)]
res$sigma.sq = best.iter$theta["sigma.sq"]
res$mixlik = best.iter$mix.lik.after.step.3
if (model=="struct") {
res$nu=best.iter$nu
}
res$A=with(res, compute.A (coefficients["logc"],coefficients["logd"],coefficients["b"],coefficients["f"], dil.r, sigma.sq, support, xvar, yvar) )
res$iterations=iterations
if (keep.history) res$histories=histories
#### compute asymptotic variance of theta_hat and xhat and yhat
# return object
res
}
# a convenience function
prcstruct=function(xvar, dil.x, yvar, dil.y,
grid.density=200, method.init=c("TLS","naive"), reltol=1e-3, max.iter=20, init=NULL,
keep.history=FALSE, verbose=FALSE)
{
method.init=match.arg(method.init)
prcsp(xvar, dil.x, yvar, dil.y, model="struct",
grid.density=grid.density, method.init=method.init, reltol=reltol, max.iter=max.iter, init=init,
keep.history=keep.history, verbose=verbose)
}
# this should match the mixture likehood obtained directly from mosek
mixlik <- function(object, ...) UseMethod("mixlik")
mixlik.prc=function(object, ...) {
#out.1=with(object, sum(log(compute.A (coefficients["logc"],coefficients["logd"],coefficients["b"],coefficients["f"], dilution.ratio, sigma.sq, support, xvar, yvar) %*%p)))
out.2=with(object, .Call("compute_mixlik", exp(coefficients["logc"]),exp(coefficients["logd"]),coefficients["b"],coefficients["f"], dilution.ratio, sigma.sq, support, xvar, yvar, p))
out.2
}
# populate A
compute.A = function(logc,logd,b,f, dil.r, sigma.sq, support, xvar, yvar) {
c=exp(logc); d=exp(logd)
.Call("compute_A", c,d,b,f, dil.r, sigma.sq, support, xvar, yvar)
# .Call is a lot faster than the R implementation below
# m.f.1=function(c,d,b,f,r,x,y,dil.r) (y - four_pl_prc(c,d,b,f, r, dil.r))^2 + (x - r)^2 # this is about four times as fast as m.f
# A=matrix(NA,n,K)
# for (i in 1:n) {
# for (k in 1:K) {
# m_ik = m.f.1(c, d, b, f, u[k], xvar[i], yvar[i], dil.r)
# A[i,k] = 1/fit$sigma.sq * exp(-m_ik/fit$sigma.sq/2)
# }
# }
}
compute.m = function(logc,logd,b,f, dil.r, sigma.sq, support, xvar, yvar) {
c=exp(logc); d=exp(logd)
n=length(xvar)
m.f.1=function(c,d,b,f,r,x,y,dil.r) (y - four_pl_prc(c,d,b,f, r, dil.r))^2 + (x - r)^2 # this is about four times as fast as m.f
K=length(support)
m=matrix(NA,n,K)
for (i in 1:n) {
for (k in 1:K) {
m[i,k] = m.f.1(c, d, b, f, support[k], xvar[i], yvar[i], dil.r)
}
}
m
}
###### some code used to be at the begining of step 3
# # initialize new.theta with an approximation to the mixture likelihood, this seems quaint a
# A=compute.A (theta["logc"], theta["logd"], theta["b"], theta["f"], dil.r, theta["sigma.sq"], support, xvar, yvar)
# best.support=support[apply(A, 1, which.max)]
# s.sqrt.k.rvar = four_pl_prc (exp(theta["logc"]), exp(theta["logd"]), theta["b"], theta["f"], best.support, dil.r^.5) # an easy to make mistake is to use k*.5 instead of k^.5
# dat.stacked=data.frame(readout=c(unname(xvar),unname(yvar)), x=rep(s.sqrt.k.rvar,2), k=rep(c(dil.r^(-1/2),dil.r^(1/2)),each=n))
# formula.gnls = as.formula( "(readout) ~ log(exp(logc)+(exp(logd)-exp(logc))/(1+k^b*(((exp(logd)-exp(logc))/(exp(x)-exp(logc)))^(1/f)-1))^f)" )
# fit.1=try(gnls(formula.gnls, data=dat.stacked, start=theta[c("c",logdbf)], # key to take only cdbf from new.theta as start for gnls
# control=gnlsControl(
# nlsTol=1e-1, # nlsTol seems to be important, if set to 0.01, then often does not converge
# tolerance=1e-4,
# # msTol=1e-1, minScale=1e-1, .relStep=1e-7,
# returnObject=TRUE, # allow the return of the fit when the max iter is reached
# maxIter=5000, nlsMaxIter=50, opt="nlminb", msVerbose=T)), silent=FALSE
# )
# new.theta=coef(fit.1)
# new.theta=c(new.theta, theta["sigma.sq"]) # has to be on two lines, otherwise names are lost
# if (any(support<new.theta["logc"]) | any(support>new.theta["logd"])) {
# #some support is outside new.theta, revert new.theta to theta
# new.theta=theta
# if (verbose>=2) cat("initial theta for step 3:", new.theta, "\n")
# } else {
# if (verbose>=2) cat("initial theta for step 3:", new.theta, "\n")
# }
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.