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