Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
Certainly, T is not a number, and you do not seem to have assigned a value to it... –  J. M. Jun 7 '13 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[ ...] –  wolfies Jun 7 '13 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. –  J. M. Jun 8 '13 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. –  SNMRamli Jun 8 '13 at 4:24
1  
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. –  wolfies Jun 8 '13 at 8:58

1 Answer 1

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

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

Delta = 4/100; A = 10/100; B = 10; Theta = 90/100; T = 5;
GExpExp[x_, s_] := PDF[GE[Theta, A, B], {x, s}]

Integrate[Exp[-Delta*s]*GExpExp[x, s], {x, 0, Infinity}, {s, 0, T}]
(*
  (250 (-1 + E^(251/5)))/(251 E^(251/5))
*)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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