R/MCpriorIntFun.pb.r

In BMAmevt: Multivariate Extremes: Bayesian Estimation of the Spectral Measure

##' Wrappers for  \code{\link{MCpriorIntFun}} with argument
##' \code{prior=prior.pb} or \code{prior=prior.nl} 
##'
##' @title Generic Monte-Carlo integration  under the prior distribution in the PB and NL  models.
##' @inheritParams MCpriorIntFun
##' @param Hpar Hyper-parameters for the PB prior (in \code{MCpriorIntFun.pb}) or the  NL prior (\code{MCpriorIntFun.nl}). See
##' \code{\link{pb.Hpar}} and \code{\link{nl.Hpar}} for the required
##' formats.
##' @param dimData Only for the PB model: The dimension of model's \emph{sample} space. The PB parameter  space is of dimension \code{choose(dimData,2)+1}. The NL model implemented here is restricted to three-dimensional sample spaces. 
##' @return The list returned by function
##' \code{\link{MCpriorIntFun}}.
##' @seealso \code{\link{MCpriorIntFun}}
##' @export
MCpriorIntFun.pb <-
function(Nsim=200,
         Hpar = get("pb.Hpar"),
         dimData=3,
         FUN=function(par, ... ){
           as.vector(par)
         },
         store=TRUE,
         show.progress = floor(seq(1, Nsim, length.out = 20 ) ),
         Nsim.min=Nsim,
         precision = 0,
         ...)
  {
    mcres <- MCpriorIntFun(Nsim=Nsim,
                           prior=prior.pb,
                           Hpar=Hpar,
                           dimData=dimData,
                           FUN=FUN,
                           show.progress=show.progress,
                           Nsim.min=Nsim.min,
                           precision=precision,
                           ...)
    return(mcres)


  }    
## ############ intialize  ############

##     start.time=proc.time()
##     not.finite=0
##     param = rparameter.prior.pb(n=1,p=p, Hpar=Hpar)
##     temp.res=FUN(param,...)
##     dim.res=dim(temp.res)
    
##     if(is.null(dim.res) || (sum(dim.res!=1) ==1)  )
##       {
##         emp.mean=rep(0,length(temp.res))
##       }
##     else
##       {
##         store=FALSE
##         emp.mean=array(0,dim=dim.res)
##       }
##     emp.variance= emp.mean
##     emp.variance.unNorm=emp.variance

##     if(store)
##       {
##         stored.vals=matrix(0,nrow=Nsim,ncol=length(emp.mean))
##       }
      
## ################ start MC ########
##     n=1
##     while((n<=Nsim) &&
##           ( (n<=Nsim.min) ||
##            (max( sqrt(emp.variance/(n-1))/emp.mean) > precision) )
##           )
##       {
##         ## show progression      
##         if(any(n==be.patient) || n==50 || n==100 || n==200)
##           {
##             cat(paste((n-1), "iterations done", "\n", sep = " " ))
##           }


##         flag=TRUE
##         count = 0
##         while(flag & (count<=10))
##           {
##             param = rparameter.prior.pb(n=1,p=p, Hpar=Hpar)
##             temp.res=FUN(param,...)
##             flag = (any(sapply(as.vector(temp.res),
##                                function(x){ ! is.finite(x) } ) ) )
##             if(flag)
##               {
##                 not.finite = not.finite+1
##               }
##             count = count+1
##           }
##         if(flag)
##           stop("more than 10 non finite values produced")
        
##         cur.res=temp.res ##FUN(alpha=param[1],beta=param[2:(ll+1)],...)
        
##         new.emp.mean=emp.mean+1/n*(cur.res-emp.mean)
        
##         emp.variance.unNorm=emp.variance.unNorm +
##           (cur.res-new.emp.mean)* (cur.res- emp.mean)
##         emp.variance = emp.variance.unNorm/(n-1)
##         emp.mean = new.emp.mean
        
##         if(store)
##           {
##             stored.vals[n,]= as.vector(cur.res)
##           }
##         n=n+1
##       }
## ######### end MC #########
##     end.time = proc.time()
##     elapsed=end.time-start.time
##     print(elapsed)
##     if(store)
##       {
##         return(list( stored.vals=stored.vals[1:(n-1), ],
##                     elapsed=elapsed,
##                     nsim = n-1,
##                     emp.mean=emp.mean,
##                     emp.stdev=sqrt(emp.variance),
##                     est.error=sqrt(emp.variance/(n-1)),
##                     not.finite = not.finite))
##       }
##     else
##       {
##         return(list(  nsim=n-1,
##                     elapsed=elapsed,
##                     emp.mean=emp.mean,
##                     emp.stdev=sqrt(emp.variance),
##                     est.error=sqrt(emp.variance/(n-1)),
##                     not.finite = not.finite ))
                  
##       }
##   }