# Estimate distribution of aleatoric variable using Bayesian inference

Given a model as follows:

$$y = cx + e$$

where y is the model output, x is the model input, c is an unknown variable and e is a Gaussian model error with zero mean: $$e \sim N(0,\sigma)$$

Data is available for outputs $$y$$ and inputs $$x$$. Suppose now that we want to estimate the parameters $$c$$ and $$\sigma$$, by using a maximum a posteriori algorithm. The algorithm is a simple optimization function (fminunc in Matlab).

If basic Bayesian inference is applied, the likelihood function for the error can be constructed as follows:

$$\mathcal L = \frac{1}{(2\pi)^{N/2}det(\Sigma)^{1/2}} exp\bigg(-\frac{1}{2}x^T\Sigma ^{-1} x\bigg)$$

where $$\Sigma$$ is the covariance matrix and $$x$$ is the model error vector. For simplicity, the priors for both parameters are taken to be uniformly distributed: $$c \sim \mathscr{U}(-\infty,\infty)$$ $$\sigma \sim \mathscr{U}(0 ,\infty)$$

If we now suppose that $$c$$ and $$\sigma$$ are both epistemic variables (unknown but fixed value), of say $$10$$ and $$1$$, respectively, Bayesian inference is rather straightforward and gives us valid estimates for both parameters.

However, if $$c$$ is an aleatory variable with values drawn from $$N(\mu_c,\sigma_c)$$, parameter estimation of $$\mu_c$$, $$\sigma_c$$ and $$\sigma$$ (epistemic) is no longer straightforward. I tried to adjust the likelihood function, so that the computed model output randomly draws a value for $$c$$ from $$N(\mu_c,\sigma_c)$$ for each likelihood function evaluation, but this doesn't seem to work. The algorithm finds a valid estimate for $$\mu_c$$ and $$\sigma$$, but the estimate for $$\sigma_c$$ goes to zero, which is incorrect.

Is there a way to estimate the hyperparameters of the aleatory variable $$c$$, using the likelihood function as specified above?