R/bayesian.dw.R

Defines functions bdw.mc dw.bay.plot.density dw.HPDinterval l.tot.dw summary.bdw dw.info.Criteria dw.BF bay.laplace.calc bdw.plot plot.bdw par.bayesian.tot.dw MH.tot.dw posterior.tot.dw prior.tot.dw loglik.qbeta.dw loglik.beta.dw loglik.q.dw loglik.sim.dw

Documented in bdw.mc plot.bdw summary.bdw

# simple loglikelihood
loglik.sim.dw <-  function(lik.par , x  ,    lq.extra ,  lbeta.extra , logit)
{
  l.l = length(lik.par)
  if( l.l != 2) stop('error in the number of parameters! ',l.l,' <> T = ',2)
  q    = lik.par[1] * lq.extra
  beta = lik.par[2] * lbeta.extra

  if (all(beta >= 0) && all(q<= 1) && all(q >= 0) ){
    loglik = sum(log(q^(x^beta)-q^((x+1)^beta)))
  }else{
    loglik = NA
  }
  return(loglik)
}
# loglikelihood - regression on q
loglik.q.dw  <-  function(lik.par , x , y , lq.extra ,  lbeta.extra  , logit )
{
  l.like.par = length(lik.par)
  if( l.like.par != ncol(x)+1) stop('error in the number of parameters!',l.like.par,'<> T=',ncol(x)+1)
  beta = lik.par[ l.like.par] * lbeta.extra
  theta= lik.par[-l.like.par] * lq.extra
  if(logit){
    q    = exp (x %*% theta-log(1+exp(x %*% theta)) )
  }else{
    q    = exp (- exp ( x %*% theta ))
  }

  if (all(beta >= 0) && all(q <= 1) && all(q >= 0) ){
    loglik = sum(log(q^(y^beta)-q^((y+1)^beta)))
  }else{
    loglik = NA
  }
  return(loglik)
}
#  loglikelihood - regression on beta
loglik.beta.dw <- function(lik.par , x , y , lq.extra ,  lbeta.extra , logit)
{
  l.l = length(lik.par)
  if(l.l  != ncol(x)+1) stop('error in the number of parameters!You:',l.l,' <> T = ',ncol(x)+1)
  q      = lik.par[ 1] * lq.extra    #par[ length(par)]
  theta  = lik.par[-1] * lbeta.extra #par[-length(par)]
  beta   = exp(x %*% theta)
  #if(q %in% theta) print('Q in theta!')

  if (all(beta >= 0) && all(q <= 1) && all(q >= 0) ){
    loglik = sum(log(q^(y^beta)-q^((y+1)^beta)))
  }else{
    loglik = NA
  }
  return(loglik)
}
#  loglikelihood - regression on both beta and q
loglik.qbeta.dw<- function(lik.par , x , y , lq.extra ,  lbeta.extra , logit)
{
  ncx        = ncol(x)
  l.l        = length(lik.par)
  if( l.l != 2*ncx) stop('error in the number of parameters!',l.l,' <> T = ',2*ncol(x) )
  theta.q    = lik.par[1: ncx ]     * lq.extra
  theta.beta = tail(lik.par , ncx ) * lbeta.extra

  x1 = x #+ xxx1
  x2 = x #+ xxx2

  if(all(x[,1]==1)){
    x1[,1]=1
    x2[,1]=1
  }

  if(logit){
    q          = exp(x2 %*% theta.q-log(1+exp(x %*% theta.q)))
  }else{
    q          = exp(-exp(x2 %*% theta.q))
  }
  beta         = exp(x1 %*% theta.beta)

  if (all(beta >= 0) && all(q <= 1) && all(q >= 0) ){
    loglik = sum(log(q^(y^beta)-q^((y+1)^beta)))
  }else{
    loglik = NA
  }
  return(loglik)
}
###########   end of likelihoods ##################
#  prior distribution
prior.tot.dw <-  function(pr.par      , q.par , b.par      ,
                          dist.q      , dist.b             ,
                          lq , lb     , l.par , dist.l     ,
                          penalized   , fixed.l           ,
                          para.q      , para.b            ,
                          lq.extra    ,
                          lbeta.extra ,
                          logit       ,
                          x = x       ,
                          y = y       ,
                          prior.function = prior.function)
{
  #------ Lambda prior parameters if needed
  par.length   = length(pr.par)

  if((penalized == TRUE) && (fixed.l <= 0) ){
    if(par.length != (lb+lq+1)) stop('error in the number of parameters!',par.length,' <> T = ',(lb+lq+1))
    #lp           = pr.par[(lb+lq+1):par.length]
    lp           = pr.par[par.length]
    l.par1       = l.par[1]
    l.par2       = l.par[2]
    lambda.prior =  sum (dist.l(lp ,  l.par1  ,   l.par2  ,   log = TRUE ))
  }else if ( (penalized == TRUE) && (fixed.l > 0) ){
    if(par.length != (lb+lq)) stop('error in the number of parameters!',par.length,' <> T = ' ,(lb+lq))
    lambda.prior = 0
    lp           = fixed.l
  }else{
    if(par.length != (lb+lq)) stop('error in the number of parameters!',par.length,' <> T = ',(lb+lq))
    #lp           = rep(1,lp+lq)
    lp           = 1
    lambda.prior = 0
  }
  #------ Q prior parameters
  q.par1 =               q.par[1]
  q.par2 =    sqrt(lp) * q.par[2]
  #q.par2 =  head(lp,lq)* q.par[2]
  #------ Beta prior parameters
  b.par1 =                 b.par[1]
  b.par2 =     sqrt(lp) *  b.par[2]
  #b.par2 =  tail(lp,lb) *  b.par[2]
  #------ model parameters
  qp     =     (pr.par[1:lq]           )
  bp     =     (pr.par[ (lq+1):(lq+lb)])
  #------ priors
  theta.prior  = sum (dist.q(qp ,  q.par1  ,   q.par2 ,    log = TRUE))
  gamma.prior  = sum (dist.b(bp ,  b.par1  ,   b.par2  ,   log = TRUE))
  prior =  ( theta.prior + gamma.prior + lambda.prior)
  return(prior)
}
#  posteriot distribution
posterior.tot.dw <- function(par     ,
                             x       ,  y      ,
                             lq.extra          ,
                             lbeta.extra       ,
                             para.q  , para.b  ,
                             prior.function    ,  ####### END OF PRIOR ####
                             dist.q  , q.par   ,
                             dist.b  , b.par   ,
                             dist.l  , l.par   ,
                             penalized  , fixed.l ,
                             logit , iter, jeffry,
                             ...
)
{
  iterTresh = 1;jm=5
  if ((para.b == FALSE) && (para.q == TRUE)) {
    loglik   = loglik.q.dw   (lik.par = par[1:(ncol(x) + 1)],
                              x = x , y = y ,
                              lq.extra = lq.extra ,
                              lbeta.extra = lbeta.extra,
                              logit = logit)
    if(jeffry==FALSE){
      prior    = prior.function(pr.par = par,
                                q.par  = q.par   , b.par  = b.par   ,
                                dist.q = dist.q  , dist.b = dist.b  ,
                                #lb = 1           , lq = ncol(x)-1  ,
                                lb = 1           , lq = ncol(x)     ,
                                l.par = l.par    , dist.l = dist.l  ,
                                penalized =  penalized , fixed.l = fixed.l, ###
                                para.q = para.q , para.b = para.b     ,
                                lq.extra = lq.extra , lbeta.extra = lbeta.extra,
                                logit = logit , x = x , y = y
      )
    }else if (jeffry==TRUE && (iter>jm*iterTresh) && (iter %% iterTresh == 0)){
      #------- jeffry Prior
      fit = optim(par =  par[1:(ncol(x) + 1)],
                  fn          = loglik.q.dw,
                  x = x , y   = y ,
                  lq.extra    = lq.extra ,
                  lbeta.extra = lbeta.extra,
                  logit       = logit,
                  #method = 'BFGS',
                  control=list("fnscale"=-1),
                  hessian=TRUE)
      fisher_info = solve(-fit$hessian)
      prior       = log(sqrt(det(fisher_info)+10^-16))
    }else{
      prior = 0
    }
    #------
    post     = loglik + prior
  }else if ((para.b == TRUE ) && (para.q == FALSE)) {
    loglik   = loglik.beta.dw(lik.par = par[1:(ncol(x) + 1)], x = x, y = y ,
                              lq.extra = lq.extra ,lbeta.extra = lbeta.extra,
                              logit = logit)
    if (jeffry==FALSE){
      prior    = prior.function(pr.par = par,
                                q.par  = q.par   , b.par  = b.par  ,
                                dist.q = dist.q  , dist.b = dist.b ,
                                dist.l = dist.l  , l.par  = l.par  ,
                                #lb = ncol(x)-1   , lq = 1  ,
                                lb = ncol(x)     , lq = 1          ,
                                penalized = penalized , fixed.l = fixed.l,
                                para.q = para.q , para.b = para.b     ,
                                lq.extra = lq.extra , lbeta.extra = lbeta.extra,
                                logit = logit , x = x , y = y)
    }else if(jeffry==TRUE && (iter>jm*iterTresh) && (iter %% iterTresh == 0)){
      #------- jeffry Prior
      fit = optim(par =  par[1:(ncol(x) + 1)],
                  fn          = loglik.beta.dw,
                  x = x , y   = y ,
                  lq.extra    = lq.extra ,
                  lbeta.extra = lbeta.extra,
                  logit       = logit,
                  #method = 'BFGS',
                  control=list("fnscale"=-1),
                  hessian=TRUE)
      fisher_info = solve(-fit$hessian)
      prior       = log(sqrt(det(fisher_info)+10^-16))
      #------
    }else{
      prior = 0
    }
    post     = loglik + prior
  }else  if ((para.b == TRUE) && (para.q == TRUE)) {
    loglik   =  loglik.qbeta.dw(lik.par  = par[1:(2*ncol(x))], x = x, y = y ,
                                lq.extra = lq.extra ,lbeta.extra = lbeta.extra,
                                logit = logit)
    if(jeffry==FALSE){
      prior    =  prior.function (pr.par    = par    ,
                                  q.par  = q.par  , b.par  = b.par  ,
                                  dist.q = dist.q , dist.b = dist.b ,
                                  dist.l = dist.l , l.par  = l.par  ,
                                  lb = ncol(x)    , lq = ncol(x)    ,
                                  penalized = penalized ,  fixed.l = fixed.l,
                                  para.q = para.q , para.b = para.b     ,
                                  lq.extra = lq.extra , lbeta.extra = lbeta.extra,
                                  logit = logit , x = x , y = y)
    }else if(jeffry==TRUE && (iter>jm*iterTresh) && (iter %% iterTresh == 0)){
      #------- jeffry Prior
      fit = optim(par =  par[1:(2*ncol(x))],
                  fn          = loglik.qbeta.dw,
                  x = x , y   = y ,
                  lq.extra    = lq.extra ,
                  lbeta.extra = lbeta.extra,
                  logit       = logit,
                  #method = 'BFGS',
                  control=list("fnscale"=-1),
                  hessian=TRUE)
      fisher_info = solve(-fit$hessian)
      prior       = log(sqrt(det(fisher_info)+10^-16))
      #------
    }else{
      prior = 0
    }
    post   = loglik + prior
  }else{
    loglik   =  loglik.sim.dw (lik.par  = par[1:2], x = y  ,
                               lq.extra = lq.extra ,lbeta.extra = lbeta.extra,
                               logit = logit)
    if(jeffry==TRUE){
      prior    =  prior.function(pr.par    = par                ,
                                 q.par  = q.par  , b.par  = b.par  ,
                                 dist.q = dist.q , dist.b = dist.b ,
                                 dist.l = dist.l , l.par = l.par   ,
                                 lb = 1          , lq = 1 ,
                                 penalized = penalized ,  fixed.l = fixed.l,
                                 para.q = para.q , para.b = para.b     ,
                                 lq.extra = lq.extra , lbeta.extra = lbeta.extra,
                                 logit = logit , x = x , y = y)

    }else if(jeffry==TRUE && (iter>jm*iterTresh) && (iter %% iterTresh == 0)){
      #------- jeffry Prior
      # fit = optim(par =  par[1:2],
      #             fn          = loglik.sim.dw,
      #             x = x , y   = y ,
      #             lq.extra    = lq.extra ,
      #             lbeta.extra = lbeta.extra,
      #             logit       = logit,
      #             #method = 'BFGS',
      #             control=list("fnscale"=-1),
      #             hessian=TRUE)
      # fisher_info = solve(-fit$hessian)
      # prior       = log(sqrt(det(fisher_info)+10^-16))
      prior         = log(((-log(par[1]))^(1/par[2]))/par[2])
    }else{
      prior = 0
    }
    #------
    post   = (loglik + prior)
  }
  return(post)

}






# Reversible jumps & Metropolis-Hasting sampling
MH.tot.dw <- function(formula, data, startvalue      ,
                      para.q  , para.b               ,
                      dist.q  , dist.b               ,
                      q.par   , b.par                ,
                      dist.l  , l.par                ,
                      fixed.l , logit                ,
                      cov.m   , penalized            ,
                      v.scale , iterations           ,
                      lq.extra, lbeta.extra          ,
                      prior.function                 ,
                      RJ = FALSE                     ,
                      jeffry = FALSE
)
{
  if ( para.q  || para.b ) {
    call <- match.call()
    mf <- match.call(expand.dots = TRUE)
    m  <- match(c("formula", "data"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    mt <- attr(mf, "terms")
    y  <- model.response(mf, "numeric")
    x  <- model.matrix(mt, mf, contrasts)
    if (para.b == TRUE) {lb = ncol(x)}else{lb = 1}
    if (para.q == TRUE) {lq = ncol(x)}else{lq = 1}
  }else{
    y  = as.vector(data)
    x  = NA
    lb = lq = 1
    RJ = FALSE
  }

  # if( is.na(startvalue) ){
  #   if(para.q || para.b){
  #     o = DWreg::dw.reg(formula = formula,data = data,
  #                       para.q1 = logit * para.q , para.q2 = para.q * !logit,
  #                       para.beta = para.b)
  #     startvalue = o$coefficients
  #     if(logit * para.q){
  #       startvalue = o$tTable$estimate[,1]
  #     }else if (!logit * para.q){
  #       startvalue = o$tTable[,1]
  #     }
  #     if(penalized){
  #       startvalue = c(startvalue,max(cov(x,y)) )
  #     }
  #   }else{
  #     if(is.matrix(data)) stop('data must have just one column!')
  #     o = DWreg::dw.parest(data = data)
  #     startvalue = c(o$q,o$beta)
  #   }
  # }

  len.par = length(startvalue)
  #### This is for RJ
  if(RJ){
    #message('Reversible-Jump in action ... ')
    penalized = FALSE ; fixed.l = -1
  }
  total.a.par = lb + lq
  prob   = rep(.5 , total.a.par) # initial  probability of existing variables
  rjmu   = rep( 0 , total.a.par)
  rjsig  = rep(.5 , total.a.par)
  model.chain       = c()
  chain             = array(dim = c(iterations + 1 , len.par));
  chain[1,]         = startvalue
  ####### Start iterations  #######
  pr2   = posterior.tot.dw( par = chain[1 , ]  ,
                            para.q = para.q , para.b = para.b ,
                            q.par = q.par   , b.par = b.par   ,
                            dist.q = dist.q , dist.b = dist.b ,
                            dist.l =  dist.l, l.par  = l.par  ,
                            x = x   , y = y , logit = logit   ,
                            lq.extra    = lq.extra            ,
                            lbeta.extra = lbeta.extra         ,
                            penalized   = penalized           ,
                            prior.function = prior.function   ,
                            fixed.l     = fixed.l , iter = 1  ,
                            jeffry = jeffry
  )
  state         = startvalue
  functionat    = pr2


  ##########################
  if(length( startvalue )>1 )
  {
    fit = optim(par =  startvalue,
                fn = posterior.tot.dw,
                para.q = para.q , para.b = para.b ,
                q.par  = q.par  , b.par = b.par   ,
                dist.q = dist.q , dist.b = dist.b ,
                dist.l = dist.l , l.par  = l.par  ,
                x = x   , y = y,  logit = logit   ,
                lq.extra = lq.extra               ,
                lbeta.extra = lbeta.extra         ,
                penalized = penalized             ,
                prior.function = prior.function   ,
                fixed.l = fixed.l, iter =1        ,
                jeffry = jeffry                   ,
                control=list("fnscale"=-1),hessian=TRUE)

    fisher_info = solve(-fit$hessian)
    prop_sigma = sqrt(diag(fisher_info)+10^-16)
    prop_sigma = diag(prop_sigma)
  }else{
    prop_sigma<- 1+ chain[1 , ]/2
  }
  mytemp = make.positive.definite(prop_sigma)



  #  Proposal distribution
  proposalfunction.tot.dw <- function(par   , v.scale  ,
                                      para.q, para.b   ,
                                      s , i , cov.m    ,
                                      lb , lq ,tmp){
    length.par = length(par)
    sm = v.scale
    if (cov.m == 1) {
      if(i >= 100 &&
         i %% 100 == 0)
      {
        sm   =    dw.cov.matrix(x = s , i = i )
      }else{
        sm = tmp
      }
      ncs =    ncol(s)
      b1  =    MASS::mvrnorm(n = 1   , mu   = par*0 ,  Sigma = sm)
      #b2  =    rnorm  (n = 1   , mean = 0     ,  sd    = rep(1/ncs , ncs))
      prop.d = par+b1 #+ v.scale*b2

    }else if (cov.m == 2) {
      prop.d     =    par + runif(length.par ,  -abs(v.scale)/2 , abs(v.scale)/2 )
    }else if (cov.m == 3) {
      prop.d = par + rlaplace(n = length.par , mu =  0 ,sigma = v.scale )
    }else{
      prop.d     =    par + rnorm(length.par,  0               , v.scale )
    }
    if (length.par > lb + lq) {
      lp = (lb + lq + 1):(length.par)
      prop.d[lp] = abs(prop.d[lp])
    }
    return(list(prop.d=prop.d,sm=sm))
  }
  #########################
  # ---- CORE
  i = 1;  j = 1;  k = 0 ; ReAc = 0; error = 0 ;
  it.interval = round(iterations/100) ; interval = 0
  while (i <= iterations) {
    if (i %% it.interval == 0) {
      acc2 = apply(chain,2, function(x){1-mean(duplicated(x))} )
      acceptance =acc2
      acceptance = max (acceptance) * 100
      info <- paste(round(i/iterations*100),
                    '% done, Acceptance = ',round(acceptance,2),'%')
      cat('\r ',  info ,rep(' ',20))
    }
    if (i < 15) { #just to make sure that the process is not stuck
      vs = 10 + v.scale
    }  else  {
      vs =      v.scale

    }


    proposal.raw = proposalfunction.tot.dw(par = chain[i , ],
                                           para.q = para.q ,para.b = para.b,
                                           v.scale = vs   ,
                                           s = chain , i = i   ,
                                           cov.m  = cov.m      ,
                                           lb = lb , lq = lq,
                                           tmp = mytemp)
    proposal=proposal.raw$prop.d
    mytemp=proposal.raw$sm

    # To make sure that proposal is not acting on zero parmeters.
    if(RJ){
      zeros = which (chain[i , ] == 0)
      if(length(zeros) > 0) {
        proposal[zeros] = 0
      }
    }
    pr1 = posterior.tot.dw( par = proposal ,
                            para.q = para.q , para.b = para.b ,
                            q.par  = q.par  , b.par = b.par   ,
                            dist.q = dist.q , dist.b = dist.b ,
                            dist.l = dist.l , l.par  = l.par  ,
                            x = x   , y = y,  logit = logit   ,
                            lq.extra = lq.extra               ,
                            lbeta.extra = lbeta.extra         ,
                            penalized = penalized             ,
                            prior.function = prior.function   ,
                            fixed.l = fixed.l , iter =i       ,
                            jeffry = jeffry
    )

    probab = min(1,exp( pr1 - pr2 )) #* correction
    if (!is.na(probab)) {
      if ((runif(1) <= probab ) ) {
        chain[i + 1 , ]  = proposal
        ReAc [i + 1   ]  = 1
        #pr4              = pr2
        pr2              = pr1
      }else{
        chain[i + 1 , ] = chain[i , ]
        ReAc [i + 1   ] = 0
        #pr4             = pr1
      }

      ###### ----- This is a maximizing procedure ---- ###
      if (tail(functionat,1)  <  pr1 )  {
        state            = rbind(state , proposal)
        functionat       =  c(functionat,pr1)
        #cat('\r OK')
      }
      ###### END of Maximizing

      #----- START Riversible jump!
      sz = 1  # one sample each step!
      if(lb==1 && lq>1){
        interval = 1:lq
        r  = sample(x = interval , size = sz)
      }else if(lb>1 && lq==1){
        interval = 2:total.a.par
        r  = sample(x = interval , size = sz)
      }else if(lb>1 && lq>1){
        interval = 1:total.a.par
        r  = sample(x =  interval, size = sz)
      }else{
        RJ = FALSE
      }
      if(RJ){
        pr4 = pr2
        tr = chain[i + 1 , r]
        newchain    = as.vector(chain[i + 1 , ])
        if ( all(tr == 0) ){
          rv          = rnorm(length(r) , rjmu[r] , rjsig[r])
          newchain[r] = rv
          pr3 = posterior.tot.dw( par = newchain ,
                                  para.q = para.q , para.b = para.b ,
                                  q.par  = q.par  , b.par = b.par   ,
                                  dist.q = dist.q , dist.b = dist.b ,
                                  dist.l = dist.l , l.par  = l.par  ,
                                  x = x , y = y   , logit = logit   ,
                                  lq.extra = lq.extra               ,
                                  lbeta.extra = lbeta.extra         ,
                                  penalized = penalized             ,
                                  prior.function = prior.function   ,
                                  fixed.l = fixed.l, iter =i        ,
                                  jeffry = jeffry
          )
          #        num   =  pr3 + dnorm(rv , chain[i + 1 , r]    , chain[i + 1 , len.par], log = TRUE) + log(prob[r])
          num   =  pr3 + sum(dist.b(rv  , b.par[1], b.par[2], log = TRUE)) + sum(log(  prob[r]))
          den   =  pr4 + sum(dnorm (rv  , rjmu[r] , rjsig[r], log = TRUE)) + sum(log(1-prob[r]))
        }else{
          rv          = 0
          newchain[r] = rv
          pr3 = posterior.tot.dw( par = newchain ,
                                  para.q = para.q , para.b = para.b ,
                                  q.par  = q.par  , b.par = b.par   ,
                                  dist.q = dist.q , dist.b = dist.b ,
                                  dist.l = dist.l , l.par  = l.par  ,
                                  x = x  , y = y  , logit = logit   ,
                                  lq.extra = lq.extra               ,
                                  lbeta.extra = lbeta.extra         ,
                                  penalized = penalized             ,
                                  prior.function = prior.function   ,
                                  fixed.l = fixed.l, iter =i        ,
                                  jeffry = jeffry
          )
          #den   =  pr4 + dnorm(tr , chain[i + 1 , r]  , chain[i + 1 , len.par] , log = TRUE) + log(prob[r])
          num   =  pr3 + sum(dnorm (tr , rjmu[r], rjsig[r], log = TRUE)) + sum(log(1-prob[r]))
          den   =  pr4 + sum(dist.b(tr ,b.par[1], b.par[2], log = TRUE)) + sum(log(  prob[r]))
        }
        if(runif(1) <= min(1,exp(num-den)) ){
          chain[i + 1 , r] = rv
          pr2              = pr3
        }else{
          chain[i + 1 , r] = tr
        }
        if(chain[i + 1 , r]==0) adre = -1 else adre = 1
        model.chain[i]   = adre*r
      }
      #----- END Riverseible jump!
      i = i + 1; j = 1
    }else{
      j = j + 1; k = k + 1; #cat('\r E:',k)
      if ( j > 1500 ) stop('error in model! Maybe you need a different set of initial values or model parameters.',call. = 0)
    }

  }
  # ---
  cat('\r\n')

  if (k > 0) {cat('\n There are ',k,' ignored values in the process!')} ; error = k
  return(list( chain = chain     , RejAcc = ReAc    ,
               error = error     , cov.m = cov.m    ,
               minf = functionat , minState = state ,
               x = x             , y = y            ,
               lb = lb           , lq = lq          ,
               logit = logit     , RJ = RJ          ,
               model.chain = checkin(model.chain,interval ) )
  )
}



# semi final function (after is the final algorithm)
par.bayesian.tot.dw <- function(formula = NA,
                                data,
                                initials    , burn.in         ,
                                q.par       , b.par           ,
                                dist.q      , dist.b          ,
                                para.q      , para.b          ,
                                dist.l      , l.par           ,
                                fixed.l     , logit           ,
                                v.scale     , iterations      ,
                                cov.m       ,
                                sampling = c('indp','syst','bin'),
                                ###
                                penalized                        ,
                                lq.extra = 1 , lbeta.extra = 1   ,
                                prior.function = prior.tot.dw,
                                RJ = FALSE, jeffry = FALSE,
                                ...
)
{

  #library('MASS')
  if(RJ & (para.q || para.b)){
    #message('Reversible-Jump in action ... ')
    penalized = FALSE ; fixed.l = -1
  }else{
    RJ = FALSE
  }
  cat('\n\n============================== Sampler configuration ============================== \n')
  cat('Iterations:',iterations,'\t|', 'Data:' ,is.data.frame(data)       ,	 '\t \t|','Length of Initials:',length(initials),'\n')
  cat('RegQ:',para.q==1        ,	 '\t \t|'   ,'RegB:',para.b==1         ,	 '\t \t|','Formula:',is(formula,"formula"),'\n')
  cat('Logit:',logit==1        ,	 '\t \t|'   ,'Scale:',v.scale          ,	 '\t\t| Rev.Jumps:',RJ==1,'\n')
  cat('Penalized:',penalized==1,	 '\t|'      ,'Fixed.penalty:',fixed.l>0,   '\t| Jeffrey.Prior:',jeffry==1,'\n')
  cat('----------------------------------------------------------------------------------   \n' )
  cat('Proposal (1=Adaptive MH,2=Uniform,3=Laplace,>3=Gaussian):',cov.m,'\n' )
  cat('* if Adaptive MH is activated, then there is no need to set proposal scale.' )
  cat('----------------------------------------------------------------------------------   \n' )
  cat('Chain summary (bin=Burn-in, syst=Systematic, indp=Independent):',sampling,'\n' )
  cat('----------------------------------------------------------------------------------   \n' )
  cat('* If Penalized=TRUE then you need to set all distributions.                \n' )
  #cat('----------------------------------------------------------------------------------   \n' )
  cat('* If RJ=TRUE then Penalized is automatically set to FALSE and fixed.l diactivates. \n' )
  cat('__________________________________________________________________________________  \n' )


  # Start the clock!
  ptm <- proc.time()

  m.chain=MH.tot.dw(formula = formula , data=data,
                    startvalue = initials,
                    para.b = para.b , para.q = para.q  ,
                    dist.q = dist.q , dist.b = dist.b  ,
                    q.par  = q.par   , b.par = b.par   ,
                    dist.l =  dist.l, l.par  = l.par   ,
                    fixed.l = fixed.l                  ,
                    cov.m = cov.m                      ,
                    v.scale = v.scale                  ,
                    iterations = iterations            ,
                    penalized =  penalized             ,
                    lq.extra = lq.extra                ,
                    lbeta.extra = lbeta.extra          ,
                    prior.function = prior.function    ,
                    logit  = logit , RJ = RJ ,
                    jeffry = jeffry
  )


  unused.sample = round(iterations * burn.in)
  #----- SAMPLING
  if(sampling =='syst'){
    samp.seq   = seq(2,iterations,length.out = iterations-unused.sample)
    samp.seq   = round(samp.seq)
    chain      = m.chain$chain      [samp.seq,]
    model.chain= m.chain$model.chain[samp.seq ]
  }else if(sampling == 'indp'){
    samp.seq   = sample(x = 2:iterations,size = iterations-unused.sample)
    chain      = m.chain$chain      [samp.seq,]
    model.chain= m.chain$model.chain[samp.seq ]
  }else{
    chain      = m.chain$chain      [-(1:unused.sample),]
    model.chain= m.chain$model.chain[-(1:unused.sample) ]
  }

  acceptance = apply(chain,2, function(x){1-mean(duplicated(x[x!=0]))}) * 100

  # Stop the clock
  cat('\n Procedure finished in ',( proc.time() - ptm  )[1],' seconds. \n')

  out = list(chain           =  chain              ,
             formula         = formula                    ,
             acceptance.rate = acceptance                 ,
             data = data        , iterations = iterations ,
             q.par = q.par      , b.par = b.par           ,
             initials = initials, burn.in = burn.in       ,
             dist.q =  dist.q   , dist.b =  dist.b        ,
             para.b = para.b    , para.q = para.q         ,
             dist.l =  dist.l   , l.par  = l.par          ,
             fixed.l = fixed.l                            ,
             v.scale = v.scale                            ,
             sampling = sampling                          ,
             RejAccChain = m.chain$RejAcc                 ,
             cov.m = cov.m                                ,
             jeffry = jeffry                              ,
             error    = m.chain$error                     ,
             x  =m.chain$x , y = m.chain$y                ,
             minf = m.chain$minf , minState = m.chain$minState,
             lq.extra = lq.extra               ,
             lbeta.extra = lbeta.extra         ,
             penalized = penalized             ,
             prior.function = prior.function   ,
             lb = m.chain$lb , lq = m.chain$lq ,
             model.chain = model.chain         ,
             logit = logit, RJ = m.chain$RJ    ,
             duration =  proc.time() - ptm
  )
  return(out)
}

plot.bdw = function(x,est = Mode, prob = 0.95,
                    adj = 2, r.outliers = TRUE,
                    density = FALSE, exc.tun = FALSE ,
                    ...){
  dw.object = x
  bdw.plot( dw.object$chain,
            estimate.statistic = est , adj = adj ,
            remove.outliers = r.outliers, sampling = TRUE ,
            ...)
  cat('\n =======',round(prob*100),'% Confidence interval ======= \n')
  dw.HPDinterval   ( dw.object,sub='Parameters',prob = prob,
                     bw = adj, remove.outliers = r.outliers,
                     density = density,
                     exclude.tuning = exc.tun,...)
}


# Draw a plot for bayesian object
bdw.plot <- function(dw.object,
                     estimate.statistic = Mode,
                     adj = 1 , truepar = NA, remove.outliers = TRUE,
                     sampling = TRUE  ,lw = 3, PlotAdj=c(0,0),
                     lableSize=2, ...)
{
  if(requireNamespace('coda') == FALSE){
    stop ('Package Coda is needed in order to execute this function!' ,call. = FALSE)
  }
  truepar= as.matrix(truepar)
  chain  = dw.object$chain
  if(dw.object$fixed.l > 0) {chain = chain[,-ncol(chain)]}
  #### initials ####
  acc    = dw.object$acceptance.rate; acc=round(acc,2)
  ncc    = ncol(chain)
  nrc    = nrow(chain)
  nctp   = ncol(truepar)
  col.names = colnames(truepar,do.NULL = FALSE)
  names     = chain.name(ncc,dw.object$lq,dw.object$lb)
  names2    = extract.vars(dw.object,reg=TRUE)
  symb.p    = names$symb.p
  symb.p2   = names2$symb.p
  scl    = dw.object$v.scale;  scl=round(scl,3)
  lw     = lw
  col    = 3
  lty    = 3
  lsize  = lableSize
  #int.shift = 1.1

  if(all(col.names!="")) {cname.chi =col.names}else{cname.chi='E/T'}
  if(length(scl) == 1) scl = rep(scl,ncc)
  if(exists('n.repeat',where = dw.object) == TRUE) {nrep=dw.object$n.repeat}else{nrep=1}

  par(mfrow = c(2,3),cex=.75)
  for(i in 1:ncc){
    chi = chain[,i]
    chi = na.omit(chi)
    if((sampling == TRUE) && (nrep>1) ){
      schi = sample(chi,nrc/nrep)
      pchi = chi[seq(1,nrc,length.out = nrc/nrep )]
      rnd.rep = round(runif(1,0,nrep-1)) /nrep * nrc
      achi =  (rnd.rep+1):(rnd.rep + nrc/nrep)
    }else{
      schi = pchi = achi = chi
    }
    if(dw.object$RJ){
      schi = schi[schi!=0];if(length(schi)==0){schi=c(0,0)}
      pchi = pchi[pchi!=0];if(length(pchi)==0){pchi=0}
      achi = achi[achi!=0];if(length(achi)==0){achi=0}
      chi  =chi  [chi !=0];if(length(chi) ==0){chi =0}
    }
    m = estimate.statistic(schi);#m=round(m,2)
    if(remove.outliers == TRUE) {
      schi = remove_outliers(schi)
      schi = na.omit(schi)
    }
    #schi = c(min(schi)*4/5,schi,max(schi)*6/5)
    #pchi = c(min(pchi)*4/5,pchi,max(pchi)*6/5)
    ac=acc[i]
    hist(schi, add=0, probability = TRUE,
         #breaks = 30,
         main = symb.p[i],
         xlab=paste('Bayesian estimation = ',round(m,2), ', AcR = ', round(ac,5)),
         xlim=c(min(schi,na.rm = TRUE)-PlotAdj[1],max(schi,na.rm = TRUE)+PlotAdj[2])         ,
         ...
    )
    lines(density(schi,adjust  = adj),col='black')
    if(nrow(truepar) >= i && all(is.na(truepar) == FALSE) ){
      abline(v=truepar[i,],
             col=col:(col+nctp-1),
             lty=lty:(lty+nctp-1),
             lw=lw
      )
      legend('topleft',
             legend = c(
               paste(cname.chi,':',round(truepar[i,],2)),
               paste('Bayesian:',round(m,2))
             ),
             seg.len = lsize,
             #fill = c(col:(col+nctp-1),2),
             lty  = c(lty:(lty+nctp-1),2),
             col  = c(col:(col+nctp-1),2),
             bg   = 'transparent' ,
             lwd=lw-1)
    }
    abline(v = m , col = 2 , lw=lw , lty = 2)
    #abline(...)
    if (nrc > 50000){ ptype = 'b'} else {ptype='l'}
    plot(pchi,xlab='Index',
         #ylim = c(min(pchi)[1]*int.shift,max(pchi)[1]*int.shift),
         ylab='Estimation', type = ptype ,
         #main=paste('Convergence of parameter ',i ,'/', ncc,sep = ''),
         main = symb.p[i],
         col='gray' ,
         pch ='.' ,
         ylim=c(min(pchi,na.rm = TRUE)-PlotAdj[1],max(pchi,na.rm = TRUE)+PlotAdj[2])
    )
    if(nrow(truepar) >= i && all(is.na(truepar[i,]) == FALSE) ){
      abline(h=truepar[i,],
             col=col:(col+nctp-1),
             lty=lty:(lty+nctp-1),
             lw=lw)
      legend('topleft',
             legend = c(
               paste(cname.chi,':',round(truepar[i,],2)),
               paste('Bayesian:',round(m,2))
             ),
             seg.len = lsize,
             #fill = c(col:(col+nctp-1),2),
             lty  = c(lty:(lty+nctp-1),2),
             col  = c(col:(col+nctp-1),2),
             bg   = 'transparent' , lwd=lw-1)
    }
    abline(h = m , col=2 , lw=lw,lty=2)
    #ACF must be computed from the ORIGINAL chain!
    acf=acf(achi,plot = FALSE , lag.max = ceiling(sqrt(nrc/nrep)))
    if(any(is.nan(acf$acf)== TRUE) ) acf=rep(0,ceiling(sqrt(nrc/nrep)))
    #plot(acf,type='h',xlab='Lag',ylab='ACF',main=paste('ACF of parameter ',i ,'/', ncc,sep = ''),col='gray')
    plot(acf,type='h',xlab='Lag',ylab='ACF',main = symb.p[i],col='gray')
    cat('\r',i,' of ',ncc,' plot completed.')
  }
  par(mfrow=c(1,1))
  if(dw.object$RJ){
    model.chain = dw.object$model.chain
    h=table(model.chain)
    h2=length(h)/2
    p=(h[(h2+1):(2*h2)]/(h[h2:1]+h[(h2+1):(2*h2)]))
    bplot=  barplot(-h[h2:1],col=2,ylim=c(-1.5*max(abs(h)),2.5*(max(abs(h)))),xpd=1, xaxt='n',density = 80)
    barplot(h[(h2+1):(2*h2)],add=1,col=3,xpd=1, xaxt='n')
    #axis(1, at=bplot ,labels=paste('V.',1:(h2),' | ',round(p,2),sep=''),las=3,cex.axis=.8)
    abline(v=bplot,col='gray',lty=3)
    axis(1, at=bplot ,labels=symb.p2,las=3,cex.axis=.8)
    legend('topright', horiz = 1, legend = c('Accepted','Rejected'),col = c(3,2),fill = c(3,2),density = c(100,80))
    text (x=bplot,y=1.5*max(abs(h)),labels = paste('+',round(p*100,1),'%'))
    title(main = 'Reversible-Jumps MH model selection')
  }

  #crosscorr.plot(as.mcmc(chain),main='Chain correlation diagram')
}


#############################
bay.laplace.calc <- function (logpost, mode, ...)
{
  options(warn = -1)
  fit = optim(mode, logpost, gr = NULL, ..., hessian = TRUE,
              control = list(fnscale = -1) )
  options(warn = 0)
  mode = fit$par
  h = -solve(fit$hessian)
  p = length(mode)
  int = p/2 * log(2 * pi) + 0.5 * log(det(h)) + logpost(mode,
                                                        ...)
  stuff = list(mode = mode, var = h, int = int, converge = fit$convergence ==
                 0)
  return(stuff)
}
########BF
dw.BF <- function (dw.object,  est.stat = Mode , ...)
{
  options(warn = -1)
  if (dw.object$chain$RJ){
    mode.v = apply (dw.object$chain$chain , 2 , function(x) { x=x[x!=0] ; est.stat(x)} )
  }else{
    mode.v = apply (dw.object$chain$chain , 2 ,  est.stat )
  }
  fit = optim( par = mode.v , fn = posterior.tot.dw ,
               gr = NULL, ..., hessian = TRUE      ,
               control = list(fnscale = -1)        ,
               para.b  = dw.object$chain$para.b    ,
               para.q  = dw.object$chain$para.q    ,
               dist.q  = dw.object$chain$dist.q    ,
               dist.b  = dw.object$chain$dist.b    ,
               q.par   = dw.object$chain$q.par     ,
               b.par   = dw.object$chain$b.par     ,
               dist.l  = dw.object$chain$dist.l    ,
               l.par   = dw.object$chain$l.par     ,
               x       = dw.object$chain$x         ,
               y       = dw.object$chain$y         ,
               logit   = dw.object$chain$logit     ,
               lbeta.extra = dw.object$chain$lbeta.extra,
               lq.extra = dw.object$chain$lq.extra   ,
               penalized = dw.object$chain$penalized ,
               prior.function = dw.object$chain$prior.function,
               fixed.l = dw.object$chain$fixed.l     ,
               iter = dw.object$chain$iter           ,
               jeffry = dw.object$chain$jeffry
  )
  options(warn = 0)
  mode = fit$par
  h = -solve(fit$hessian)
  p = length(mode)
  int = p/2 * log(2 * pi) + 0.5 * log(det(h)) +
    posterior.tot.dw(par = mode,
                     para.b  = dw.object$chain$para.b    ,
                     para.q  = dw.object$chain$para.q    ,
                     dist.q  = dw.object$chain$dist.q    ,
                     dist.b  = dw.object$chain$dist.b    ,
                     q.par   = dw.object$chain$q.par     ,
                     b.par   = dw.object$chain$b.par     ,
                     dist.l  = dw.object$chain$dist.l    ,
                     l.par   = dw.object$chain$l.par     ,
                     x       = dw.object$chain$x         ,
                     y       = dw.object$chain$y         ,
                     logit   = dw.object$chain$logit     ,
                     lbeta.extra = dw.object$chain$lbeta.extra,
                     lq.extra = dw.object$chain$lq.extra ,
                     penalized = dw.object$chain$penalized,
                     fixed.l = dw.object$chain$fixed.l    ,
                     prior.function = dw.object$chain$prior.function,
                     iter = dw.object$chain$iter          ,
                     jeffry = dw.object$chain$jeffry)

  stuff = list(mode = mode,
               var = h,
               bf = int,
               converge = fit$convergence == 0
  )
  return(stuff)
}

##############################################################################
# Estimate Deviance Information Criterion (DIC)
#
# References:
#   Bayesian Data Analysis.
#   Gelman, A., Carlin, J., Stern, H., and Rubin D.
#   Second Edition, 2003
#
#   Bayesian predictive information criterion for the evaluation of
#     hierarchical Bayesian and empirical Bayes models.
#   Ando, T.
#   Biometrika, 2007
#
# Input:
#   x       : matrix of posterior samples
#   lik     : vector of the likelihood of the posterior samples
#   lik.fun : function that calculates the likelihood
#   ...     : other parameters that are passed to 'lik.fun'
#
# Output:
#   list()
#     DIC   : Deviance Information Criterion
#     IC    : Bayesian Predictive Information Criterion
#     pD    : Effective number of parameters (pD = Dbar - Dhat)
#     pV    : Effective number of parameters (pV = var(D)/2)
#     Dbar  : Expected value of the deviance over the posterior
#     Dhat  : Deviance at the mean posterior estimate
##############################################################################
# calc.dic <- function(x,lik,lik.fun,...) {
#   D.bar <- -2*mean(lik)
#   theta.bar <- summary(x)$statistics[,"Mean"]
#   D.hat <- -2*lik.fun(theta.bar,...)
#   pD <- D.bar - D.hat
#   pV <- var(-2*lik)/2
#   list(DIC=pD+D.bar,IC=2*pD+D.bar,pD=pD,pV=pV,Dbar=D.bar,Dhat=D.hat)
# }
dw.info.Criteria = function (dw.object,  est.stat = Mode ,
                             prob = .95, sample.p = TRUE , ...){
  lik   =       l.tot.dw(dw.object , est.stat = NA   , chain.values = TRUE  , sample.p = sample.p)
  D.bar = -2 * mean(lik)
  D.hat = -2 * (l.tot.dw(dw.object , est.stat = mean , chain.values = FALSE , sample.p = sample.p))
  pD    = D.bar - D.hat
  pV    = var(-2*lik)/2
  #############################
  L  = ( l.tot.dw(dw.object , est.stat = est.stat , sample.p = sample.p) )
  k  = dw.object$chain$lb + dw.object$chain$lq
  if (dw.object$chain$penalized || dw.object$chain$RJ){
    zr = dw.HPDinterval(dw.object = dw.object ,prob =  prob ,draw = FALSE , density = FALSE )
    df = k - sum(zr[1:k,4])
  }else{
    df =k
  }

  n  = length(dw.object$chain$y)
  BIC  = -2*L + k *  log(n)
  CAIC = -2*L + k * (log(n)+1)
  AIC  = -2*L + k *2
  AICc = AIC + 2*k*(k+1)/(n-k+1)
  QIC  = -2 * (L - k*log(k))/n
  stuff = list(    AIC  = AIC ,
                   AICc = AICc,
                   BIC  = BIC ,
                   QIC  = QIC ,
                   CAIC = CAIC ,
                   df   = df   ,
                   DIC=pD+D.bar,
                   BPIC=2*pD+D.bar,
                   pD=pD,pV=pV,Dbar=D.bar,Dhat=D.hat,
                   loglik = L
  )
  return(stuff)

}
summary.bdw = function(object,est = Mode ,  prob =.95 , samp = TRUE , ...){
  dw.object = object
  cat('\r Please wait ...')
  BF = dw.BF (dw.object, est.stat = est , ... )
  cr = dw.info.Criteria ( dw.object ,est.stat = est, prob = prob ,sample.p = samp , ...)
  info = dw.object$chain
  cat('\r\n',
      '============================== Sampler ================================ \n' ,
      '    Iterations : ' , info$iterations    , '\t Logit : ' , info$logit==1 , '\t Scale : ' ,round(info$v.scale,5) , '\n' ,
      '    Rev.Jump   : ' , info$RJ==1         , '\t RegQ  : ' , info$para.q==1 ,'\t RegB  : ' ,info$para.b==1 , '\n' ,
      '    Penalized  : ' , info$penalized==1  , '\t Fixed.penalty : ' , info$fixed.l>0, '\n' ,
      '============================ Model Summary ============================ \n'  ,
      '    AIC : ' , cr$AIC , '\t AICc : ' , cr$AICc , '\t BIC     : ' ,  cr$BIC , '\n' ,
      '    QIC : ' , cr$QIC , '\t CAIC : ' , cr$CAIC , '\t LogPPD  : ' ,   BF$bf , '\n' ,
      '    DIC : ' , cr$DIC , '\t PBIC : ' , cr$BPIC , '\t df      : ' ,   cr$df , '\n' ,
      '=======================================================================  \n'
  )
}
######### likelihood
l.tot.dw <- function(dw.object , est.stat = Mode , chain.values = FALSE , sample.p )
{
  para.q   = dw.object$chain$para.q
  para.b   = dw.object$chain$para.b
  l.par    = dw.object$chain$lb + dw.object$chain$lq
  x        = dw.object$chain$x
  y        = dw.object$chain$y
  logit    = dw.object$chain$logit
  lq.extra = dw.object$chain$lq.extra
  lbeta.extra   = dw.object$chain$lbeta.extra


  if(chain.values == FALSE){
    par      = apply(dw.object$chain$chain[,1:l.par] , 2 , est.stat)
    par      = matrix(par,nrow=1)
  }else{
    par = dw.object$chain$chain[,1:l.par]
    ny  = nrow(par)
    smp = sample(1:ny,round(ny*sample.p))
    par = par[smp,]
  }

  if(      (para.b == FALSE) && (para.q == TRUE)){
    loglik = apply(X =  par , MARGIN = 1 , FUN = loglik.q.dw     , x=x , y=y , lq.extra=lq.extra , lbeta.extra=lbeta.extra , logit = logit)
  }else if((para.b == TRUE ) && (para.q == FALSE)){
    loglik = apply(X =  par , MARGIN = 1 , FUN = loglik.beta.dw  , x=x , y=y , lq.extra=lq.extra , lbeta.extra=lbeta.extra , logit = logit)
  }else  if((para.b == TRUE) && (para.q == TRUE)) {
    loglik = apply(X =  par , MARGIN = 1 , FUN = loglik.qbeta.dw , x=x , y=y , lq.extra=lq.extra , lbeta.extra=lbeta.extra , logit = logit)
  }else{
    loglik = apply(X =  par , MARGIN = 1 , FUN = loglik.sim.dw   , x=y       , lq.extra=lq.extra , lbeta.extra=lbeta.extra , logit = logit)
  }
  return(  loglik   )
}
########### HPDinterval
dw.HPDinterval <- function (dw.object , prob = 0.95 , bw = 2.5 ,
                            remove.outliers = TRUE ,
                            draw = TRUE , density = FALSE , exclude.tuning = FALSE,
                            var.lab = NA  ,fixedcol = 1 , ...){
  if(requireNamespace('coda') == FALSE){
    stop ('Package Coda is needed in order to execute this function!' ,call. = FALSE)
  }
  lq = dw.object$chain$lq
  lb = dw.object$chain$lb
  chain = dw.object$chain$chain
  nc.chain = ncol(chain)
  if((nc.chain> lq+lb) && (exclude.tuning==TRUE)){
    chain    = chain[,1:(lb+lq)]
    nc.chain = nc.chain -1
  }

  if (remove.outliers == TRUE){
    chi = apply(chain,2,remove_outliers)
  }else{
    chi = chain
  }
  if(dw.object$chain$RJ){
    est   = apply(chi,2,function(x)   {x=x[x!=0] ; x=na.omit(x); Mode(x = x)})
    res   = apply(chi,2,function(x)   {x=x[x!=0] ; x=na.omit(x); if(length(x)>5) {coda::HPDinterval(coda::as.mcmc(x),prob = prob)}else{c(-10^-8,10^-8)} })
    if(density==TRUE) {d     = apply(chi,2,function(x,dw){x=x[x!=0] ; x=na.omit(x); density(na.omit(x),adjust  = bw)})}
  }else{
    est   = apply(chi,2,function(x)   {x=na.omit(x); Mode(x = x)})
    res   = apply(chi,2,function(x)   {x=na.omit(x); if(length(x)>5) {coda::HPDinterval(coda::as.mcmc(x),prob = prob)}else{c(-10^-8,10^-8)} })
    if(density==TRUE) {d     = apply(chi,2,function(x,dw){x=na.omit(x); density(na.omit(x),adjust =  bw)})}
  }
  #acc.r = round(dw.object$chain$acceptance.rate,2)
  chname = chain.name(nc.chain,lq,lb,var.lab)
  symb   = chname$symb
  symb.p = chname$symb.p
  las    = chname$las
  marb   = chname$marb

  colnames(res)           = symb
  rownames(res)  = c('lower','upper')
  res=t(res)
  lower = res[,1]
  upper = res[,2]
  Zero.included = apply(res,1,function(x){prod(sign(x))<=0})


  if (draw){
    if(lb+lq > 20){cexv= 1/log(nrow(res)+1)+.5}else{cexv=1}
    par(mar =c(marb,4,3,2),cex=cexv)
    nx = nrow(res)
    matplot( cbind ( res , est)                     ,
             col=c(2,3,1) , lty = c ( 3 , 3 , 1 )   ,
             type = c('n' , 'n' , 'p')              ,
             xaxt = 'n'                             ,
             lwd = 1                               ,
             pch = c( 0 , 0 , 1 )                   ,
             ylim = c(min(res,est,na.rm = 1)-.5,max(res,est,na.rm = 1)+.5),
             xlim = c(1, ncol(chain)+1)             ,
             #xlab =  ' \r\n Variables'                   ,
             ylab = paste(round(prob*100,3),'% HPD interval ')
             #xlab = paste('Variables, Acceptance rate = ', round(dw.object$chain$acceptance.rate[1],2),'%' )
    )
    abline ( h = 0 , v =1:nx, col = 'cornsilk2' , lty = 6)
    title(...)
    col.count = 1
    for(k in 1:nx){
      if(abs(res[k,1]-res[k,2]) > 10^-6){
        arrows(k                     ,
               res[k,1]              ,
               k                     ,
               res[k,2]              ,
               length = 0.1          ,
               angle  = 90           ,
               code   = 3            ,
               #col   = c ( col.count%%2+2 , col.count%%2+2 )    ,
               col    =  c(fixedcol*Zero.included[k]+2,fixedcol*Zero.included[k]+2) ,
               lwd    = 2,
               lty    = fixedcol*Zero.included[k]*2+1
        )
        col.count = col.count + 1
      }
    }
    if(density == TRUE){
      for(i in 1:ncol(chain)){
        #with(d,expr = lines(d[[i]]$y/2.65+i,(d[[i]]$x) , lwd =1 ,col=(i)%%2+2 ,lty = i%%2+1))
        with(d,expr = lines(d[[i]]$y/sd(d[[i]]$y)/4+i,(d[[i]]$x) , lwd =1 ,col= fixedcol*Zero.included[i]+2 ,lty = fixedcol*Zero.included[i]*2+1 ) )
        abline(v = i  , lty = 3 , col = 'gray')
      }
    }
    axis(1, 1:nx , labels = symb.p , las = las )
    if(dw.object$chain$fixed.l <= 0  && dw.object$chain$penalized==TRUE){abline (v = ncol(chain)  , lty = 5 , col = 'gray') }
  }

  res=cbind(lower,est,upper,Zero.included)
  return(res)
}


#### Plot density for a fixed x
dw.bay.plot.density <- function (x, dw.object , xmax = 15 ,
                                 smooth = .5 , draw = TRUE , var.lab = NA , ...){

  ncx = ncol(dw.object$chain$x)
  logit = dw.object$chain$logit
  if(     dw.object$chain$para.q == TRUE && dw.object$chain$para.b == FALSE){
    a = 1:ncx
    b = ncx+1
    q.t    = cbind(1 , x) %*% as.matrix(dw.object$res[a])
    b.t    = dw.object$res[b]
    if(logit){
      q      = exp(q.t - log ( 1 + exp(q.t) ) )
    }else{
      q      = exp(-exp(q.t ))
    }
    beta   = rep(b.t , 1)
  }else if(dw.object$chain$para.q == FALSE && dw.object$chain$para.b == TRUE){
    a = 1
    b = 2:(ncx+1)
    q.t    = dw.object$res[a]
    b.t    = cbind( 1 , x ) %*% as.matrix(dw.object$res[b])
    q      = rep( q.t , 1)
    beta   = exp(b.t)
  }else if(dw.object$chain$para.q == TRUE && dw.object$chain$para.b == TRUE){
    a = 1:ncx
    b = (ncx+1):(2*ncx)
    q.t    = cbind( 1 , x ) %*% as.matrix(dw.object$res[a])
    b.t    = cbind( 1 , x ) %*% as.matrix(dw.object$res[b])
    if(logit){
      q      = exp(q.t - log ( 1 + exp(q.t) ) )
    }else{
      q      = exp(-exp(q.t ))
    }
    beta   = exp(b.t)
  }else{
    a = 1
    b = 2
    q     = rep (dw.object$res[a] , 1)
    beta  = rep (dw.object$res[b] , 1)
  }

  res = c()
  d = seq(from = 0 , to = xmax , by = 1 )
  for(i in d ){
    res[i+1] = ddw(i , q = q , beta = beta )
  }
  res = smooth.spline (res , spar = smooth )
  if(draw == TRUE){
    plot( res$x , res$y , type = 'l' ,main = 'PMF' , xlab ='x' , ylab = 'pmf' , ...)
  }
  return(list ( y = res$y , range = res$x , q = q , beta = beta ) )
}




# multicore
bdw.mc <- function(dw.object     ,
                   n.repeat = 10 ,
                   cores = 0   )
{
  if(requireNamespace('foreach')    == FALSE |
     requireNamespace('doParallel') == FALSE |
     requireNamespace('BDWreg')        == FALSE ) stop ('Package needed!')
  #require('foreach')
  # Register all cores
  cat('\r Processing ...')
  if (cores<1) cores = parallel::detectCores()
  doParallel::registerDoParallel(cores = cores)
  ini   = dw.object$chain$initials
  dw.temp.object = dw.object
  ###################
  stime = proc.time()
  i=1
  ch4   = foreach(i = 1:n.repeat , .inorder = FALSE  ) %dopar% {
    dw.temp.object$chain$initials   =  ini/i
    do.call(par.bayesian.tot.dw  ,  dw.temp.object$chain)
  }
  message('Algorithm executed in ',round((proc.time()-stime)[3]) ,' seconds.')


  chain           = ch4[[1]]$chain
  acceptance.rate = ch4[[1]]$acceptance.rate
  minf            = ch4[[1]]$minf
  minState        = ch4[[1]]$minState
  RejAccChain     = ch4[[1]]$RejAccChain
  error           = ch4[[1]]$error
  initials        = ch4[[1]]$initials
  model.chain     = ch4[[1]]$model.chain
  ###########
  foreach::foreach(i = 2:n.repeat,.inorder = TRUE) %do% {
    chain           = rbind(chain        , ch4[[i]]$chain)
    acceptance.rate = (acceptance.rate   + ch4[[i]]$acceptance.rate)
    minf            = c(minf             , ch4[[i]]$minf)
    minState        = rbind(minState     , ch4[[i]]$minState)
    RejAccChain     = c(RejAccChain      , ch4[[i]]$RejAccChain)
    error           = c(error            + ch4[[i]]$error)
    initials        = rbind(initials     , ch4[[i]]$initials)
    model.chain     = c(model.chain      , ch4[[i]]$model.chain)
    cat('\r merging ...',i,' of ' ,n.repeat)
  }
  ### close connection here
  message('Closing connections ...')
  base::closeAllConnections()
  message('End.')

  dw.object$chain$chain=chain
  dw.object$chain$acceptance.rate=acceptance.rate/n.repeat
  dw.object$chain$minf=minf
  dw.object$chain$minState=minState
  dw.object$chain$RejAccChain=RejAccChain
  dw.object$chain$error=error
  dw.object$chain$initials=initials
  dw.object$chain$n.repeat=n.repeat
  dw.object$chain$model.chain=model.chain
  dw.object$chain$iterations=dw.object$chain$iterations*n.repeat
  dw.object$chain$duration = (proc.time()-stime)
  dw.object$chain$cores=cores
  dw.object$all = ch4
  class(dw.object) = 'bdw'
  return(dw.object  )

}

Try the BDWreg package in your browser

Any scripts or data that you put into this service are public.

BDWreg documentation built on May 30, 2017, 12:31 a.m.