2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Are the values of $y$ exact or are they measurements that might deviate from the true values at random? $\endgroup$ – whuber May 1 at 13:38
  • 2
    $\begingroup$ 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 $\endgroup$ – raffamaiden May 1 at 13:40
  • 1
    $\begingroup$ 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? $\endgroup$ – whuber May 1 at 13:55
  • 1
    $\begingroup$ @whuber it appears to be that there are parts of the perimeter where the density is not near to 0 $\endgroup$ – raffamaiden May 1 at 15:08
  • $\begingroup$ I know that the covariance matrix is diagonal, if it simplifies things $\endgroup$ – raffamaiden May 1 at 15:14
4
$\begingroup$

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)
$\endgroup$
  • 1
    $\begingroup$ Thank you very much, I will read your explanation carefully $\endgroup$ – raffamaiden May 1 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.