# How to construct the likelihood in this case where I believe errors not to be Gaussian

I have a black box computer code, which takes in input frequency $\omega$, a vector of other predictors $\mathbf{x}$ and a vector of calibration parameters $\boldsymbol{\theta}$, and predicts a vector output $\boldsymbol{\eta}$ having four components. I also have measurements $\mathbf{y}$ for these outputs, which differ from $\boldsymbol{\eta}$ because of some error:

$\mathbf{y}=\boldsymbol{\eta}(\omega,\mathbf{x},\boldsymbol{\theta})+\boldsymbol{\epsilon}$

I have a data set $D=\{(\omega_i,\mathbf{x}_i),\mathbf{y}_i\}_{i=1}^N$ and I want to calibrate my computer code via Bayesian inference, i.e., I want to compute a posterior density for them, given some prior. I can write down the Bayes' formula:

$p(\boldsymbol{\theta}|D)=\frac{p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})}{\int p(D|\boldsymbol{\theta})p(\boldsymbol{\theta})}$

The problem now is that I don't expect the error to be Gaussian, so I'm having troubles writing out the likelihood. The code computes the rigidity and the damping of an advanced seal for a gas turbine. You can get more information here, here and here, but basically just think of a rotor which spins inside a stator. You fit a particular ring on the stator which has the property of making your rotor less sensitive to vibration instabilities, if some design conditions are met when designing your ring. Now, the reasons why I think the error is not Gaussian (and maybe not even iid!) are two:

1. the system being tested is a dynamic linear system. As explained by linear system theory, you can compute the response of a linear system to frequencies $\omega_1,\dots,\omega_m$ by exciting it at all the frequencies in one go, instead than one at a time, because of the superposition principle. This means that my $N$ data points are actually derived from $n=\frac{N}{m}$ separate experiments, each one measuring the response of the system to all $m$ frequencies. I can assume independence among experiments, but I don't think I can consider the $m$ results of a single experiment as independent.
2. The vectors $F=\{\mathbf{y}(\omega_1),\dots,\mathbf{y}(\omega_m)\}$ which are the output of a given experiment are not directly measured. Each vector has 4 components, two rigidities and two dampings. These are not measured, but the displacement in time of the center of mass of the rotor $x_G(t)$ is measured. Then $X_G(\omega)=\text{fft}(x_G(t))$ is computed, and a system identification procedure (basically, a linear system is solved) is performed for each frequency, which leads to the computation of $\{\mathbf{y}(\omega_1),\dots,\mathbf{y}(\omega_m)\}$. In other words, I can say that

$F=\mathcal{L}(x_G)$

where $\mathcal{L}$ is a linear operator. I can safely assume that the measurement error on $x_G(t)$ is Gaussian (the sensor used for the measurement is well-characterized and accurately calibrated), but then I have no idea what becomes of the error on $\mathbf{y}$ after all these transformations!

My problem is, how do I write the likelihood now? To have an idea of the conditional error distribution, I could have a look ar the residuals between the experiments and my code predictions on $D$. But to run my code, I need the values of the calibration parameters, which is what I want to compute! How can I proceed here?

You can use a mixed model to take into account the iid-problem. You make another predictor, $z_{i}$, $i=1, \dots, n$ denoting the experiment number. $z$ should be the random effects. The variables $\omega$ and $x$ will remain fixed effects.

If you have $y=\mathcal L (x_G)$ where $\mathcal L$ is a linear operator and $X_G$ is normally distributed, then $y$ will also be normally distributed.

If you do my suggested choices the likelihood will be the density of Matrix Normal distribution with

$$n=N\\ p=4\\ M=\begin{pmatrix}\eta(\omega_1, x_1, \theta)\\ \vdots\\ \eta(\omega_N, x_N, \theta)\end{pmatrix}\\ U=Z\Psi Z^T+I_N\\ V=\Sigma$$

where $Z$ is a matrix with

$$Z_{ij}=1(\text{observation i belongs to experiment j})$$

and $\Sigma$ is the covariance matrix of the error and $\Psi$ is the covariance matrix of the random effects. If you assume linearity of $\eta$, that is $\eta(\omega,x,\theta)=X(\omega,x)\cdot \theta$ where $X(\omega,x)$ is a matrix, you should be able to use mixed models software straight out of the box.

• Hi, thanks for the answer. Two notes: 1. with the notation $1(\dots)$ you mean the indicator function, right? I.e., $Z_{ij}$ is 1 if observation $i$ belongs to experiment $j$, otherwise it's 0. 2. Linearity of $\eta$ w.r.t. $\theta$ is absolutely an untenable assumption in my case. However, I think I should be able to perform Bayesian inference using MCMC, even if $\eta$ is nonlinear in $\theta$ , as long as I can use your expression for the likelihood. What do you think? Would this be a sound approach? Feb 10, 2017 at 9:39
• @DeltaIV Yes, I mean the indicator function. You can definitely still fit the model even though $\eta$ is nonlinear, it just requires a technique like MCMC. So your approach sounds good. Feb 10, 2017 at 17:09
• sorry, forgot to accept. One last question: in my specific case, which is the random matrix whose distribution is being modeled as Random Normal? It's the difference $\mathbf{X}=\mathbf{y}_i-\boldsymbol{\eta}(\omega_i,\mathbf{x}_i,\theta)$ for all experimental points ($N$ rows) and for the 4 outputs for each point ($p$ columns), so that $\mathbf{X}$ is $n\times p$, right? The Wiki article you linked to claims that the Random Matrix model is not identifiable...doesn't this prevent the use of MCMC? I've advised against using MCMC with non-identifiable models...don't know if it's good advice Feb 13, 2017 at 13:12
• @DeltaIV In my notation I assume that the $N\times 4$ matrix $Y=\{y_i\}_{i=1,...,N}$ are normally distributed with mean $M$ but you can just say that $Y-M$ is normally distributed with mean $0$. The non-identifiability is a problem when $U$ and $V$ are unknown. If $U$ is completely fixed, the model parameters are identifiable(as long as N>4). Perhaps having $\Psi$ unknown could be a problem, but you can always give it a known value. Feb 13, 2017 at 22:49
• @DeltaIV I suggest $\Psi=\sigma^2 I_m$, and then let $\sigma^2$ be a parameter in the model. Feb 15, 2017 at 23:01