Nothing
#####
mcMAM<-
function(EDdata, ncomp=-1, addsigma=0, iflog=TRUE,
nsim=5e4, inis=list(), control.args=list()) {
UseMethod("mcMAM")
} #
### 2023.09.10.
mcMAM.default<-
function(EDdata, ncomp=-1, addsigma=0, iflog=TRUE,
nsim=5e4, inis=list(), control.args=list()) {
### Stop if not.
stopifnot(ncol(EDdata)==2L, nrow(EDdata)>=5L,
all(EDdata[,2L,drop=TRUE]>0),
length(ncomp)==1L, ncomp %in% c(-1L,-2L,-3L,-4L),
length(addsigma)==1L, is.numeric(addsigma),
length(iflog)==1L, is.logical(iflog),
length(nsim)==1L, is.numeric(nsim), nsim>=100L, nsim<=2e5,
is.list(inis), is.list(control.args),
all(names(inis) %in% c("p","gamma","mu","sigma")),
all(names(control.args) %in% c("w","m","nstart")))
###
ed1<-as.numeric(EDdata[,1L,drop=TRUE])
sed1<-as.numeric(EDdata[,2L,drop=TRUE])
nED<-nrow(EDdata)
if (ncomp %in% c(-1L,-2L)) {
inverse <- FALSE
} else {
inverse <- TRUE
} # end if.
###
if (iflog==TRUE && any(ed1<=0)) {
stop("Error: minus(zero) ED canot be logged!")
} # end if
### Bounds of dose population.
rangeED<-range(ed1)
if(all(ed1>0)) {
lowerGamma<-rangeED[1L]*0.999
upperGamma<-rangeED[2L]*1.001
} else if (all(ed1<=0)) {
lowerGamma<-rangeED[1L]*1.001
upperGamma<-rangeED[2L]*0.999
} else {
lowerGamma<-rangeED[1L]*1.001
upperGamma<-rangeED[2L]*1.001
} # end if
###
### Default inis.
args.inis<-list("p"=0.5,"gamma"=quantile(ed1,probs=0.25,names=FALSE),
"mu"=quantile(ed1,probs=0.5,names=FALSE),
"sigma"=ifelse(iflog==TRUE,0.3,3.0))
args.inis[names(inis)]<-inis
stopifnot(args.inis[["p"]]>0.0, args.inis[["p"]]<1.0,
args.inis[["gamma"]]>lowerGamma,
args.inis[["gamma"]]<upperGamma,
args.inis[["mu"]]>lowerGamma,
args.inis[["mu"]]<upperGamma,
args.inis[["sigma"]]>0.0,
args.inis[["sigma"]]<ifelse(iflog==TRUE,5.0,50.0))
###
### Default arguments of slice sampling.
args.control<-list(w=1, m=-100, nstart=1L)
args.control[names(control.args)]<-control.args
stopifnot(args.control[["w"]]>=1e-2,
args.control[["w"]]<=1e3,
args.control[["m"]]<=1e9,
args.control[["nstart"]]>=1L,
args.control[["nstart"]]<=1000L)
###
w<-args.control$w
m<-args.control$m
nstart<-args.control$nstart
inis<-if (ncomp==-1L) {
c(args.inis[["p"]],
args.inis[["gamma"]],
args.inis[["sigma"]])
} else {
c(args.inis[["p"]],
args.inis[["gamma"]],
args.inis[["mu"]],
args.inis[["sigma"]])
} # end if
iflag<-0
chains<-matrix(0,nrow=nsim,ncol=ifelse(ncomp %in% c(-1L,-3L),3L,4L))
###
if (ncomp %in% c(-1L,-3L)) {
res <- .Fortran("mcMAM3",as.integer(nED),as.integer(nsim),as.double(ed1),
as.double(sed1),as.double(addsigma),as.double(inis),
as.integer(iflog),as.integer(nstart),as.double(w),
as.double(m),chains=as.double(chains),as.integer(inverse),
iflag=as.integer(iflag),PACKAGE="numOSL")
} else if (ncomp %in% c(-2L,-4L)) {
res <- .Fortran("mcMAM4",as.integer(nED),as.integer(nsim),as.double(ed1),
as.double(sed1),as.double(addsigma),as.double(inis),
as.integer(iflog),as.integer(nstart),as.double(w),
as.double(m),chains=as.double(chains),as.integer(inverse),
iflag=as.integer(iflag),PACKAGE="numOSL")
} # end if.
###
if (res$iflag!=0) {
niter<-sum(res$chains[seq(nsim)]>0)
stop(paste("Error: the simulation failed at the ",
niter+1L,"th iteration!",sep=""))
} # end if
###
if (ncomp==-1L) {
myNAMES <- c("Prop","MAM3.De","Sigma")
} else if (ncomp==-2L) {
myNAMES <- c("Prop","MAM4.De","Mu","Sigma")
} else if (ncomp==-3L) {
myNAMES <- c("Prop","MXAM3.De","Sigma")
} else if (ncomp==-4L) {
myNAMES <- c("Prop","MXAM4.De","Mu","Sigma")
} # end if.
chains<-matrix(res$chains,ncol=ifelse(ncomp %in% c(-1L,-3L),3L,4L),
dimnames=list(NULL,myNAMES))
###
if (iflog==TRUE) {
if (ncomp==-1L) {
chains[,2L]<-exp(chains[,2L])
model <- "MAM3"
} else if (ncomp==-3L) {
chains[,2L]<-exp(-chains[,2L])
model <- "MXAM3"
} else if (ncomp==-2L) {
chains[,2L:3L]<-exp(chains[,2L:3L])
model <- "MAM4"
} else if (ncomp==-4L) {
chains[,2L:3L]<-exp(-chains[,2L:3L])
model <- "MXAM4"
} # end if
} else {
if (ncomp==-1L) {
model <- "MAM3"
} else if (ncomp==-2L) {
model <- "MAM4"
} else if (ncomp==-3L) {
chains[,2L]<- -chains[,2L]
model <- "MXAM3"
} else if (ncomp==-4L) {
chains[,2L:3L]<- -chains[,2L:3L]
model <- "MXAM4"
} # end if
} # end if.
###
output<-list("EDdata"=EDdata,
"addsigma"=addsigma,
"model"=model,
"iflog"=iflog,
"nsim"=nsim,
"chains"=chains)
class(output)<-"mcAgeModels"
invisible(output)
###
} # end function mcMAM.
#####
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.