How to minimize the weighted sum of squared errors with bivariate Gaussian model in R in order to estimate model's parameter?

Question

I would like to estimate 7 parameters of the following bivariate Gaussian model:

So I have the following 7 parameters to estimate, so that the Gaussian model best fit my data X:

The method I would like to use in order to fit the Gaussian bivariate model to my data X is the following:

Because the model has 7 parameters that need to be estimated and because my dataset has not many observations when the inputs "i" and "pi" are big, this could result in imprecise estimates.

So, for estimation purposes, I would like to create a 10 x 10 grid by dividing the "i" / "pi" domain into equally sized squares with the same boundaries (-2%, -1%, 0%, 1%, 2%, 3%, 4%, 5%, 6%) along both dimentions.

Then I want to minimize the weighted sum of squared errors by using the median values of X, "i", "pi" withing each of the regimes delineated by these boundaries:

The weights are directly proportional to the squared root of the number of observations and inversely proportional to the standard deviation of X, both measured within the confines of each regime. This choice of weights forces optimization to pay more attention to areas on the grid that have more observations and less variability.

I don't know how to program this optimization in R, is there any package that can do it?

Answer 1

If $f(i, \pi)$ is a probability density function it has to integrate to 1, so a=0 otherwise it integrates to infinity, then the equation will match Multivariate normal distribution.

The matrix in your equation is two times the covariance of your data, and $\mu_i$ and $\mu_\pi$ is the mean of your data. Here I generate some random gaussian data with preset parameters

original = matrix(rnorm(200), 2, 100)
mu = matrix(c(3,5), 2, 1)
sigma = matrix(c(4, -3, -3, 10), 2, 2)
X = sigma %*% original + rep(mu, ncol(original))

print('Reference means')
print(mu)
print('Reference covariance')
print(sigma %*% t(sigma))

And recover the average and covariance from the data.

mu_ = rowMeans(X)
mu_
sigma2_ = (X - rep(mu_, ncol(X))) %*% t(X - rep(mu_, ncol(X))) / ncol(X)
sigma2_

print('Empirical means')
print(mu_)
print('Empirical covariance')
print(sigma2_)

Comparing your equation to the standard bivariate normal distribution one gets

a = 0
b = (2*pi)*sqrt(det(sigma2_))
mu_i = mu_[1]
mu_pi = mu_[2]
sigma_i = 2*sqrt(sigma2_[1,1])
sigma_pi = 2*sqrt(sigma2_[2,2])
rho = 2 * sigma2_[1,2] / (sigma_i * sigma_pi)