# Optimizing $\chi^2$ using MCMC

I have measurements of an object.

Let's say I have its length $$L$$, mass $$M$$, and age $$t$$: $$\mathbf y = (10~\text{m},\ 0.01~\text{g},\ 5~\text{s}).$$ I also have the uncertainties on my measurements $$\boldsymbol \sigma = (0.1~\text{m},\ 0.001~\text{g},\ 2~\text{s}).$$ The measurements aren't independent so I actually have a covariance matrix $$\boldsymbol \Sigma$$.

I am numerically simulating these objects using some complex non-linear physical theory, we can call it $$\mathbf f$$.

Given some initial conditions $$\mathbf X = (X_1, X_2, X_3)$$, where $$\mathbf X$$ are parameters of my model, I can generate $$\mathbf f(\mathbf X)=(L, M, t)$$.

Now I want to find the $$\mathbf X$$ (and their uncertaintes) that generated my observed $$\mathbf y$$, a common problem. Specifically, I think I should optimize the $$\chi^2$$:

$$$$\chi^2(\mathbf X) = \sum_i \frac{(y_i - f_i(\mathbf X))^2}{\boldsymbol \sigma^2}$$$$ ...but since the measurements aren't independent, it's actually: $$$$\chi^2(\mathbf X) = \mathbf R'\boldsymbol\Sigma^{-1}\mathbf R \qquad \text{where} \qquad \mathbf R = \mathbf y - \mathbf f(\mathbf X).$$$$

Since $$\chi^2(\mathbf X)$$ may be multi-modal with many solutions, I think I should use MCMC to find the posterior distributions of $$\mathbf X$$.

Thus, I need to minimize the negative log likelihood.

My question is, do I minimize

$$-\log \frac{\chi^2}{2} \qquad \text{or} \qquad -\chi^2/2\text{ ?}$$

Or something else?

I am going to have a go at answering this, but unfortunately it is quite difficult to follow the question due to the excessive amount of notation. Nevertheless, if I understand correctly, this is essentially a problem where you have observed values that are coming from a multivariate normal distribution (which gives the objective function you have specified) and you want to estimate the parameters that determine the expected value of the outcome, based on some complex non-linear function.

Re-specifying the model: First I am going to try re-specifying your problem using less notation. In particular, I am going to refer to the parameter vector as $$\boldsymbol{\alpha}$$ to avoid adding another variable name, and I will assume that the variance is built into the covariance matrix $$\mathbf{\Sigma}$$ so that no additional parameter is needed. This will help to clarify the model. Suppose that you specify your problem as:

$$\mathbf{Y}_i = \mathbf{f}(\boldsymbol{\alpha}) + \boldsymbol{\varepsilon}_i \quad \quad \quad \boldsymbol{\varepsilon}_i \sim \text{N}(\mathbf{0}, \mathbf{\Sigma}).$$

In this model form the values $$\mathbf{y}_1,...,\mathbf{y}_n$$ are your observed vectors (containing measurements of length, mass and age) and the parameters $$\boldsymbol{\alpha}$$ and $$\boldsymbol{\Sigma}$$ are the objects of interest in the inference. In this case you have log-likelihood function given by:

$$\ell_\mathbf{y}(\boldsymbol{\alpha}, \boldsymbol{\Sigma}) = - \frac{n}{2} \det (2 \pi \boldsymbol{\Sigma}) -\frac{1}{2} \sum_{i=1}^n (\mathbf{y}_i - \mathbf{f}(\boldsymbol{\alpha}))^\text{T} \mathbf{\Sigma}^{-1} (\mathbf{y}_i - \mathbf{f}(\boldsymbol{\alpha})).$$

If $$\boldsymbol{\Sigma}$$ is known then this essentially reduces to the form you specified as your objective function, which you have denoted as chi-squared. (This form is where I am getting the idea that you are essentially asking about a multivariate normal model.) Now, it is quite simple to deal with this problem analytically, without using MCMC methods.

Estimating the parameters: Minimising the log-likelihood function (and then adjusting the MLE for the covariance matrix to correct for bias) gives the estimators:

$$\mathbf{f}(\hat{\boldsymbol{\alpha}}) = \frac{1}{n} \sum_{i=1}^n \mathbf{y}_i \quad \quad \quad \hat{\mathbf{\Sigma}} = \frac{1}{n-1} \sum_{i=1}^n \Big( \mathbf{y}_i - \frac{1}{n} \sum_{i=1}^n \mathbf{y}_i \Big) \Big( \mathbf{y}_i - \frac{1}{n} \sum_{i=1}^n \mathbf{y}_i \Big)^\text{T}.$$

(I note that you have asked whether you should minimise the log-likelihood or the original likelihood function. Since the logarithm is a monotonic transformation, these minimisation problems yield the same minimising values.) The estimator for the covariance matrix is the sample covariance, which is a simple closed-form estimator. The estimator $$\hat{\boldsymbol{\alpha}}$$ is potentially a bit more complicated, since it involves inversion of the vector function $$\mathbf{f}$$ (which you have stated is a complex non-linear function). Solving this part depends on how much you know about this function.

• Hi @Ben, thanks for this reply. I agree with your "respecification" of the model, it all seems to line up fine. I am a bit confused however about $f(\hat{\alpha})$. The function/simulation $f$ returns multiple values, but here it seems that it is now returning the mean of the observed values? Can you elaborate on that please? – rhombidodecahedron Nov 18 '18 at 18:14
• I have edited to call this $\mathbf{f}$ instead, to stress that this function is producing a vector (of the same length as $\mathbf{y}_i$). You can also see this from the form of its estimator, since this is a sample mean of $\mathbf{y}_1,...,\mathbf{y}_n$ and is therefore also a vector. – Ben Nov 18 '18 at 22:33
• OK, here is one more point of confusion. I only have one observed vector $\mathbf{y}$ (and an estimate of the covariance matrix), not many. Since you have $1/(n-1)$ where $n$ is the number of observed vectors, that seems to imply division by 0. – rhombidodecahedron Nov 18 '18 at 22:44
• If you only have one observed vector $\mathbf{y}$ then for starters, your specification of the problem should not give it a subscript $i$ (since that suggests that you have many of them with different indexes). In any case, that means you now have a situation where you have only one observation from a normal distribution; you can estimate the mean (which is trivially estimated by the observed value) but you cannot estimate the covariance. – Ben Nov 18 '18 at 22:49
• The subscript was because $y_1 = 10~\text{m}$, $y_2=0.01~\text{g}$, $y_3=5~\text{s}$. – rhombidodecahedron Nov 18 '18 at 22:50

This paper gives an answer. The answer is in order to minimize $$\chi^2$$, one can maximize a log likelihood function $$-\chi^2/2$$.