R: Why does the sum of multi-variate normal densities not equal to 1

Question

library(mvtnorm)
sigma <- matrix(c(4,2,2,3), ncol=2)
x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
> sum(dmvnorm(x, mean = c(1,2), sigma = sigma))
[1] 14.07509

I'm simulating 500 draws of a bivariate normal with mean vector 1 2 and variance sigma using rmvnorm.

Then I used the dvmnorm to obtain the densities and summed across all of the draws. However, the sum is > 1. Should the sum be 1 instead? Is there a way to have the densities sum up to 1?

Answer 1

The area under a density is always 1. Therefore the following sums to 1 where we took squares of area .01 around the mean of the distribution, approximating the volume above each square as about the density times the base of the square.

x=seq(-19, 11, by=.1) #a region of 
y=seq(-18,22,by=.1)
s=0
for (i in 1:length(x)) {
  for (j in 1:length(y)) {
    s=s+dmvnorm(c(x[i], y[j]), mean = c(1,2), sigma = sigma)*.01
  }
}
s
0.9999997

There is no clear relationship between random draws from a multivariate normal density, evaluating the pdf at each of those points, and summing to 1.