Mathematica Asked on September 4, 2021
I need to maximize a function numerically, thus the speed of this function is very important. Here is the most time consuming part of this function.
pPoisson[lambda_, mu_, x_,
y_] := (E^(-lambda - mu)*lambda^x*mu^y)/(x!*y!)
list[a1_, b1_, a2_, b2_, c_, upper_ : 15] :=
Module[{lambda = Exp[a1 - b2 + c], mu = Exp[a2 - b1], m},
m = Table[pPoisson[lambda, mu, i, j], {j, 0, upper}, {i, 0, upper}];
{Total[UpperTriangularize[m, 1], Infinity], Total[Diagonal[m]],
Total[LowerTriangularize[m, -1], Infinity]}]
Do[list[1.3, 0.6, 0.2, 0.2, 0.17, 15], {i, 1000}]; // AbsoluteTiming
(* {1.06158, Null} *)
The pPoisson
is the product of 2 Poisson Distribution probability mass function. The function list
generates a matrix of pPoisson
and calculates the sum of the matrix upper, diagnal and lower part. Running list
1000 times takes 1.06 s. This is not the speed I would like to have. Therefore, I tried to compile it, but the compiled version seems to be 2x slower.
clist = Compile[{{a1, _Real}, {b1, _Real}, {a2, _Real}, {b2, _Real},
{c, _Real}, {upper, _Integer}},
Module[{lambda = Exp[a1 - b2 + c], mu = Exp[a2 - b1], i, j},
{Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Product[x, {x, 1, i}]*
Product[x, {x, 1, j}]), {i, 1, upper}, {j, 0, i - 1}],
Sum[(Exp[(-lambda - mu)]*lambda^i*
mu^i)/(Product[x, {x, 1, i}]^2), {i, 0, upper}],
Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Product[x, {x, 1, i}]*
Product[x, {x, 1, j}]), {j, 1, upper}, {i, 0, j - 1}]}],
CompilationTarget -> "C",
Parallelization -> True, RuntimeOptions -> "Speed"]
Do[clist[1.3, 0.6, 0.2, 0.2, 0.17, 15], {i, 1000}]; // AbsoluteTiming
(* {1.9038489`, Null} *)
Any idea why and how to accelarate it?
The slowness is due to several instances of MainEvaluate
. I replaced the Product
factorials with Gamma
, which turns out to be compilable.
clist = Compile[{{a1, _Real}, {b1, _Real}, {a2, _Real}, {b2, _Real},
{c, _Real}, {upper, _Integer}},
Module[{
lambda = Exp[a1 - b2 + c],
mu = Exp[a2 - b1],
i, j},
{Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Gamma[1 + i]*
Gamma[1 + j]), {i, 1, upper}, {j, 0, i - 1}]
, Sum[(Exp[(-lambda - mu)]*lambda^i*mu^i)/(Gamma[1 + i]^2), {i, 0,
upper}]
, Sum[(Exp[(-lambda - mu)]*lambda^i*mu^j)/(Gamma[1 + i]*
Gamma[1 + j]), {j, 1, upper}, {i, 0, j - 1}]
}]
, CompilationTarget -> "C"
(*,Parallelization[Rule]True*) (* useless without Listable *)
, RuntimeOptions -> "Speed"]
Do[clist[1.3, 0.6, 0.2, 0.2, 0.17, 15], {i, 1000}]; // RepeatedTiming
(* {0.0067, Null} *)
The OP's clist
took about 2 sec.
Correct answer by Michael E2 on September 4, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP