Description Usage Arguments Details Value See Also Examples

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

1 2 |

`FUN` |
log likelihood function of the parameters to be estimated.
Defaults to |

`X` |
Matrix of quadrature points, see |

`...` |
Additional arguments passed on to FUN. |

`W` |
Vector of weights, or |

`forcePD` |
Logical, should the returned estimate be forced to the nearest positive definite matrix - if not already PD? If TRUE (Default: TRUE), |

`debug` |
Logical, should we return the results of FUN? |

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.

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.

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)
``` |

```
Loading required package: mvtnorm
Loading required package: Matrix
[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
```

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.