# Set up a distribution in LogLikelihood that depends on data

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

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jan 16 at 17:38

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 #[]^2, #[]] & /@ data,
{#[]} & /@ data, Plus][]
(* -(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 - Log - Log + 2 (-Log - 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]}];]

• Thank you for your suggestion, I haven't thought to normalize the residual and use a NormalDistribution[0,1]. Jan 17 at 13:26
• But is there a more general way to do it? This could be a non-standard case, but if I want to use for example a Poisson distribution I can't use this method Jan 17 at 13:37
• Editing your question to include that need would get you more targeted help.
– JimB
Jan 17 at 16:40
• @Andrea I added something a bit more general. And, of course, if your function is linear, then GeneralizedLinearModelFit should be considered. And if your error structure consists of an additive normal distribution, then NonlinearModelFit with weights might be the way to go.
– JimB
Jan 17 at 19:59