# Compute moments Compute the moment of any function ff using AGHQ. ```r compute_moment(obj, ...) ## S3 method for class 'list' compute_moment( obj, ff = function(x) 1, gg = NULL, nn = NULL, type = c("raw", "central"), method = c("auto", "reuse", "correct"), ... ) ## S3 method for class 'aghq' compute_moment( obj, ff = function(x) 1, gg = NULL, nn = NULL, type = c("raw", "central"), method = c("auto", "reuse", "correct"), ... ) ## Default S3 method: compute_moment( obj, ff = function(x) 1, gg = NULL, method = c("auto", "reuse", "correct"), ... ) ``` ## Arguments - `obj`: Object of class `aghq` output by `aghq::aghq()`. See `?aghq`. - `...`: Used to pass additional argument `ff`. - `ff`: Any R function which takes in a numeric vector and returns a numeric vector. Exactly one of `ff` or `gg` must be provided. If both are provided, `aghq::compute_moment()` will use `gg`, without warning. - `gg`: The output of, or an object which may be input to `aghq::make_moment_function()`. See documentation of that function. Exactly one of `ff` or `gg` must be provided. If both are provided, `aghq::compute_moment()` will use `gg`, without warning. - `nn`: A numeric scalar. Compute the approximate moment of this order, `E(theta^nn|Y)`. See details. If `nn` is provided, `compute_moment` will use it over `ff` or `gg`, without warning. - `type`: Either `'raw'` (default) or `'central'`, see details. - `method`: Method for computing the quadrature points used to approximate moment. One of `'reuse'` (default) or `'correct'`. See details. The default SHOULD be `'correct'`; it is currently set to `'reuse'` to maintain compatibility of results with previous versions. This will be switched in a future major release. ## Returns A numeric vector containing the moment(s) of ff with respect to the joint distribution being approximated using AGHQ. ## Details If multiple of `nn`, `gg`, and `ff` are provided, then `compute_moment` will use `nn`, `gg`, or `ff`, in that order, without warning. There are several approximations available. The "best" one is obtained by specifying `gg` and using `method = 'correct'`. This recomputes the mode and curvature for the function `g(theta)posterior(theta)`, and takes the ratio of the AGHQ approximation to this function to the AGHQ approximation to the marginal likelihood. This obtains the same relative rate of convergence as the AGHQ approximation to the marginal likelihood. It may take a little extra time, and only works for positive, scalar-valued functions `g`. `method = 'reuse'` re-uses the AGHQ adapted points and weights. It's faster than the correct method, because it does not involve any new optimization, it's just a weighted sum. No convergence theory. Seems to work ok in "practice". "Works" for arbitrary `g`. Specifying `ff` instead of `gg` automatically uses `method = 'reuse'`. This interface is provided for backwards compatibility mostly. However, one advantage is that it allows for vector-valued functions, in which case it just returns the corresponding vector of approximate moments. Also, it only requires the adapted nodes and weights, not the ability to evaluate the log-posterior and its derivatives, although this is unlikely to be a practical concern. Specifying a numeric value `nn` will return the moment `E(theta^nn|Y)`. This automatically does some internal shifting to get the evaluations away from zero, to avoid the inherent problem of multi-modal "posteriors" that occurs when the posterior mode is near zero, and account for the fact that some of the new adapted quadrature points may be negative. So, the actual return value is `E(theta^nn + a|Y) - a` for a cleverly-chosen value `a`. Finally, `type='raw'` computes raw moments `E(g(theta)|Y)`, where `type='central'` computes central moments, `E(g(theta - E(g(theta)|Y))|Y)`. See examples. ## Examples ```r logfteta2d <- function(eta,y) { # eta is now (eta1,eta2) # y is now (y1,y2) n <- length(y) n1 <- ceiling(n/2) n2 <- floor(n/2) y1 <- y[1:n1] y2 <- y[(n1+1):(n1+n2)] eta1 <- eta[1] eta2 <- eta[2] sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 + sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2 } set.seed(84343124) n1 <- 5 n2 <- 5 n <- n1+n2 y1 <- rpois(n1,5) y2 <- rpois(n2,5) objfunc2d <- function(x) logfteta2d(x,c(y1,y2)) funlist2d <- list( fn = objfunc2d, gr = function(x) numDeriv::grad(objfunc2d,x), he = function(x) numDeriv::hessian(objfunc2d,x) ) quad <- aghq(funlist2d,7,c(0,0)) ```

compute_moment function