How to define and compute the marginal densify function from the joint Gaussian distribution?

I try to numerically compute the marginal densify function from the joint Gaussian probability density function given a specific value in the joint density. I use dmvnorm from mvtnorm package:

library(mvtnorm)

myfun = function(x, y, mu, ht) {
  dmvnorm(c(x, y), mean = mu, sigma = ht)
}

But the integrate shows that x and mean have non-conforming size when I run the codes:

integrate(myfun,
          lower = -10,
          upper = 2,
          y = -2,
          mu = c(0.06, 0.03),
          ht = diag(2))

When I test myfun(-1,-2, mu = c(0.06, 0.03), ht = diag(2)), the function works. But I am not sure why the integrate function does not work.

The way your function is written, it can only handle a single x value at a time. If you pass multiple values of x at once (as integrate does), the single y value is just stuck on to the end of it in your call to dmvnorm rather than being interpreted as a fixed y value.

You can use cbind instead of c to fix this, since this will create a matrix of x values in the first column, with the second column being filled with the fixed y value:

library(mvtnorm)

myfun = function(x, y, mu, ht) {
    dmvnorm(cbind(x, y), mean = mu, sigma = ht)
}

integrate(myfun,
          lower = -10,
          upper = 2,
          y = -2,
          mu = c(0.06, 0.03),
          ht = diag(2))
#> 0.04949283 with absolute error < 7.5e-09