6
$\begingroup$

I have been trying to integrate and compute moments of LogMultinormalDistributionin Mathematica, but the integrals don't seem to evaluate. Say for the first moment, $\langle x \rangle $

Integrate[
 PDF[LogMultinormalDistribution[
  {Subscript[μ, 1], Subscript[μ, 2]},
  {
   {Subscript[σ, 1]^2, Subscript[σ, 1] Subscript[σ, 2] ρ},
   {Subscript[σ, 1] Subscript[σ, 2] ρ, Subscript[σ, 2]^2}
   }
  ]
 , {x, y}
 ]*x
, {x, 0, Infinity}
, {y, 0, Infinity}
]

I also tried typing out the PDF,

Integrate[(2*Pi*σx*σy Sqrt[1 - ρ^2])^-1*
  Exp[-1/2*(1/(
     1 - ρ^2)) (((Log[x] - μx)/σx)^2 - 
      2*ρ*((Log[x] - μx)/σx)*((
        Log[y] - μy)/σy) + ((
       Log[y] - μy)/σy)^2)]*x, {x, 0, Infinity}, {y, 0, Infinity}]

But no luck.

The integrals don't seem to evaluate and either I end up aborting them or the program seems to crash.

In general, I would like compute $\langle x^n y^m \rangle$ for bivariate log normal distribution in mathematica.

Is there a way to speed up the integrals/evaluate them at all.

I am running Mathematica 11.3 on Mac OS Mojave if that helps.

$\endgroup$
  • $\begingroup$ This may not be the problem, but in general you should avoid using Subscript while defining symbols (variables). Subscript[x, 1] is not a symbol, but a composite expression where Subscript is an operator without built-in meaning. You expect to do $x_1=2$ but you are actually doing Set[Subscript[x, 1], 2] which is to assign a DownValues to the operator Subscript and not an OwnValues to an indexed x as you may intend. Read how to properly define indexed variables here $\endgroup$ – rhermans Jul 19 at 17:50
  • $\begingroup$ BTW, Welcome to Mma.SE. Start by taking the tour now and learning about asking and what's on-topic. Always edit if improvable, show due diligence, give brief context, include minimal working example of code and data in formatted form. By doing all this you help us to help you and likely you will inspire great answers. The site depends on participation, as you receive give back: vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$ – rhermans Jul 19 at 17:52
  • $\begingroup$ Thanks for the link to defining indexed variables. I will keep that in mind. removing the subscripts however, doesn't help either to be honest. $\endgroup$ – jcp Jul 19 at 17:53
  • 1
    $\begingroup$ There is nothing wrong with using subscript notation with integration in Mma, provided you don't mix using $x$ and subscript $x$ (say $x_1$) in the same expression. $\endgroup$ – wolfies Jul 19 at 18:34
  • $\begingroup$ @wolfies, I know you think that, you made your point in the Q&A I link, but many people disagree with that, including me. For a new user I think my advice is beneficial, there is nothing to gain and much to loose by using Subscript for anything other than display. $\endgroup$ – rhermans Jul 19 at 18:49
10
$\begingroup$

Life is easier when one can go straight to the answer without going through the integration.

Given a mean vector and covariance matrix

μ = {μ1, μ2};
Σ = {{σ1^2, σ1 σ2 ρ}, {σ1 σ2 ρ, σ2^2}};

the moment generating function for a bivariate normal is

mgf[t1_, t2_] := Exp[{t1, t2}.μ + {t1, t2}.Σ.{t1, t2}/2]

Because the definition of the moment generating function is $E[e^{t_1 X_1+t_2 X_2}]=E[Y_1^{t_1} Y_2^{t_2}]$, then mgf[t1_, t2_] is the function that gives the expectation of $Y_1^{t_1} Y_2^{t_2}$.

mgf[m, n]
(* E^(m μ1 + n μ2 + 1/2 (m (m σ1^2 + n ρ σ1 σ2) + n (m ρ σ1 σ2 + n σ2^2))) *)

You can see that this easily generalizes to more than two dimensions. Here is a reference with the steps of integration: The Lognormal Random Multivariate

$\endgroup$
  • $\begingroup$ There can be no better answer than this. $\endgroup$ – wolfies Jul 19 at 20:10
  • $\begingroup$ @wolfies. I agree. But I wish I would have immediately thought of it. Good thing Google is around. $\endgroup$ – JimB Jul 19 at 20:38
  • $\begingroup$ Ah this makes a lot of sense. Can't believe I did not think of doing this earlier. $\endgroup$ – jcp Jul 19 at 21:45
6
$\begingroup$
Format[μ[n_]] := Subscript[μ, n];
Format[σ[n_]] := Subscript[σ, n]

dist = LogMultinormalDistribution[{μ[1], μ[2]},
   {{σ[1]^2, ρ σ[1] σ[2]}, {ρ σ[1] σ[2], σ[2]^2}}];

The Moments are built-in for this distribution

Moment[dist, {m, n}]

enter image description here

$\endgroup$
  • $\begingroup$ I guess there is a better answer. $\endgroup$ – JimB Jul 20 at 2:40
  • $\begingroup$ @JimB Some of those formulas are baked into the code. $\endgroup$ – Michael E2 Jul 20 at 4:02
4
$\begingroup$

For Mathematica the integral can be very complicated if all general cases are considered for each parameter. One option to speed and simplify things is to analyze more specific cases using Assuming.

You would need to decide if these assumptions are valid, and define your own if not, but these seem quite reasonable to me.

pdf = Assuming[
 And[x > 0, y > 0, ρ^2 < 1, {σ1, σ2, μ1, μ2} ∈ Reals, σ2 > 0],
 Simplify@PDF[
   LogMultinormalDistribution[
    {μ1, μ2},
    {
     {σ1^2, ρ σ1 σ2},
     {ρ σ1 σ2, σ2^2}
     }
    ]
   , {x, y}]
 ]

enter image description here

Notice I have avoided the use of Subscript and used Simplify and indentation to make the code more readable.

Assuming[
 And[x > 0, y > 0, ρ^2 < 1, {σ1, σ2, μ1, μ2} ∈ Reals, σ2 > 0],
 Integrate[
  pdf
  , {x, 0, Infinity}
  , {y, 0, Infinity}
  ]
 ]

(* 1 *)

So the integral is $1$, as expected for a PDF.

Assuming[
 And[x > 0, y > 0, ρ^2 < 1, {σ1, σ2, μ1, μ2} ∈ Reals, σ2 > 0],
 Integrate[
  pdf x
  , {x, 0, Infinity}
  , {y, 0, Infinity}
  ]
 ]

enter image description here

E^(μ1 + σ1^2/2)

Using AbsoluteTiming the solution took 123 seconds in my computer (Mathematica 12 i7-4770 3.4GHz)

$\endgroup$
  • $\begingroup$ Hi. thanks this definitely helps. Running absolute timing on my computer however yields 165.78 and 203.667 s for the first and second computation respectively. Is there anyway to speed this up further? I need to repeat this computation several times. $\endgroup$ – jcp Jul 19 at 18:34
  • $\begingroup$ @jcp see the answer by JimB $\endgroup$ – rhermans Jul 19 at 18:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.