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

Question

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?

Answer 1

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.