From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Innaccuracy in a gammainc and , how best to override?
Date: Thu, 22 Sep 2011 05:16:26 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 28
Message-ID: <j5egba$da$>
References: <> <j5dedo$1rj$> <>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1316668586 426 (22 Sep 2011 05:16:26 GMT)
NNTP-Posting-Date: Thu, 22 Sep 2011 05:16:26 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:743877

DOD <> wrote in message <>...
> ......
> Integrate[Exp(t x ) (-((L x^(-1 - L/ Log[d] ))/ Log[d] )) ,{x,0,1}]
> ......
- - - - - - - - -
  Dennis, the MGF you are attempting to evaluate as a function of t is:

 MGF(t) = int('exp(t*x)*A*x^(A-1)','x',0,1)

and you apparently wish to use positive values of t as high as 50 combined with an 'A' value of .0267 or thereabouts.  You hope to use gammainc(-t,A) for this purpose.

  In my earlier posting I cast doubts as to the appropriateness of such negative values as inputs to matlab's 'gammainc'.  Actually 'gammainc' does accept them but appears to suffer catastrophic rounding errors for values of -t more negative than -30 or so, and this would explain the large "fluctuations" you observed.  Apparently the algorithm used in 'gammainc' does not envision being called upon by such negative values in its first argument along with a parameter as small as this A.

  However, there is a simple numerical method of evaluating MGF(t) and its derivative using matlab's quadrature routines.  Make the substitution

 u = x^A

so as to obtain

 MGF(t) = int('exp(t*u^(1/A))','u',0,1)


 dMGF(t)/dt = int('u^(1/A)*exp(t*u^(1/A))','u',0,1)

Each of these, and therefore their ratio, has no singularities and can readily be evaluated using various of the quadrature routines, 'quad', 'quadl', etc.

Roger Stafford