Manual maximum likelihood estimation of a mixture with no closed form

Question

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`}

Answer 1

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"}]