Skip to main content
1 of 2
Miguel
  • 201
  • 1
  • 5

NIntegrate appears to give incorrect results

I am trying to specify a bivariate probability density function in Mathematica. As a check, I would like to confirm that it integrates to one. Here is the function:

f[x1_, x2_, u1_, u2_, v11_, v22_, v12_] := Det[2*Pi*{{v11, v12}, {v12, v22}}]^(-0.5)*(x1*x2*(1 - x1 - x2) )^(-1)*Exp[-0.5*(Log[{x1, x2}/(1 - x1 - x2)] - {u1, u2}).Inverse[{{v11, v12}, {v12, v22}}].(Log[{x1, x2}/(1 - x1 - x2)] - {u1, u2})]

Here, $X_1$ and $X_2$ are the random variables, where $X_1 > 0$, $X_2 > 0$ and $X_1 + X_2 < 1$, and $U_1 \in \mathbb{R}$, $U_2 \in \mathbb{R}$, $V_{11}>0$, $V_{22}>0$ and $V_{12}\in \mathbb{R}$ are parameters.

Using NIntegrateto integrate over the $(X_1,X_2)$ space with $U_1=U_2=V_{12}=0$ and $V_{11}=V_{22}=1$, I get:

NIntegrate[f[x1, x2, 0, 0, 1, 1, 0], {x1, 0, 1}, {x2, 0, 1 - x1}]
(* 0.364031 *)

The answer should be 1 (not 0.364 as above) for any choice of the parameters. I have specified the same function in the R language with only syntax changes. The function in R gives exactly the same values as in Mathematica when supplied with the same arguments (confirming that it is not a programming error). However, I get an integral of 1 in R using adaptive quadrature (i.e. a different answer to Mathematica!).

Any suggestions?

Miguel

Miguel
  • 201
  • 1
  • 5