Description Usage Arguments Details Value Author(s) See Also Examples

The `GIV`

function generates initial values for use with the
`IterativeQuadrature`

, `LaplaceApproximation`

,
`LaplacesDemon`

, `PMC`

, and
`VariationalBayes`

functions.

1 |

`Model` |
This required argument is a model specification
function. For more information, see |

`Data` |
This required argument is a list of data. For more
information, see |

`n` |
This is the number of attempts to generate acceptable initial values. |

`PGF` |
Logical. When |

Initial values are required for optimization or sampling algorithms. A
user may specify initial values, or use the `GIV`

function for
random generation. Initial values determined by the user may fail to
produce a finite posterior in complicated models, and the `GIV`

function is here to help.

`GIV`

has several uses. First, the
`IterativeQuadrature`

, `LaplaceApproximation`

,
`LaplacesDemon`

, and `VariationalBayes`

functions use `GIV`

internally if unacceptable initial values are
discovered. Second, the user may use `GIV`

when developing their
model specification function, `Model`

, to check for potential
problems. Third, the user may prefer to randomly generate acceptable
initial values. Lastly, `GIV`

is recommended when running
multiple or parallel chains with the `LaplacesDemon.hpc`

function (such as for later use with the `Gelman.Diagnostic`

) for
dispersed starting locations. For dispersed starting locations,
`GIV`

should be run once for each parallel chain, and the results
should be stored per row in a matrix of initial values. For more
information, see the `LaplacesDemon.hpc`

documentation for
initial values.

It is strongly recommended that the user specifies a Parameter-Generating Function (PGF), and includes this function in the list of data. Although the PGF may be specified according to the prior distributions (possibly considered as a Prior-Generating Function), it is often specified with a more restricted range. For example, if a user has a model with the following prior distributions

*beta_j ~ N(0,
1000), j=1,…,5*

*sigma ~ HC(25)*

then the PGF, given the prior distributions, is

`PGF <- function(Data) return(c(rnormv(5,0,1000),rhalfcauchy(1,25)))`

However, the user may not want to begin with initial values that could be so far from zero (as determined by the variance of 1000), and may instead prefer

`PGF <- function(Data) return(c(rnormv(5,0,10),rhalfcauchy(1,5)))`

When `PGF=FALSE`

, initial values are attempted to be constrained
to the interval *[-100,100]*. This is done to prevent numeric
overflows with parameters that are exponentiated within the model
specification function. First, `GIV`

passes the upper and lower
bounds of this interval to the model, and any changed parameters are
noted.

At this point, it is hoped that a non-finite posterior is not
found. If found, then the remainder of the process is random and
without the previous bounds. This can be particularly problematic in
the case of, say, initial values that are the elements of a matrix
that must be positive-definite, especially with large matrices. If a
random solution is not found, then `GIV`

will fail.

If the posterior is finite and `PGF=FALSE`

, then initial values
are randomly generated with a normal proposal and a small variance at
the center of the returned range of each parameter. As `GIV`

fails to find acceptable initial values, the algorithm iterates toward
its maximum number of iterations, `n`

. In each iteration, the
variance increases for the proposal.

Initial values are considered acceptable only when the first two
returned components of `Model`

(which are `LP`

and
`Dev`

) are finite, and when initial values do not change through
constraints, as returned in the fifth component of the list:
`parm`

.

If `GIV`

fails to return acceptable initial values, then it is
best to study the model specification function. When the model is
complicated, here is a suggestion. Remove the log-likelihood,
`LL`

, from the equation that calculates the logarithm of the
unnormalized joint posterior density, `LP`

. For example, convert
`LP <- LL + beta.prior`

to `LP <- beta.prior`

. Now, maximize
`LP`

, which is merely the set of prior densities, with any
optimization algorithm. Replace `LL`

, and run the model with
initial values that are in regions of high prior density (preferably
with `PGF=TRUE`

. If this fails, then the model specification
should be studied closely, because a non-finite posterior should
(especially) never be associated with regions of high prior density.

The `GIV`

function returns a vector equal in length to the
number of parameters, and each element is an initial value for the
associated parameter in `Data$parm.names`

. When `GIV`

fails
to find acceptable initial values, each returned element is `NA`

.

Statisticat, LLC. [email protected]

`as.initial.values`

,
`Gelman.Diagnostic`

,
`IterativeQuadrature`

,
`LaplaceApproximation`

,
`LaplacesDemon`

,
`LaplacesDemon.hpc`

,
`PMC`

, and
`VariationalBayes`

.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | ```
library(LaplacesDemon)
############################## Demon Data ###############################
data(demonsnacks)
y <- log(demonsnacks$Calories)
X <- cbind(1, as.matrix(log(demonsnacks[,c(1,4,10)]+1)))
J <- ncol(X)
for (j in 2:J) X[,j] <- CenterScale(X[,j])
######################### Data List Preparation #########################
mon.names <- c("LP","sigma")
parm.names <- as.parm.names(list(beta=rep(0,J), sigma=0))
pos.beta <- grep("beta", parm.names)
pos.sigma <- grep("sigma", parm.names)
PGF <- function(Data) {
beta <- rnorm(Data$J)
sigma <- runif(1)
return(c(beta, sigma))
}
MyData <- list(J=J, PGF=PGF, X=X, mon.names=mon.names,
parm.names=parm.names, pos.beta=pos.beta, pos.sigma=pos.sigma, y=y)
########################## Model Specification ##########################
Model <- function(parm, Data)
{
### Parameters
beta <- parm[Data$pos.beta]
sigma <- interval(parm[Data$pos.sigma], 1e-100, Inf)
parm[Data$pos.sigma] <- sigma
### Log-Prior
beta.prior <- sum(dnormv(beta, 0, 1000, log=TRUE))
sigma.prior <- dhalfcauchy(sigma, 25, log=TRUE)
### Log-Likelihood
mu <- tcrossprod(Data$X, t(beta))
LL <- sum(dnorm(Data$y, mu, sigma, log=TRUE))
### Log-Posterior
LP <- LL + beta.prior + sigma.prior
Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
yhat=rnorm(length(mu), mu, sigma), parm=parm)
return(Modelout)
}
######################## Generate Initial Values ########################
Initial.Values <- GIV(Model, MyData, PGF=TRUE)
``` |

LaplacesDemonR/LaplacesDemon documentation built on Aug. 15, 2018, 4:34 a.m.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.