# Numerical integration involving Inverse Normal CDF

I'd like to evaluate the following numerical integration using Mathematica

$$\ \int_0^T\int_0^\infty xe^{-0.04 s}g(x,s) dxds\$$

where g(x,s) is a Gaussian copula function with say, marginal exponential distributions with parameter A and parameter B (i.e. X~exp(A) and S~exp(B)). In the following MMA attempt, the g(x,s) is defined by GE[Theta,A,B].

This is what I've tried so far:

GE[Theta_, A_, B_] := CopulaDistribution[{"Binormal", Theta}, {ExponentialDistribution[A],ExponentialDistribution[B]}];

Delta = 0.04; A = 0.10; B = 10; Theta = 0.9; T=5;
GExpExp[x_,s_] := PDF[GE[Theta, A, B], {x, s}]

NIntegrate[Exp[-Delta*s]*x GExpExp, {x, 0, Infinity}, {s, 0, T}]


and the following is the error I obtain:

*NIntegrate::inumr: The integrand 1.1547 E^(-10 s-0.1 x-s [Delta]+2 InverseErfc[2 \ Plus[<<2>>]]^2-1.33333 (0.707107 InverseErfc[<<1>>]-Power[<<2>>] \ InverseErfc[<<1>>])^2) x has evaluated to non-numerical values for all sampling points in the region with boundaries {{0.05,36},{0.5,1}}. >>*

I then try to integrate the expression with respect to only one variable, in this case I chose variable x, so that I'd end up with a univariate expression (and hence I might be able to work on a simpler problem from then on)

  Integrate[Exp[-Delta*s]*GExpExp1[x, s], {x, 0, Infinity}]

(* 11.547 E^(-10.04 s + 1.33333 InverseErfc[4 - 4 E^(-10 s)]^2) *)


It seems that eventually I still have to deal with the Inverse[NormalDistribution] function, which I think is the source of this issue. I then continued with the following:

  T = 5;
Assuming[s > 0, LaplaceTransform[NIntegrate[11.547005383792516 E^(-10.04 s + 1.3333333333333337 InverseErfc[4 - 4 E^(-10 s)]^2), {s, 0, T}], s, T]]

(* NIntegrate::inumr: The integrand 11.547 E^(-10.04 s+1.33333 InverseErfc[4+Times[<<2>>]]^2) has evaluated to non-numerical values for all sampling points in the region with boundaries {{1.,2}}. >>  *)


It could be because the function does not converge, but even after reducing the limit interval from [0,T] to [2,3] (and of course changing T to 5 instead of 2) I still can't obtain a number. Does anyone know the correct syntax that I should use so that I can obtain a numerical answer eventually? Your assistance is much appreciated.

• Certainly, T is not a number, and you do not seem to have assigned a value to it... Jun 7, 2013 at 19:31
• Instead of GExpExp = PDF[GE[\[Theta], A, B, {x, s}] .. use: GExpExp = PDF[GE[Theta, A, B], {x, s}] ... note the different placement of the closing brackets ] ... and then you will need to choose a value for T too, or set up K[T] := NIntegrate[ blah] .... not .... K[T] = NIntegrate[ ...] Jun 7, 2013 at 19:34
• You still haven't addressed the thing @wolfies was pointing out. In any event, K is a built-in reserved symbol, so please do not use it to denote your function. Jun 8, 2013 at 4:17
• Thanks for pointing out the errors guys. I have edited the question by including the value of parameter T as well as have written the correct expression of PDF[GE[Theta,A,B]] like how it appears in my code. Jun 8, 2013 at 4:24
• There does indeed appear to be a problem here now (after everything has been fixed up) ... mma 9 hangs or kernel crashes on my Mac. You can force it to produce an answer by adding an option such as NIntegrate[blah, Method -> "MonteCarlo"] ... but, in any event, it would appear that the integral does not converge. Jun 8, 2013 at 8:58

If I use exact coefficients, I get an exact answer with Integrate after a couple of minutes:

GE[Theta_, A_, B_] :=
`