#Discussion
This is maybe not the main problem but the frequent calls to LogLikelihood
cause a lot of symbolic computation that tends to be slow. Indeed, we are in a situation in which the outcome of LogLikelihood
has a simple symbolic expression. Thus, we may exploit that to perform the symbolic computation only once by defining:
cuts the computation time in half.
Moreover, by setting
g[λ1_][u_] = Simplify[LogLikelihood[PoissonDistribution[λ1], {u}], u >= 0];
which exploits that g
is never fed a negative second argument, we obtain a vectorized functions g[λ1]
which can be applied to lists without Map
or MapThread
.
#Updated Implementation
With some further minor changes, this is the improved code.
fitBivariatePoissonModel2[x_?ListQ, init3_: 1, maxIteration_: 300,
prec_: 10^-8] :=
Block[{n, s, like, zeroQ, lambda1, lambda2, lambda3, var1, var2,
difllike = 1000., loglike0 = 1000., i, lbp1, lbp2, loglike,
loglikeli, f, g, χ, meanvar1, meanvar2},
i = 0;
lambda1 = 0.;
lambda2 = 0.;
lambda3 = N[init3];
{var1, var2} = Transpose[x];
meanvar1 = Mean[N[var1]];
meanvar2 = Mean[N[var2]];
n = Length[x];
s = ConstantArray[0., n];
like = ConstantArray[0., n];
zeroQ = 1. - Unitize[var1 var2];
lambda1 = Max[0.1, Mean[N[var1]] - lambda3];
lambda2 = Max[0.1, Mean[N[var2]] - lambda3];
Block[{u, v, λ1, λ2, λ3},
f[λ1_, λ2_, λ3_][u_, v_] = LogLikelihood[ MultivariatePoissonDistribution[λ3, {λ1, λ2}], {{u, v}}];
g[λ1_][u_] = Simplify[LogLikelihood[PoissonDistribution[λ1], {u}], u >= 0];
];
loglikeli = Reap[
While[difllike > prec && i <= maxIteration, i++;
lbp1 = MapThread[f[lambda1, lambda2, lambda3], {var1 - 1, var2 - 1}];
lbp2 = MapThread[f[lambda1, lambda2, lambda3], {var1, var2}];
χ = N[Unitize[var1 var2]];
s = χ Exp[Log[lambda3] + lbp1 - lbp2];
like = χ lbp2 + (1. - χ) (g[lambda1][var1] + g[lambda2][var2] - lambda3);
Sow[loglike = Total[like]];
difllike = Abs[1. - loglike/loglike0];
loglike0 = loglike;
lambda3 = Mean[s];
lambda1 = meanvar1 - lambda3;
lambda2 = meanvar2 - lambda3;
]
][[2, 1]];
If[i == maxIteration + 1, Print["Maximum iterations reached"]];
Print[ListLinePlot[loglikeli]];
{lambda1, lambda2, lambda3, loglike0}
]
Timing and accuracy comparison:
SeedRandom[1234];
l1 = RandomInteger[{0, 10}, {100, 2}];
aa = fitBivariatePoissonModel[l1, 0.01, 500, 10^-8]; //
AbsoluteTiming // First
bb = fitBivariatePoissonModel2[l1, 0.01, 500, 10^-8]; //
AbsoluteTiming // First
Max[Abs[aa - bb]]
18.7093
0.863204
2.27374*10^-13
There is certainly more potential for improvement.