# R/mukernel.R In NegBinBetaBinreg: Negative Binomial and Beta Binomial Bayesian Regression Models

mukernel <- function(y, x, z, betas.now,betas.old,gammas.ini,bpri,Bpri,model,m,ni) {

if (model=="NB1"){
Etas_1 = x%*%betas.old
Etas_2 = z%*%gammas.ini
Mu = exp(Etas_1)
Alpha = exp(Etas_2)

Ym_1 = Etas_1 + y/Mu-1 #Var de trabajo

Q1=as.vector((Mu+Alpha)/Alpha*1/Mu)
SIGMA1=diag(Q1)

} else if (model=="NB2") {

Etas_1 = x%*%betas.old
Etas_2 = z%*%gammas.ini
Mu = exp(Etas_1)
sigma2 = exp(Etas_2)

Ym_1 = Etas_1 + y/Mu-1 #Var de trabajo

Q1=as.vector(sigma2/Mu^2)
SIGMA1=diag(Q1)

} else {

Etas_1 = x%*%betas.old
Etas_2 = z%*%gammas.ini

Mu=inv.logit(Etas_1)
phi=exp(Etas_2)

if (is.null(m)){
m=length(y)
}

if (is.null(ni)){

nro <- length(y)

ni <- rep(m,nro)

}

Ym_1 = Etas_1 + ((y/ni)-Mu)/(Mu*(1-Mu)) #Var de trabajo

Q1=as.vector((ni*Mu*(1-Mu))^(-1)*(phi+ni)/(phi+1))
SIGMA1=diag(Q1)

}

Bpos <- qr.solve(qr.solve(Bpri,tol = 1e-100)+ t(x)%*%qr.solve(SIGMA1,tol = 1e-100)%*%x,tol = 1e-100)
Bpos <- as.matrix(forceSymmetric(as.matrix(Bpos)))
bpos <- Bpos%*%(qr.solve(Bpri,tol = 1e-100)%*%bpri + t(x)%*%qr.solve(SIGMA1,tol = 1e-100)%*%Ym_1)
dmvnorm(t(betas.now),bpos,Bpos) #These functions provide the density function for the multivariate normal
}

## Try the NegBinBetaBinreg package in your browser

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

NegBinBetaBinreg documentation built on May 2, 2019, 10:52 a.m.