Estimate parameters of a normal knowing the density function

Question

I have a matrix with 3 columns. The first column contains values of a variable $x_1 \in [-1,1]$. The second column contains values of a variable $x_2 \in [-1,1]$. The third column contains a variable $y$ s.t. $y = p(x_1,x_2)$, where $p$ is the probability density function of a multivariate normal distribution $p = \mathcal{N}(x_1, x_2 | \mu,\Sigma)$, with $\mu \in \mathcal{R}^2$ is the mean of the distribution and $\Sigma \in \mathcal{R}^{2 \times 2}$ is the variance-covariance matrix.

How can I estimate $\mu$ and $\Sigma$?

I thought that I could estimate $\mu$ by taking the mode of the distribution; that is, $\mu = [\bar{x_1}, \bar{x_2}] = \arg_{\text{max}} p(x_1, x_2)$. This is because the mode of a Gaussian coincides with the mean.

But I have no idea for the variance-covariance matrix.

Answer 1

It is essentially the definition of a multivariate normal distribution that the logarithm of its density is a quadratic function of the variables. Thus, all we need do is regress $\log(y)$ on the $x_i, x_i^2,$ and $x_ix_j$ and compute $\mu$ and $\Sigma$ from the parameter estimates.

The quadratic coefficients in the regression are the coefficients of

$$-\frac{1}{2}x^\prime \Sigma^{-1} x.$$

Thus, if we form the matrix of these coefficients, multiply it by -2, and invert that, we will have $\Sigma.$ Since the linear coefficients $\beta$ are those of

$$\mu^\prime \Sigma^{-1}x,$$

we need only solve

$$\Sigma^{-1}\mu = \beta$$

to find $\mu.$

If some values of $y$ are recorded as zeros, you may either ignore them or you could use censored regression methods.

---

To illustrate, here is commented R code showing an implementation, including generation of random data in a data frame X, least-squares fitting of the parameters, and testing the method by comparing its output to the known parameters used in the data-generation process. It is coded to work with $d$-dimensional data where $d\ge 1.$ (Don't try it in more than about $d=50$ dimensions unless you have appreciable computing power!)

#
# Compute the multivariate normal log density.
#
logmvdnorm <- function(x, mu, sigma) {
  d <- length(mu)
  y <- t(x) - mu
  z <- solve(sigma, y)
  (colSums(matrix(z*y, nrow=d)) + log(2*pi)*d + c(determinant(sigma, TRUE)$modulus))/(-2)
}
#
# Specify the data.
#
d <- 2               # Dimensions
n <- d*(d+3)/2 + 1   # Data points
#
# Obtain a random multivariate Normal distribution.
#
Q <- qr.Q(qr(matrix(rnorm(d^2), d)))
rho <- rexp(d, rate=3)
sigma <- crossprod(Q, diag(rho, d, d)%*%Q)
mu <- rnorm(d)
#
# Generate data points.
#
vars <- paste0("x.", 1:d)
X <- as.data.frame(matrix(runif(d*n, -1, 1), ncol=d, dimnames=list(NULL, vars)))
X$y <- logmvdnorm(X, mu, sigma) # These already are the *logarithms* of the density
#------------------------------------------------------------------------------#
#
# Fit a density to the data.
# Begin by constructing the formula for a quadratic fit.
#
d <- length(mu) # Dimension
vars <- paste0("x.", 1:d)
vars2 <- outer(vars, vars, paste, sep="*")
s <- c("1", vars, paste0("I(", vars2[!lower.tri(vars2)], ")"))
f <- as.formula(paste("y ~", paste(s, collapse="+")))
#
# Estimate the parameters using ordinary least squares.
#
fit <- lm(f, subset(X, !is.na(y)))
#
# Convert the parameter estimates into mean and covariance.
#
xform <- function(beta) {
  #
  # Format the coefficients appropriately.
  #
  k <- length(beta)
  d <- round((sqrt(1+8*k)-3)/2)
  Sigma.m1 <- matrix(0, d, d)
  Sigma.m1[!lower.tri(Sigma.m1)] <- -beta[-(1:(d+1))]
  Sigma.m1 <- Sigma.m1 + t(Sigma.m1)
  #
  # Do the calculations.
  #
  Sigma <- solve(Sigma.m1)
  mu <- solve(Sigma.m1, beta[1:d+1])
  return(list(mu=mu, Sigma=Sigma))
}
parameters <- xform(coef(fit))
#
# Check the result.
#
all.equal(parameters$mu, mu)
all.equal(parameters$Sigma, sigma)