# R/bayesian.dw.R In BDWreg: Bayesian Inference for Discrete Weibull Regression

#### Documented in bdw.mcplot.bdwsummary.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
}
#----- 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)),
...
)
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 ='.' ,
)
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)
#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.