Skip to main content
Commonmark migration
Source Link

#Discussion

Discussion

#Updated Implementation

Updated Implementation

#Discussion

#Updated Implementation

Discussion

Updated Implementation

added 114 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322

There is certainly more potential for improvement. Probably the algorithm is not implemented correctly; this comes to my mind because I observed that the algorithm takes rather many interations to "converge" and that the result is very sensitive to decreasing prec.

There is certainly more potential for improvement.

There is certainly more potential for improvement. Probably the algorithm is not implemented correctly; this comes to my mind because I observed that the algorithm takes rather many interations to "converge" and that the result is very sensitive to decreasing prec.

added 2298 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322

#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.

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.

#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.

added 661 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322
Loading
added 3 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322
Loading
added 340 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322
Loading
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322
Loading