Set up a distribution in LogLikelihood that depends on data

Question

I want to define a loglikelihood function to perform a fit of some data points to a function, like the following, but using the native LogLikelihood[] function.

$ likelihood=\sum_{i}(Log{(\frac{1}{\sqrt{2*\pi*}*\sigma_i}\cdot e^{-\frac{(y_i-f(x_i))^2}{2\sigma_i^2}})}) $

where the function f depends also on the parameters I want to estimate by maximization.

If the $\sigma_i$ are constant I can create the loglikelihood in this way (e.g. using $\sigma=0.2$):

likelihood = LogLikelihood[NormalDistribution[0, 0.2], ydata - f[xdata]];

But if the $\sigma_i$ of the normal distribution are given by the experimental data, they can be different for each data point and I don't if it is possible to set a different distribution parameter for each data point

Obviously, creating the LogL "by hand" with the explicit PDF expression works perfectly, but I saw that if the LogL function is defined using the native function of Mathematica, the computation time is much lower, so I'd like to know if there is a way to solve this issue.

Thanks

Answer 1

Using the NonlinearModelFit example with weights (and I'm assuming that the $\sigma_i$ values are known), one can use the LogLikelihood function in the following manner:

data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}};
σ = data[[All, 1]] + 1;

logL = LogLikelihood[NormalDistribution[0, 1], 
  (data[[All, 2]] - Log[a + b data[[All, 1]]^2])/σ];

In other words you use a Normal[0, 1] distribution and divide the difference between the observed value and the model value by the standard deviation.

Maximizing the log of the likelihood results in

FindMaximum[{logL1, a > 0 && b > 0}, {{a, 2.3}, {b, 0.4}}]
(* {-3.84479, {a -> 2.23705, b -> 0.435855}} *)

This gets you the same estimates of the parameters as using NonlinearModelFit with weights:

nlm = NonlinearModelFit[data, Log[a + b x^2], {a, b}, x, Weights -> 1/σ^2];
nlm["BestFitParameters"]
(* {a -> 2.23705, b -> 0.435855} *)

Something more general

Suppose your data consists of {x, σ, y}:

data = {{0, 1, 1}, {1, 2, 0}, {3, 4, 2}, {5, 6, 4}};

You can use Inner to essentially "vectorize" the LogLikelihood function:

logL = Inner[LogLikelihood[#, {#2}] &, 
  NormalDistribution[a + b #[[1]]^2, #[[2]]] & /@ data, 
  {#[[3]]} & /@ data, Plus][[1]]
(* -(1/2) (-1 + a)^2 - 1/8 (a + b)^2 - 1/32 (-2 + a + 9 b)^2 - 
   1/72 (-4 + a + 25 b)^2 - Log[2] - Log[4] - Log[6] + 2 (-Log[2] - Log[π]) *)

And there are probably cleaner ways to use Inner. But a brute-force sum of log likelihoods seems just about as fast and more straightforward:

logL = Sum[LogLikelihood[NormalDistribution[a + b data[[i, 1]]^2, data[[i, 2]]],
  {data[[i, 3]]}], {i, Length[data]}];]