Skip to contents

Calculate the log-potential (log-likelihood plus log-density of prior)

Source: R/glm_utils.R
log_potential_from_betaj.Rd

Calculates the log-potential as a function of a new coordinate of the beta parameter vector. Done like this to use the unexporte

Usage

log_potential_from_betaj(
  new_beta_j,
  j,
  current_beta,
  current_eta,
  Y,
  X,
  family,
  beta_prior,
  linear_predictor_calc = "update",
  ...
)

Arguments

new_beta_j

a numeric new value of the j'th component of the beta parameter vector

j

a numeric with the index of the parameter vector

current_beta

the current value of the beta parameter vector in the sampling procedure

current_eta

the current value of the linear predictor corresponding to the current_beta value

Y

a numeric vector of the response variable in which to evaluate the density

X

the design matrix

family

a description of the error distribution and link function to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function. (See family for details of family functions.)

beta_prior

a distribution object created by a function from the distributional package. Could fx. be distributional::dist_normal(mean = 0, sd = 1).

linear_predictor_calc

a character specifying the method used to calculate the linear predictor in each step of the gibbs algorithm. Default is "update", which uses the CGGibbs procedure as described at the start of this section in a vignette. Other option is "naive", which does the usual

...

arguments passed on to the S3 generic log_density

Value

The value of the log-potential having changed the j'th component of the current_beta to new_beta_j

Examples

# Create a test data
n <- 100
x1 <- rnorm (n)
x2 <- rbinom (n, 1, .5)
b0 <- 1
b1 <- 1.5
b2 <- 2
lin_pred <- b0+b1*x1+b2*x2
known_sigma <- 1

y_norm <- rnorm(n, mean = lin_pred, sd = known_sigma)
model_matrix_norm <- as.matrix(
   data.frame(int = 1, X1 = x1, X2 = x2))
b_prior <- distributional::dist_normal(mean = 0, sd = 1)
b_prior_init <- distributional::generate(
    b_prior,
    ncol(model_matrix_norm)
)[[1]]

eta_init <- model_matrix_norm %*% b_prior_init

j <- 1
new_beta_j <- 4

log_potential_from_betaj(new_beta_j = new_beta_j,
                         j = j, current_beta = b_prior_init,
                         current_eta = eta_init,
                         Y = y_norm,
                         X = model_matrix_norm,
                         family = gaussian,
                         beta_prior = b_prior,
                         sd = known_sigma)
#> [1] -526.7973