## Description

Evaluates the posterior expectation of a (set of) parameters with a given likelihood function and multivariate normal prior distribution by Gauss-Hermite quadrature.

## Usage

 ```1 2``` ```eval.quad(FUN = function(x) 1, X = NULL, ..., W = NULL, forcePD = TRUE, debug = FALSE) ```

## Arguments

 `FUN` log likelihood function of the parameters to be estimated. Defaults to `function(x) 1`, in which case only the prior likelihood is evaluated. `X` Matrix of quadrature points, see `init.quad`. Alternatively, the list of quadrature points and weights produced by `init.quad`. `...` Additional arguments passed on to FUN. `W` Vector of weights, or `NULL` (the default) if provided by `X`. `forcePD` Logical, should the returned estimate be forced to the nearest positive definite matrix - if not already PD? If TRUE (Default: TRUE), `nearPD` is used to arrive at the closest PD matrix. `debug` Logical, should we return the results of FUN?

## Details

The evaluated function is assumed to have a multivariate normal prior distribution, with a mean vector and covariance matrix specified in `init.quad`. The log-likelihood function defaults to an identity function `FUN(x) = 1`, which reduces the distribution under evaluation to specified prior distribution.

The integral under evaluation is;

Int g(X | mu, Sigma) * f(X) * X dX

where g(X | mu, Sigma) is the prior likelihood of X, and f(X) is the likelihood of X in function `FUN`. If left default, the result is the expectation E(X), where X ~ N(mu, Sigma).

Note: FUN is evaluated in a loop, vectorisation is a future possibility. FUN must return a single scalar on the natural log-scale.

## Value

A vector with the evaluated integrals, with attribute `variance` containing the (co)variance (matrix) of the estimate(s), or the positive definite matrix closest to the estiamted covariance matrix.

`init.quad` for creating quadrature points.

## Examples

 ``` 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 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71``` ```### Basic example; E(X), X ~ N(0,1) grid <- init.quad(Q = 1, prior = list(mu = 0, Sigma = diag(1))) eval.quad(X = grid) ### Example; Rasch model person parameter # E(theta), theta ~ N(0,1) * P(X = 1 | theta, beta), P is simplified rasch model # set up rasch model with fixed beta, returns LL rasch <- function(theta, beta, responses){ p <- exp(theta - beta)/(1 + exp(theta - beta)) q <- 1 - p return(log(p) * sum(responses == 1) + log(q) * sum(responses == 0)) } # when theta == beta, P(X = 1) = .5, generate some bernoulli trials with p = .5 responses <- rbinom(5, 1, .5) # get EAP estimate for theta, prior N(0,1) eval.quad(rasch, grid, beta = 0, responses = responses) # with more data, the estimate becomes more accurate, and variance decreases eval.quad(rasch, grid, beta = 0, responses = rbinom(20, 1, .5)) eval.quad(rasch, grid, beta = 0, responses = rbinom(50, 1, .5)) eval.quad(rasch, grid, beta = 0, responses = rbinom(100, 1, .5)) ### problem; the result starts to 'snap' to the closest quadrature point when # the posterior distribution is too dissimilar to the prior. evals <- eval.quad(rasch, grid, beta = 0, responses = rbinom(100, 1, .5), debug = TRUE) evals.values <- attr(evals, "values") # posterior density after 40 items p <- plot(function(x) exp(dnorm(x, log = TRUE) + rasch(x, beta = 1, responses = rbinom(100, 1, .5))), from = -3, to = 3) # quadrature points used points(grid\$X, exp(grid\$W)*max(p\$y), pch = 20) # the evaluation relies almost completely on one quadrature point, # which causes results to 'snap' to that point. # we could add more quadrature points... grid2 <- init.quad(Q = 1, ip = 20) points(grid2\$X, exp(grid2\$W)*max(p\$y), pch = 20, col = "grey") # but if the posterior is not centered on the prior, this quickly fails: p <- plot(function(x) exp(dnorm(x, log = TRUE) + rasch(x, beta = 2, responses = rbinom(100, 1, .5))), from = -3, to = 3) points(grid2\$X, exp(grid2\$W)*max(p\$y), pch = 20, col = "grey") # additionally, adding extra quadrature points in a multidimensional # problem quickly grows out of control. ### a better solution; adaptive quadrature grid. # say we have an idea of where our parameter is located, through another estimator, # or a previous estimate. # we can then use this to adapt where our quadrature grid should be. # get an estimate; responses <- rbinom(10, 1, .5) est <- eval.quad(rasch, grid, beta = 2, responses = responses) print( est ) # adapt the grid; grid3 <- init.quad(Q = 1, adapt = est) # grid is now much closer to posterior p <- plot(function(x) exp(dnorm(x, log = TRUE) + rasch(x, beta = 2, responses = rep(c(0,1), each = 20))), from = -3, to = 3) points(grid3\$X, exp(grid3\$W)*max(p\$y), pch = 20, col = "grey") est <- eval.quad(rasch, grid3, beta = 2, responses = responses) print(est) ```

### Example output

```Loading required package: mvtnorm
[1] 1.534474e-18
attr(,"variance")
[,1]
[1,]    1
[1] 0.2471169
attr(,"variance")
[,1]
[1,] 0.4852284
[1] 2.324137e-17
attr(,"variance")
[,1]
[1,] 0.3810174
[1] 0.6078907
attr(,"variance")
[,1]
[1,] 0.0107994
[1] -0.6159535
attr(,"variance")
[,1]
[1,] 0.0009283237
[1] 0.6587773
attr(,"variance")
[,1]
[1,] 0.3370857
[1] 0.6965349
attr(,"variance")
[,1]
[1,] 0.3815439
```

MultiGHQuad documentation built on May 2, 2019, 3:42 a.m.