# Estimate parameters of a normal knowing the density function

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.

• Are the values of $y$ exact or are they measurements that might deviate from the true values at random?
– whuber
May 1, 2019 at 13:38
• They are exact apart from numerical approximation. $y$ is a posterior probability in a Bayesian estimation setting. So y = prior * likelihood. I want to "test" the theoretical result for the parameters of y. By multiplying directly, I should obtain the same parameters that I obtain from the theory May 1, 2019 at 13:40
• Thank you. One more question: if the values of $y$ become essentially zero all around the perimeter of the domain $[-1,1]\times[-1,1],$ there is a particularly simple solution available (namely, compute the first two bivariate moments of the y-weighted x-coordinates), but otherwise a substantially more complicated algorithm is required. Which situation are you in?
– whuber
May 1, 2019 at 13:55
• @whuber it appears to be that there are parts of the perimeter where the density is not near to 0 May 1, 2019 at 15:08
• I know that the covariance matrix is diagonal, if it simplifies things May 1, 2019 at 15:14

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)

• Thank you very much, I will read your explanation carefully May 1, 2019 at 20:53