#' profile
#'
#' @description
#' Performs a profile using residual sum of squares, fixes some parameters for a range of values
#' and then estimate the others
#'
#' @param fitted: an \code{aspic} object
#' @param which: \code{character} giving the parameters to do the profile for, i.e. to fix.
#' @param range; \code{numeric} relative values by which to vary parameter values, default seq(0.5,1.5,length.out=21).
#' @param fn: \code{function} that gives values to be profiled.
#' @param run: \code{logical} if \code{TRUE} then returns profile, otherwise it just sets the control object-
#'
#' @return a \code{data frame} with results turned by \code{fn} by values in \code{which}.
#' @seealso \code{\link{biodyn},\link{fit}}
#'
#' @export
#' @docType methods
#' @rdname profile
#'
#' @examples
#' /dontrun{
#' data(asp)
#' res=profile(asp,which="msy",range=seq(0.75,1.25,length.out=21))
#' ggplot(res)+geom_line(aes(k,rss))
#'
#' fn =function(x) cbind( data.frame(model.frame(params(x)),rss=sum(diags(x)$residual^2,na.rm=T)))
#' res=profile(asp,which="msy",length.out=3,range=c(0.5,1.1),fn=fn)
#'
#' fn=function(x) ddply(x@diags, .(name), with, sum(residual^2,na.rm=T)/sum(count(!is.na(residual))))
#'
#' }
setMethod("profile", signature(fitted="aspic"),
function(fitted,which,
range=seq(0.5,1.5,length.out=21),
fn =function(x) cbind(model.frame(params(x)),
model.frame(refpts(x)),
model.frame(x@objFn)[,-3],
model.frame(x@ll[,2]/(x@ll[,3]-1))),
run=TRUE,...){
if (dims(fitted)$iter>1) stop("can only be done for a single iter")
if (length(range)==1) range=c(range,2-range)
if (dim(fitted@control)[3]==1){
sq=list(range)
sq=do.call("expand.grid",sq[rep(1,length(which))])
names(sq)=which
fitted@control=propagate(fitted@control,dim(sq)[1])
#fix min bounds for msy & k
for (i in which){
fitted@control[i,"val"]=params(fitted)[i,]*sq[,i]
fitted@control[i,"min"]=min(fitted@control[i,"val"])*.9
fitted@control[i,"max"]=max(fitted@control[i,"val"])*1.1
fitted@control[i,"fit"]=0
if (("k" == i) & all(c("msy","k")%in%dimnames(fitted@control)$params))
fitted@control["msy","min"]=min(fitted@control["k","min"],fitted@control["msy","min"])*0.9
if (("msy" == i) & all(c("msy","k")%in%dimnames(fitted@control)$params))
fitted@control["msy","min"]=min(fitted@control["k","min"],fitted@control["msy","min"])*0.9
}
if (!run) return(fitted)
fitted=fit(fitted)
fitted=fn(fitted)}
else
fitted@control=profileGrid(fitted@control,which,range)
return(fitted)})
refpts<-function(object){
msy=fmsy(object)*bmsy(object)
dimnames(msy)$params="msy"
rbind(msy,
bmsy(object),
fmsy(object))}
setMethod('profile', signature(fitted='aspics'),
function(fitted,which,
range=seq(0.5,1.5,length.out=21),
fn =function(x) cbind(model.frame(params(x)),
model.frame(refpts(x)),
model.frame(x@objFn)[,-3],
t(unlist(x@ll)[1,drop=T])),
run=TRUE,
.multicombine=T,.maxcombine=10,.packages=c("aspic","plyr")){
if (run) .combine=rbind.fill else .combine=list
res=foreach(i=names(fitted), .combine=.combine,
.multicombine=.multicombine,
.maxcombine =.maxcombine,
.packages =.packages) %dopar% {
profile(fitted[[i]],which=which,range=range,fn=fn,run=run)}
if (!run) {
res =aspics(res)
names(res)=names(fitted)}
else {
ns =ldply(fitted, function(x) length(unique(index(x)[,"name"])))[,2]
ns =ns*length(range)*length(which)
ns =rep(seq(length(ns)),ns)
res=cbind(.id=names(fitted)[ns],res)
}
res})
#' @description
#' Outputs as a data.frame a summary of parameters and RSS etc by data component
#' generated when do a profile
#'
#' @param x: an \code{aspic} object
#'
#' @return a \code{data frame} with results by data component.
#' @seealso \code{\link{biodyn},\link{profile} \link{fit}}
#'
#' @export
#' @docType methods
#' @rdname fnPiner
#'
#' @examples
#' /dontrun{
#' data(asp)
#' dcK=profile(asp,which=c("k"),range=seq(0.2,2.0,length.out=21),fn=fnPiner)
#' ggplot(dcK)+geom_line(aes(k,value,group=variable,col=variable))+
#' theme_ms(12,legend.position="bottom")+
#' ylab("Residual Sum of Squares")+xlab("K")
#'
#' }
fnPiner=function(x) {
res=cbind(model.frame(params(x),drop=T)[,-(dim(params(x))[1]+1)],
model.frame(refpts(x),drop=T)[,-c(1,4)],
model.frame(FLQuants(stock =stock( x)[,ac(range(x)["maxyear"])]%/%bmsy(x),
harvest=harvest(x)[,ac(range(x)["maxyear"])]%/%fmsy(x)),drop=T)[,-1],
model.frame(x@ll[,"ss"])[,seq(dim(x@ll)[1])]
)
nms=dimnames(params(x))$params
res=melt(res,id=c(dimnames(refpts(x))$refpts[-1],nms,c("stock","harvest")))
res=transform(res,name=unique(index(x)$name)[res$variable])
res}
# ### debugging stuff
# data(bd)
# fitted=biodyn(factor("pellat"),params(bd),catch=catch(bd))
# cpue=rlnorm(1,log(stock(bd)),.2)[,-60]
# setParams(fitted) =cpue
#
#
# attach(list(length.out=11, range=0.5, ci=c(0.25, 0.5, 0.75, 0.95),
# plot=TRUE,fixed=c()))
# which="r"
# fixed=c("p","b0")
# ###
# rtn=profile(swon[[1]],which="k",fixed="b0",length.out=31,range=c(.75,1.5))
# ggplot(rtn)+geom_line(aes(msy,rss))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.