Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!15g2000vbb.googlegroups.com!not-for-mail
From: DOD <dcodea@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Innaccuracy in a gammainc and , how best to override?
Date: Thu, 22 Sep 2011 08:03:19 -0700 (PDT)
Organization: http://groups.google.com
Lines: 84
Message-ID: <0963bec6-df81-4dd9-857e-6c255fc489d4@15g2000vbb.googlegroups.com>
References: <e6d2bf35-7011-4fc8-baf4-45138f5e8798@a7g2000yqb.googlegroups.com>
 <j5dedo$1rj$1@newscl01ah.mathworks.com> <7583ba5c-88b2-48a7-8f24-f341f096c335@n40g2000yqb.googlegroups.com>
 <j5egba$da$1@newscl01ah.mathworks.com>
NNTP-Posting-Host: 128.174.250.146
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1316703800 15520 127.0.0.1 (22 Sep 2011 15:03:20 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Thu, 22 Sep 2011 15:03:20 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: 15g2000vbb.googlegroups.com; posting-host=128.174.250.146; posting-account=yN5uLwoAAAD-9ZO1c3IjHfbc8c1Nakp-
User-Agent: G2/1.0
X-Google-Web-Client: true
X-Google-Header-Order: HUALESNKC
X-HTTP-UserAgent: Mozilla/5.0 (Windows NT 5.1; rv:6.0.2) Gecko/20100101 Firefox/6.0.2,gzip(gfe)
Xref: news.mathworks.com comp.soft-sys.matlab:743941

On Sep 22, 12:16 am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> DOD <dco...@gmail.com> wrote in message <7583ba5c-88b2-48a7-8f24-f341f096c...@n40g2000yqb.googlegroups.com>...
> > ......
> > 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)
>
> and
>
>  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

I had previously tried things with quad, quadl and quadgk, but I don't
think it actually helps.  I think you made an error in your
substitution, which renders the method less useful.  The MGF is

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

With the substitution u = x^A, the integral becomes

int (A * exp ( t u ^ (1/A) ) * u ^( (A-1) / A ) ,0,1)

IE, the x term outside the exp() does not disappear.  Even if you use
U = x^(A-1), you still have

int (A * exp (t u^(1/(A-1)) * u  ,0,1)

Unfortunately, quad, quadl and quadgk all have very tough times with
integrals of this form.  I had a question on this very topic some
weeks ago!  For example, let A = .0267, t=10, a number the above
method has no trouble with:

quad(@(x)exp( 10 .* x.^(1./(A-1))) .* x,0,1)

for those values returns NaN, while quadl and quadgk both return
Inf.   So using the ratios of the actual integrals in question is not
a promising avenue.  I think that is the whole reason these sorts of
integrals are implemented as special functions.

I certainly agree that gammainc is perhaps not full prepared to
calculated the number I want it to, but given that this is what I need
to do, I don't see any better method than what I have above.  I would
be fairly content with that if what I wanted simply wasn't possible,
but for the fact that Mathematica has no trouble calculating the above
function out to t~500000, much large than my  needs(For me, t is a
time parameter that I am fairly sure need not be looked at above
t=100;  I'm not sure about t~50, which is why I want to calculate for
these values), so I think it must be possible, especially when the
function clearly goes to 1.

Dennis