| Title: | Bayesian Distribution Regression |
|---|---|
| Description: | Implements Bayesian Distribution Regression methods. This package contains functions for three estimators (non-asymptotic, semi-asymptotic and asymptotic) and related routines for Bayesian Distribution Regression in Huang and Tsyawo (2018) <doi:10.2139/ssrn.3048658> which is also the recommended reference to cite for this package. The functions can be grouped into three (3) categories. The first computes the logit likelihood function and posterior densities under uniform and normal priors. The second contains Independence and Random Walk Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithms as functions and the third category of functions are useful for semi-asymptotic and asymptotic Bayesian distribution regression inference. |
| Authors: | Emmanuel Tsyawo [aut, cre], Weige Huang [aut] |
| Maintainer: | Emmanuel Tsyawo <[email protected]> |
| License: | GPL-2 |
| Version: | 0.1.0 |
| Built: | 2024-10-29 06:44:19 UTC |
| Source: | CRAN |
asymcnfB obtains asymmetric bayesian distribution confidence bands
asymcnfB(DF, DFmat, alpha = 0.05, scale = FALSE)asymcnfB(DF, DFmat, alpha = 0.05, scale = FALSE)
DF |
the target distribution/quantile function as a vector |
DFmat |
the matrix of draws of the distribution, rows correspond to
elements in |
alpha |
level such that |
scale |
logical for scaling using the inter-quartile range |
cstar - a constant to add and subtract from DF to create confidence bands if no scaling=FALSE else a vector of length DF.
set.seed(14); m=matrix(rbeta(500,1,4),nrow = 5) + 1:5 DF = apply(m,1,mean); plot(1:5,DF,type="l",ylim = c(min(m),max(m)), xlab = "Index") asyCB<- asymcnfB(DF,DFmat = m) lines(1:5,DF-asyCB$cmin,lty=2); lines(1:5,DF+asyCB$cmax,lty=2)set.seed(14); m=matrix(rbeta(500,1,4),nrow = 5) + 1:5 DF = apply(m,1,mean); plot(1:5,DF,type="l",ylim = c(min(m),max(m)), xlab = "Index") asyCB<- asymcnfB(DF,DFmat = m) lines(1:5,DF-asyCB$cmin,lty=2); lines(1:5,DF+asyCB$cmax,lty=2)
distreg draws randomly from the density of F(yo) at a threshold value yo
distreg(thresh, data0, MH = "IndepMH", ...)distreg(thresh, data0, MH = "IndepMH", ...)
thresh |
threshold value that is used to binarise the continuous outcome variable |
data0 |
original data set with the first column being the continuous outcome variable |
MH |
metropolis-hastings algorithm to use; default:"IndepMH", alternative "RWMH" |
... |
any additional inputs to pass to the MH algorithm |
fitob a vector of fitted values corresponding to the distribution at threshold thresh
data0=faithful[,c(2,1)]; qnt<-quantile(data0[,1],0.25) distob<- distreg(qnt,data0,iter = 102, burn = 2); plot(density(distob,.1),main="Kernel density plot")data0=faithful[,c(2,1)]; qnt<-quantile(data0[,1],0.25) distob<- distreg(qnt,data0,iter = 102, burn = 2); plot(density(distob,.1),main="Kernel density plot")
distreg draws randomly from the density of counterfactual of F(yo) at a threshold
value yo
distreg_cfa(thresh, data0, MH = "IndepMH", cft, cfIND, ...)distreg_cfa(thresh, data0, MH = "IndepMH", cft, cfIND, ...)
thresh |
threshold value that is used to binarise the continuous outcome variable |
data0 |
original data set with the first column being the continuous outcome variable |
MH |
metropolis-hastings algorithm to use; default:"IndepMH", alternative "RWMH" |
cft |
column vector of counterfactual treatment |
cfIND |
the column index(indices) of treatment variable(s) to replace with |
... |
any additional inputs to pass to the MH algorithm |
robj a list of a vector of fitted values corresponding to random draws from F(yo), counterfactual F(yo), and the parameters
data0=faithful[,c(2,1)]; qnt<-quantile(data0[,1],0.25) cfIND=2 #Note: the first column is the outcome variable. cft=0.95*data0[,cfIND] # a decrease by 5% dist_cfa<- distreg_cfa(qnt,data0,cft,cfIND,MH="IndepMH",iter = 102, burn = 2) par(mfrow=c(1,2)); plot(density(dist_cfa$counterfactual,.1),main="Original") plot(density(dist_cfa$counterfactual,.1),main="Counterfactual"); par(mfrow=c(1,1))data0=faithful[,c(2,1)]; qnt<-quantile(data0[,1],0.25) cfIND=2 #Note: the first column is the outcome variable. cft=0.95*data0[,cfIND] # a decrease by 5% dist_cfa<- distreg_cfa(qnt,data0,cft,cfIND,MH="IndepMH",iter = 102, burn = 2) par(mfrow=c(1,2)); plot(density(dist_cfa$counterfactual,.1),main="Original") plot(density(dist_cfa$counterfactual,.1),main="Counterfactual"); par(mfrow=c(1,1))
distreg_cfa.sas takes input object from dr_asympar() for counterfactual semi
asymptotic bayesian distribution. This involves taking random draws from the normal
approximation of the posterior at each threshold value.
distreg_cfa.sas(ind, drabj, data, cft, cfIND, vcovfn = "vcov", iter = 100)distreg_cfa.sas(ind, drabj, data, cft, cfIND, vcovfn = "vcov", iter = 100)
ind |
index of object in list |
drabj |
object from dr_asympar() |
data |
dataframe, first column is the outcome |
cft |
column vector of counterfactual treatment |
cfIND |
the column index(indices) of treatment variable(s) to replace with |
vcovfn |
a string denoting the function to extract the variance-covariance. Defaults at "vcov". Other variance-covariance estimators in the sandwich package are usable. |
iter |
number of draws to simulate |
fitob vector of random draws from density of F(yo) using semi-asymptotic BDR
y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) cfIND=2 #Note: the first column is the outcome variable. cft=0.95*data[,cfIND] # a decrease by 5% cfa.sasobj<- distreg_cfa.sas(ind=2,drabj,data,cft,cfIND,vcovfn="vcov") par(mfrow=c(1,2)); plot(density(cfa.sasobj$original,.1),main="Original") plot(density(cfa.sasobj$counterfactual,.1),main="Counterfactual"); par(mfrow=c(1,1))y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) cfIND=2 #Note: the first column is the outcome variable. cft=0.95*data[,cfIND] # a decrease by 5% cfa.sasobj<- distreg_cfa.sas(ind=2,drabj,data,cft,cfIND,vcovfn="vcov") par(mfrow=c(1,2)); plot(density(cfa.sasobj$original,.1),main="Original") plot(density(cfa.sasobj$counterfactual,.1),main="Counterfactual"); par(mfrow=c(1,1))
distreg.asymp takes input object from dr_asympar() for asymptotic bayesian distribution.
distreg.asymp(ind, drabj, data, vcovfn = "vcov", ...)distreg.asymp(ind, drabj, data, vcovfn = "vcov", ...)
ind |
index of object in list |
drabj |
object from dr_asympar() |
data |
dataframe, first column is the outcome |
vcovfn |
a string denoting the function to extract the variance-covariance. Defaults at "vcov". Other variance-covariance estimators in the sandwich package are usable. |
... |
additional input to pass to |
a mean Fhat and a variance varF
y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) (asymp.obj<- distreg.asymp(ind=2,drabj,data,vcovfn="vcov"))y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) (asymp.obj<- distreg.asymp(ind=2,drabj,data,vcovfn="vcov"))
distreg.sas takes input object from dr_asympar() for semi asymptotic bayesian
distribution. This involves taking random draws from the normal approximation of the
posterior at each threshold value.
distreg.sas(ind, drabj, data, vcovfn = "vcov", iter = 100)distreg.sas(ind, drabj, data, vcovfn = "vcov", iter = 100)
ind |
index of object in list |
drabj |
object from dr_asympar() |
data |
dataframe, first column is the outcome |
vcovfn |
a string denoting the function to extract the variance-covariance. Defaults at "vcov". Other variance-covariance estimators in the sandwich package are usable. |
iter |
number of draws to simulate |
fitob vector of random draws from density of F(yo) using semi-asymptotic BDR
y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) drsas1 = lapply(1:5,distreg.sas,drabj=drabj,data=data,iter=100) drsas2 = lapply(1:5,distreg.sas,drabj=drabj,data=data,vcovfn="vcovHC",iter=100) par(mfrow=c(3,2));invisible(lapply(1:5,function(i)plot(density(drsas1[[i]],.1))));par(mfrow=c(1,1)) par(mfrow=c(3,2));invisible(lapply(1:5,function(i)plot(density(drsas2[[i]],.1))));par(mfrow=c(1,1))y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) drsas1 = lapply(1:5,distreg.sas,drabj=drabj,data=data,iter=100) drsas2 = lapply(1:5,distreg.sas,drabj=drabj,data=data,vcovfn="vcovHC",iter=100) par(mfrow=c(3,2));invisible(lapply(1:5,function(i)plot(density(drsas1[[i]],.1))));par(mfrow=c(1,1)) par(mfrow=c(3,2));invisible(lapply(1:5,function(i)plot(density(drsas2[[i]],.1))));par(mfrow=c(1,1))
dr_asympar computes a normal approximation of the likelihood at a vector of threshold
values
dr_asympar(y, x, thresh, ...)dr_asympar(y, x, thresh, ...)
y |
outcome variable |
x |
matrix of covariates |
thresh |
vector of threshold values on the support of outcome y |
... |
additional arguments to pass to |
a list of glm objects corresponding to thresh
y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus) lapply(drabj,coef); lapply(drabj,vcov) # mean and covariance at respective threshold valuesy = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus) lapply(drabj,coef); lapply(drabj,vcov) # mean and covariance at respective threshold values
fitdist function generates a vector of mean fitted probabilities that constitute the
distribution. This involves marginalising out covariates.
fitdist(Matparam, data)fitdist(Matparam, data)
Matparam |
an M x k matrix of parameter draws, each being a 1 x k vector |
data |
dataframe used to obtain Matparam |
dist fitted (marginalised) distribution
fitlogit obtains a vector of fitted logit probabilities given parameters (pars)
and data
fitlogit(pars, data)fitlogit(pars, data)
pars |
vector of parameters |
data |
data frame. The first column of the data frame ought to be the binary dependent variable |
vec vector of fitted logit probabilities
IndepMH computes random draws of parameters using a specified proposal distribution.
IndepMH(data, propob = NULL, posterior = NULL, iter = 1500, burn = 500, vscale = 1.5, start = NULL, prior = "Uniform", mu = 0, sig = 10)IndepMH(data, propob = NULL, posterior = NULL, iter = 1500, burn = 500, vscale = 1.5, start = NULL, prior = "Uniform", mu = 0, sig = 10)
data |
data required for the posterior distribution |
propob |
a list of mean and variance-covariance of the normal proposal distribution (default:NULL) |
posterior |
the posterior distribution. It is set to null in order to use the logit posterior. The user can specify log posterior as a function of parameters and data (pars,data) |
iter |
number of random draws desired (default: 1500) |
burn |
burn-in period for the MH algorithm (default: 500) |
vscale |
a positive value to scale up or down the variance-covariance matrix in the proposal distribution |
start |
starting values of parameters for the MH algorithm. It is automatically generated but the user can also specify. |
prior |
the prior distribution (default: "Normal", alternative: "Uniform") |
mu |
the mean of the normal prior distribution (default:0) |
sig |
the variance of the normal prior distribution (default:10) |
val a list of matrix of draws pardraws and the acceptance rate
y = indicat(faithful$waiting,70) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x); propob<- lapl_aprx(y,x) IndepMH_n<- IndepMH(data=data,propob,iter = 102, burn = 2) # prior="Normal" IndepMH_u<- IndepMH(data=data,propob,prior="Uniform",iter = 102, burn = 2) # prior="Uniform" par(mfrow=c(3,1));invisible(apply(IndepMH_n$Matpram,2,function(x)plot(density(x)))) invisible(apply(IndepMH_u$Matpram,2,function(x)plot(density(x))));par(mfrow=c(1,1))y = indicat(faithful$waiting,70) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x); propob<- lapl_aprx(y,x) IndepMH_n<- IndepMH(data=data,propob,iter = 102, burn = 2) # prior="Normal" IndepMH_u<- IndepMH(data=data,propob,prior="Uniform",iter = 102, burn = 2) # prior="Uniform" par(mfrow=c(3,1));invisible(apply(IndepMH_n$Matpram,2,function(x)plot(density(x)))) invisible(apply(IndepMH_u$Matpram,2,function(x)plot(density(x))));par(mfrow=c(1,1))
This function creates 0-1 indicators for a given threshold y0 and vector y
indicat(y, y0)indicat(y, y0)
y |
vector y |
y0 |
threshold value y0 |
val
jdpar.asymp takes input object from dr_asympar() for asymptotic bayesian distribution.
It returns objects for joint mutivariate density of parameters across several thresholds.
Check for positive definiteness of the covariance matrix, else exclude thresholds yielding
negative eigen values.
jdpar.asymp(drabj, data, jdF = FALSE, vcovfn = "vcovHC", ...)jdpar.asymp(drabj, data, jdF = FALSE, vcovfn = "vcovHC", ...)
drabj |
object from dr_asympar() |
data |
dataframe, first column is the outcome |
jdF |
logical to return joint density of F(yo) across thresholds in drabj |
vcovfn |
a string denoting the function to extract the variance-covariance. Defaults at "vcov". Other variance-covariance estimators in the sandwich package are usable. |
... |
additional input to pass to |
mean vector Theta and variance-covariance matrix vcovpar of parameters across
thresholds and if jdF=TRUE,
a mean vector mnF and a variance-covariance matrix vcovF of F(yo)
y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) (drjasy = jdpar.asymp(drabj=drabj,data=data,jdF=TRUE))y = faithful$waiting x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) qtaus = quantile(y,c(0.05,0.25,0.5,0.75,0.95)) drabj<- dr_asympar(y=y,x=x,thresh = qtaus); data = data.frame(y,x) (drjasy = jdpar.asymp(drabj=drabj,data=data,jdF=TRUE))
jntCBOM implements calibrated symmetric confidence bands (algorithm 2)
in Montiel Olea and Plagborg-Moller (2018).
jntCBOM(DF, DFmat, alpha = 0.05, eps = 0.001)jntCBOM(DF, DFmat, alpha = 0.05, eps = 0.001)
DF |
the target distribution/quantile function as a vector |
DFmat |
the matrix of draws of the distribution, rows correspond to
indices elements in |
alpha |
level such that |
eps |
steps by which the grid on 1-alpha:alpha/2 is searched. |
CB - confidence band, zeta - the optimal level
set.seed(14); m=matrix(rbeta(500,1,4),nrow = 5) + 1:5 DF = apply(m,1,mean); plot(1:5,DF,type="l",ylim = c(min(m),max(m)), xlab = "Index") jOMCB<- jntCBOM(DF,DFmat = m) lines(1:5,jOMCB$CB[,1],lty=2); lines(1:5,jOMCB$CB[,2],lty=2)set.seed(14); m=matrix(rbeta(500,1,4),nrow = 5) + 1:5 DF = apply(m,1,mean); plot(1:5,DF,type="l",ylim = c(min(m),max(m)), xlab = "Index") jOMCB<- jntCBOM(DF,DFmat = m) lines(1:5,jOMCB$CB[,1],lty=2); lines(1:5,jOMCB$CB[,2],lty=2)
This function generates mode and variance-covariance for a normal proposal distribution for the bayesian logit.
lapl_aprx(y, x, glmobj = FALSE)lapl_aprx(y, x, glmobj = FALSE)
y |
the binary dependent variable y |
x |
the matrix of independent variables. |
glmobj |
logical for returning the logit glm object |
val A list of mode variance-covariance matrix, and scale factor for proposal draws from the multivariate normal distribution.
y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) gg<- lapl_aprx(y,x)y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) gg<- lapl_aprx(y,x)
lapl_aprx2 is a more flexible alternative to lapl_aprx. This creates
glm objects from which joint asymptotic distributions can be computed.
lapl_aprx2(y, x, family = "binomial", ...)lapl_aprx2(y, x, family = "binomial", ...)
y |
the binary dependent variable y |
x |
the matrix of independent variables. |
family |
a parameter to be passed |
... |
additional parameters to be passed to |
val A list of mode variance-covariance matrix, and scale factor for proposal draws from the multivariate normal distribution.
y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) (gg<- lapl_aprx2(y,x)); coef(gg); vcov(gg)y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) (gg<- lapl_aprx2(y,x)); coef(gg); vcov(gg)
logit is the logistic likelihood function given data.
logit(start, data, Log = TRUE)logit(start, data, Log = TRUE)
start |
vector of starting values |
data |
dataframe. The first column should be the dependent variable. |
Log |
a logical input (defaults to |
like returns the likelihood function value.
y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x) logit(rep(0,3),data)y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x) logit(rep(0,3),data)
This is the link function for logit regression
LogitLink(x)LogitLink(x)
x |
Random variable |
val Probability value from the logistic function
par_distreg uses parallel computation to compute bayesian distribution regression for a given
vector of threshold values and a data (with first column being the continuous outcome variable)
par_distreg(thresh, data0, fn = distreg, no_cores = 1, type = "PSOCK", ...)par_distreg(thresh, data0, fn = distreg, no_cores = 1, type = "PSOCK", ...)
thresh |
vector of threshold values. |
data0 |
the original data set with a continous dependent variable in the first column |
fn |
bayesian distribution regression function. the default is distreg provided in the package |
no_cores |
number of cores for parallel computation |
type |
|
... |
any additional input parameters to pass to fn |
mat a G x M matrix of output (G is the length of thresh, M is the number of draws)
data0=faithful[,c(2,1)]; qnts<-quantile(data0[,1],c(0.05,0.25,0.5,0.75,0.95)) out<- par_distreg(qnts,data0,no_cores=1,iter = 102, burn = 2) par(mfrow=c(3,2));invisible(apply(out,1,function(x)plot(density(x,30))));par(mfrow=c(1,1))data0=faithful[,c(2,1)]; qnts<-quantile(data0[,1],c(0.05,0.25,0.5,0.75,0.95)) out<- par_distreg(qnts,data0,no_cores=1,iter = 102, burn = 2) par(mfrow=c(3,2));invisible(apply(out,1,function(x)plot(density(x,30))));par(mfrow=c(1,1))
parLply uses parlapply from the parallel package with
a function as input
parLply(vec, fn, type = "FORK", no_cores = 1, ...)parLply(vec, fn, type = "FORK", no_cores = 1, ...)
vec |
vector of inputs over which to parallel compute |
fn |
the function |
type |
this option is set to "FORK", use "PSOCK" on windows |
no_cores |
the number of cores to use. Defaults at 1 |
... |
extra inputs to |
out parallel computed output
posterior computes the value of the posterior at parameter values pars
posterior(pars, data, Log = TRUE, mu = 0, sig = 25, prior = "Normal")posterior(pars, data, Log = TRUE, mu = 0, sig = 25, prior = "Normal")
pars |
parameter values |
data |
dataframe. The first column must be the binary dependent variable |
Log |
logical to take the log of the posterior.(defaults to TRUE) |
mu |
mean of prior of each parameter value in case the prior is Normal (default: 0) |
sig |
standard deviation of prior of each parameter in case the prior is Normal (default: 25) |
prior |
string input of "Normal" or "Uniform" prior distribution to use |
val value function of the posterior
y = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x) posterior(rep(0,3),data,Log = FALSE,mu=0,sig = 10,prior = "Normal") # no log posterior(rep(0,3),data,Log = TRUE,mu=0,sig = 10,prior = "Normal") # log posterior(rep(0,3),data,Log = TRUE) # use default valuesy = indicat(faithful$waiting,mean(faithful$waiting)) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x) posterior(rep(0,3),data,Log = FALSE,mu=0,sig = 10,prior = "Normal") # no log posterior(rep(0,3),data,Log = TRUE,mu=0,sig = 10,prior = "Normal") # log posterior(rep(0,3),data,Log = TRUE) # use default values
This normal prior distribution is a product of univariate N(mu,sig)
prior_n(pars, mu, sig, Log = FALSE)prior_n(pars, mu, sig, Log = FALSE)
pars |
parameter values |
mu |
mean value of each parameter value |
sig |
standard deviation of each parameter value |
Log |
logical to take the log of prior or not (defaults to FALSE) |
val Product of probability values for each parameter
prior_n(rep(0,6),0,10,Log = TRUE) #log of prior prior_n(rep(0,6),0,10,Log = FALSE) #no logprior_n(rep(0,6),0,10,Log = TRUE) #log of prior prior_n(rep(0,6),0,10,Log = FALSE) #no log
This uniform prior distribution proportional to 1
prior_u(pars)prior_u(pars)
pars |
parameter values |
val value of joint prior =1 for the uniform prior
quant_bdr converts a bayesian distribution regression matrix from par_distreg()
output to a matrix of quantile distribution.
quant_bdr(taus, thresh, mat)quant_bdr(taus, thresh, mat)
taus |
a vector of quantile indices |
thresh |
a vector of threshold values used in a |
mat |
bayesian distribution regression output matrix |
qmat matrix of quantile distribution
RWMH computes random draws of parameters using a specified proposal distribution.
The default is the normal distribution
RWMH(data, propob = NULL, posterior = NULL, iter = 1500, burn = 500, vscale = 1.5, start = NULL, prior = "Normal", mu = 0, sig = 10)RWMH(data, propob = NULL, posterior = NULL, iter = 1500, burn = 500, vscale = 1.5, start = NULL, prior = "Normal", mu = 0, sig = 10)
data |
data required for the posterior distribution. First column is the outcome |
propob |
a list of mean and variance-covariance of the normal proposal distribution (default: NULL i.e. internally generated) |
posterior |
the posterior distribution. It is set to null in order to use the logit posterior. The user can specify log posterior as a function of parameters and data (pars,data) |
iter |
number of random draws desired |
burn |
burn-in period for the Random Walk MH algorithm |
vscale |
a positive value to scale up or down the variance-covariance matrix in the proposal distribution |
start |
starting values of parameters for the MH algorithm. It is automatically generated from the proposal distribution but the user can also specify. |
prior |
the prior distribution (default: "Normal", alternative: "Uniform") |
mu |
the mean of the normal prior distribution (default:0) |
sig |
the variance of the normal prior distribution (default:10) |
val a list of matrix of draws Matpram and the acceptance rate
y = indicat(faithful$waiting,70) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x); propob<- lapl_aprx(y,x) RWMHob_n<- RWMH(data=data,propob,iter = 102, burn = 2) # prior="Normal" RWMHob_u<- RWMH(data=data,propob,prior="Uniform",iter = 102, burn = 2) par(mfrow=c(3,1));invisible(apply(RWMHob_n$Matpram,2,function(x)plot(density(x)))) invisible(apply(RWMHob_u$Matpram,2,function(x)plot(density(x))));par(mfrow=c(1,1))y = indicat(faithful$waiting,70) x = scale(cbind(faithful$eruptions,faithful$eruptions^2)) data = data.frame(y,x); propob<- lapl_aprx(y,x) RWMHob_n<- RWMH(data=data,propob,iter = 102, burn = 2) # prior="Normal" RWMHob_u<- RWMH(data=data,propob,prior="Uniform",iter = 102, burn = 2) par(mfrow=c(3,1));invisible(apply(RWMHob_n$Matpram,2,function(x)plot(density(x)))) invisible(apply(RWMHob_u$Matpram,2,function(x)plot(density(x))));par(mfrow=c(1,1))
simcnfB obtains symmetric bayesian distribution confidence bands
simcnfB(DF, DFmat, alpha = 0.05, scale = FALSE)simcnfB(DF, DFmat, alpha = 0.05, scale = FALSE)
DF |
the target distribution/quantile function as a vector |
DFmat |
the matrix of draws of the distribution, rows correspond to
elements in |
alpha |
level such that |
scale |
logical for scaling using the inter-quartile range |
cstar - a constant to add and subtract from DF to create confidence bands if no scaling=FALSE else a vector of length DF.
set.seed(14); m=matrix(rbeta(500,1,4),nrow = 5) + 1:5 DF = apply(m,1,mean); plot(1:5,DF,type="l",ylim = c(0,max(m)), xlab = "Index") symCB<- simcnfB(DF,DFmat = m) lines(1:5,DF-symCB,lty=2); lines(1:5,DF+symCB,lty=2)set.seed(14); m=matrix(rbeta(500,1,4),nrow = 5) + 1:5 DF = apply(m,1,mean); plot(1:5,DF,type="l",ylim = c(0,max(m)), xlab = "Index") symCB<- simcnfB(DF,DFmat = m) lines(1:5,DF-symCB,lty=2); lines(1:5,DF+symCB,lty=2)