# Manual maximum likelihood estimation of a mixture with no closed form

Let $X$ be a random variable that is from a mixture of two Lognormals (Lognormals(mu1,sigma) and Lognormals(mu2,sigma)) and set $Y=X_1+X_2+...+X_{10}$ a sum of 10 random variables from this mixture. There is no closed form for $Y$, but according to the Fenton approximation, the sum of Lognormals is approximately Lognormal distributed. I cannot use FindDistributionParameters in this case and have to use NMaximize instead.

First let me show you how I did this in an easy case, a simple Mixture of two Lognormals (no sums here):

Generate 100 data points with arbitrary mu and sigma:

dat2 = RandomVariate[
MixtureDistribution[{0.2,
1 - 0.2}, {LogNormalDistribution[-2, 0.3],
LogNormalDistribution[1.1, 0.3]}], 100]


define the function to be maximized:

loglikeli2[x_, p_, μ1_, μ2_, σ_] :=
Total[Log[
PDF[MixtureDistribution[{p,
1 - p}, {LogNormalDistribution[μ1, σ],
LogNormalDistribution[μ2, σ]}], x]]]


this function will be maximized:

toMaximize2[p_, μ1_, μ2_, σ_] :=
loglikeli2[dat2, p, μ1, μ2, σ]


NMaximize works like a charm:

NMaximize[{toMaximize2[
p, μ1, μ2, σ], μ1 ∈ Reals, σ >
0, μ2 ∈ Reals, p >= 0,
1 - p >= 0}, {p, μ1, μ2, σ}]


However in the case of the sum of Lognormals this is getting a bit more difficult.

this is the function that has to be maximized: $f(y|p,\mu_1,\mu_2,\sigma) = \sum_{j=0}^{10} \binom{10}{j} p^j (1-p)^{10-j} g_{(j,10-j)}(y)$

where $g(y)$ is the density sum of a random variable from the mixture of two lognormals.

First, generate the data:

mix[p_, m1_, m2_, s_] :=
MixtureDistribution[{p, 1 - p}, {NormalDistribution[m1, s],
LogNormalDistribution[m2, s]}];

dat3 =
Plus @@ RandomVariate[mix[0.75, 0.5, -1.5, 0.2], {2, 100}];


This function approximates the density of a sum of lognormally distributed paramters (Fenton approximation). These are the formulas you can find on the wiipedia page of the Lognormal distribution.

dSumOfLogN[y_, μ_, σ_] :=
Module[{ExpV, V, s2, mu}, ExpV = Total[Exp[μ + 0.5*σ^2]];
(* Expectation *)
V = Total[Exp[μ + 0.5*σ^2]^2*(Exp[σ^2] - 1)];
(*Varinace *)
s2 = Log[1 + V/(ExpV^2)];
(* the new sigma^2 *)
mu = Log[ExpV] - 0.5*s2;
(* the new mu *)
PDF[LogNormalDistribution[mu, Sqrt[s2]], y]]
(* returns the pdf with the new mu and sigma *)


this function builds 2 vectors of length 10. jvector is a vector of two variables that sum to 10. For example with jvector = {3, 7}, mu = {1.2, 4.5} and sigma = {2, 2}, you get a fullmuvector of 3x 1.2 and 7x 4.5 = {1.2,1.2,1.2,4.5,...,4.5} same for fullsigmavector.

dSumOfTypes[y_, jvector_, μ_, σ_] :=
Module[{fullMuVector, fullSigmaVector},
fullMuVector =
Join @@ (ConstantArray @@@ Transpose[{#2, #1}]) &[jvector, μ];
fullSigmaVector =
Join @@ (ConstantArray @@@ Transpose[{#2, #1}]) &[jvector, σ];
dSumOfLogN[y, fullMuVector, fullSigmaVector]]


this is the sum of mixtures $f(x)$:

dSumOfMixtures[y_, pvector_, μ_, σ_] :=
Module[{mysum, jCombis},
jCombis =
Join @@ (Permutations /@
IntegerPartitions[10, {2}, Range[0, 10]]);
(* jCombis generates a vector {{0,10},{1,9},..,{10,0}} *)
mysum =
Sum[PDF[MultinomialDistribution[10, pvector], jCombis[[i]]]*
dSumOfTypes[y, jCombis[[i]], μ, σ], {i, 1,
Length[jCombis]}]; mysum]


this is the loglikelihood function

loglikeli3[y_, p_, μ_, σ_] := Total[Log[dSumOfMixtures[y, p, μ, σ]]]


and this has to be maximized:

toMaximize3[p_?NumericQ, μ_?NumericQ, σ_?NumericQ] := loglikeli3[dat3, p, μ, σ]


but NMaximize fails due to numerical problems and I have no clue how to interpret them

NMaximize[{toMaximize3[{p, 1 - p}, {m1, m2}, {s, s}],
m1 ∈ Reals, s > 0, m2 ∈ Reals, p >= 0, 1 - p >= 0}, {p, m1, m2, s}]


I'm also looking for advice to improve the code, make it shorter, faster etc.

Also, can I somehow make my own distribution out of this with ProbabilityDistribution ?

this is the error I get:

NMaximize::nnum: "The function value -Hold[toMaximize3[{0.731302,1-
0.731302},{0.676536,0.189867},{1.96476,1.96476}]] is not a number at {p,\
[Mu]1,\   [Mu]2,
\[Sigma]} = {0.7313017730309139,0.6765358037113711,0.18986655930905658,1.9647629127390598}

• For starters, your mix[ ] function is supposed to be a mix of two Lognormals, but you have actually set it up as a mix of a Normal and a Lognormal. Another case of using blackbox functions, and not checking intermediate output. Second, there should be no need to use approximating functions in Mma. Third, your question is far too long and convoluted to get to the essence ... – wolfies Aug 2 '14 at 19:19
• @wolfies thanks, you are right, I pasted the wrong mix function. I'd love to avoid having to approximate the mixture function, but everything I tried using TransformedDistribution failed with LogNormals, it works with Normals because there is a closed form for sums. I'll probably delete this post and make a new shorter one. thanks again for taking a glance at my problem, I spent days on this problems now. – spore234 Aug 2 '14 at 19:51

dist = MixtureDistribution[{0.2,
1 - 0.2}, {LogNormalDistribution[-2, 0.3],
LogNormalDistribution[1.1, 0.3]}];

SeedRandom[1];

dat2 = RandomVariate[dist, 100];

estDist =
EstimatedDistribution[dat2,
MixtureDistribution[{p,
1 - p}, {LogNormalDistribution[μ1, σ],
LogNormalDistribution[μ2, σ]}]]


MixtureDistribution[{0.81, 0.19}, {LogNormalDistribution[1.15042, 0.276632], LogNormalDistribution[-1.88732, 0.276632]}]

dat2e = RandomVariate[estDist, 100];

Plot[
Evaluate[PDF[#, x] & /@ {dist, estDist}],
{x, -1, 7},
Frame -> True,
Axes -> False,
FrameLabel -> {x, PDF},
PlotRange -> All,
PlotLegends -> {"dist", "estDist"}]


SmoothHistogram[
{dat2, dat2e}, Automatic, "PDF",
Frame -> True,
Axes -> False,
FrameLabel -> {x, PDF},
PlotLegends -> {"dat2", "dat2e"}]


• thanks, but this is not a MLE for a sum of mixtures but a simple mixture. I like your graphs though. – spore234 Aug 3 '14 at 7:11