# Probability of data of being generated by multiple models

Given one or more observations and multiple models, can I compute the probability of some observed data of being generated by each model?

More specifically, I estimated $N$ dynamic models, each one representing a different regime of a dynamic system, in the form $\boldsymbol{x}_{k+1} = f(\boldsymbol{x}_k) + \epsilon$, where $\boldsymbol{x}$ is the system's state, $f(\cdot)$ is a matrix and $\epsilon \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ represents the noise.

Is it possible to compute the probability (or likelihood) of one or more observations of being generated by each model, in order to detect the most likely regime?

• I'm not sure that that is worth doing as one can just perform a minimization to obtain a best model. – Carl Aug 30 '18 at 19:35
• Minimization is a solution, but leaving aside the specific application, assuming to have $N$ probability density functions (pdfs), because each model is affected by gaussian noise, and to observe $M$ samples ($y_1, y_2, ..., y_M$) may I know the most likely pdf that could have generated that samples? – Pas Aug 31 '18 at 8:54
• Yes, for example, in Mathematica you could do a survey with the FindDistribution command, and then refine the parameter values with a FindDistributionParameters command. – Carl Sep 3 '18 at 19:32

Residuals

Given the data $\lbrace \mathbf{x_1}, \mathbf{x_2}, ... , \mathbf{x_n} \rbrace$ you can compute for each model the $n-1$ residuals:

$$\mathbf{\epsilon_{i+1}} = f(\mathbf{x_i}) - \mathbf{x_{i+1}}$$

and the likelihood is computed according to $\mathbf{\epsilon} \sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$

Simpler expression to use in least sum of squares

Possibly you have multiple models, or you want to optimize some parameter.

When you set $\mathbf{y_i} = \mathbf{x_{i+1}}$ then your problem can be expressed more simple like:

$$\mathbf{y_i} = f(\mathbf{x_{i}})+\mathbf{\epsilon_i}$$

This you can solve with any of the multitude of least sum of squares techniques.

Sidenote

Your model does assume that the randomness in the observations is not in a measurement error (or that those errors are negligible) but only in the generation process of the next $\mathbf{x_{i+1}}$ from the past $\mathbf{x_{i}}$.

If your process is related to true variables $\mathbf{x_{i+1}}$ and $\mathbf{x_{i}}$

$$\mathbf{x_{i+1}} = f(\mathbf{x_{i}}) + \mathbf{\epsilon_{i+1}}$$

but your observed variables $\mathbf{y_{i}}$ have errors in relation to those true (underlying) variables $\mathbf{y_{i}} = \mathbf{x_{i}} + \mathbf{\epsilon^\prime_{i}}$ then your model to fit (which relates the $\mathbf{y_{i}}$ and not the $\mathbf{x_{i}}$) should be

$$\mathbf{y_{i+1}} + \mathbf{\epsilon^\prime_{i+1}} = f(\mathbf{y_{i}}+\mathbf{\epsilon^\prime_{i}}) + \mathbf{\epsilon_{i+1}}$$

since you say that $f(\cdot)$ is a matrix you can solve this (it would be more difficult when $f(\cdot)$ is a non-linear function), and you get a model like before $\mathbf{y_i} = f(\mathbf{x_{i}})+\mathbf{\epsilon_i}$ but the succeeding errors $\mathbf{\epsilon_{i}}$ and $\mathbf{\epsilon_{i+1}}$ are correlated

• Thank you for your detailed reply. It is interesting to also consider a measurement noise. However, I have some doubts about the Residuals section: since each model has already been estimated what I know are $f_{j}(\cdot)$ and $\epsilon_j \sim (\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j), j = 1, ..., N$ where $N$ is the number of models. If I observe $x_1, ..., x_M$, I can compute for each model $y_1, ..., y_{M-1}$ but how do I handle the noise $\epsilon_j$? – Pas Aug 31 '18 at 17:07
• Should I consider $L$ possible outcomes for each model, i.e. $L$ values for $y_1$, $L$ for $y_2$, and so on...?. And then, how could I estimate the likelihood of such data of being generated by each model? In other words, finally I would like something like this: 50% Model 1, 15% Model 2, ...., 0.5% Model N. – Pas Aug 31 '18 at 17:07