Thread Subject:

Innaccuracy in a gammainc and , how best to override?

Here is a function

function out = exact_d_cgf(t,lambda,delta)
A = - lambda ./ log(delta);

out = real( (-A ./ t) - (exp(t) .* (-t) .^ (A-1) ) ./
( gamma(A)*gammainc(-t,A,'lower'))) ;


end

I know, for theoretical reasons, that lim t->inf, exact_d_cgf -> 1.
However, for t>~ 30, there is inaccuracy in this calculation; I
imagine it has to do with gammainc(-t,A,'lower') ->inf and (-A ./ t)
- (exp(t) .* (-t) .^ (A-1) ) ->inf at the same time with t.

For large t, this begins to fluctuate and introduces errors. What is
the best way to override this? I tried somethings like

out = min([ real(...) 1])

but the problem is the innaccuracy goes either way; sometimes it
returns arbirarilty large answers, and sometimes it returns small
answers. What is my best bet? Just force out = 1 for t>tbar, where
I just pick tbar by hand?

Dennis

DOD <dcodea@gmail.com> wrote in message <e6d2bf35-7011-4fc8-baf4-45138f5e8798@a7g2000yqb.googlegroups.com>...
> .....
> out = real( (-A ./ t) - (exp(t) .* (-t) .^ (A-1) ) ./
> ( gamma(A)*gammainc(-t,A,'lower'))) ;
> ......
> I know, for theoretical reasons, that lim t->inf, exact_d_cgf -> 1.
> ......
- - - - - - - - -
  You say "lim t->inf". If you mean that, your first 'gammainc' argument would lie in the negative range! I am not sure it is defined for negative arguments, and it certainly becomes non-integrable as this argument approaches minus infinity. Also your numerator, (-t).^(A-1), would be rather ill-defined as a possibly irrational power of a negative quantity.

  On the other hand if you actually meant lim t->-inf with -t a positive quantity, I believe you will find that your expression should approach zero, not one - the value of 'gammainc' would approach 1, while the above numerator would approach 0, and -A./t would also approach 0, giving you a limiting difference of 0.

  Also you speak of 'gammainc' approaching infinity, but my understanding of it is that it is normalized so as to always lie between 0 and 1.

  Furthermore you speak of the single quantity

 "(-A./t)-(exp(t).*(-t).^(A-1))"

approaching infinity but your parentheses for the 'out' expression are inconsistent with that as a part of the calculation.

  In other words you have me thoroughly confused, DOD. Perhaps you could explain your problem again in a more detailed manner.

Roger Stafford

On Sep 21, 2:37 pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> DOD <dco...@gmail.com> wrote in message <e6d2bf35-7011-4fc8-baf4-45138f5e8...@a7g2000yqb.googlegroups.com>...
> > .....
> > out = real( (-A ./ t)  - (exp(t) .* (-t) .^ (A-1) ) ./
> > ( gamma(A)*gammainc(-t,A,'lower'))) ;
> > ......
> > I know, for theoretical reasons, that lim t->inf,  exact_d_cgf -> 1.
> > ......
>
> - - - - - - - - -
>   You say "lim t->inf".  If you mean that, your first 'gammainc' argument would lie in the negative range!  I am not sure it is defined for negative arguments, and it certainly becomes non-integrable as this argument approaches minus infinity.  Also your numerator, (-t).^(A-1), would be rather ill-defined as a possibly irrational power of a negative quantity.
>
>   On the other hand if you actually meant lim t->-inf with -t a positive quantity, I believe you will find that your expression should approach zero, not one - the value of 'gammainc' would approach 1, while the above numerator would approach 0, and -A./t would also approach 0, giving you a limiting difference of 0.
>
>   Also you speak of 'gammainc' approaching infinity, but my understanding of it is that it is normalized so as to always lie between 0 and 1.
>
>   Furthermore you speak of the single quantity
>
>  "(-A./t)-(exp(t).*(-t).^(A-1))"
>
> approaching infinity but your parentheses for the 'out' expression are inconsistent with that as a part of the calculation.
>
>   In other words you have me thoroughly confused, DOD.  Perhaps you could explain your problem again in a more detailed manner.
>
> Roger Stafford

Allow me to walk you through every single step that leads to this
expression, then.

Consider a random variable X with the following PDF:


-((L x^(-1 - L/ Log[d] ))/ Log[d] )

0<L, d<1

This is the simply the pdf of d^tau, where tau is exponential with
parameter L,

Consider next the moment generating function of X

Integrate[Exp(t x ) (-((L x^(-1 - L/ Log[d] ))/ Log[d] )) ,{x,0,1}]

Here I use Mathematica notation, because that's how I calculated it.
x is only defined from 0 to 1, since d^tau, for d <1, can only take
that range.

Finally, the culmulant generating function of X, is Log(MGF)

Now, if I plug this into Mathematica, I get the following expression
for the MGF. Let A = - L / log d. then the MGF is

- A (-t)^-A (Gamma[A] - Gamma[A, -t])

However, due to differences in how Mathematica and Matlab implement
the incomplete gamma function, this is calculated in matlab as

- A (-t)^-A (gamma(A)*gammainc(-t,A,'lower')


So that's the MGF. I don't know what to tell you other than, "it is".
what I am interested in the the CGF, and what this function calculates
is the derivative of the CGF. Happly, the derivative of Log(MGF) is
just MGF' / MGF. It simplifies to

(-A + (exp(t).* (-t).^A)/ (gamma(A).*gammainc(-t,A,'lower'))./t

So, that is the function I am interested in. I suppose you could tell
me some more that I shouldn't be, that there must be some other
function I should be interested in, but let's assume that this is in
fact what I am interested in. Maybe if you ran it for a few values of
t>0, you would be convinced that it's real part actually goes to 1?
  Give it a try with A = .0267, t= 20. I get 0.9471 + 0.0000i with
R2011b. With t = 50, I get -5.5700e-004 +1.9981e-025i, and I am
convinced this is do to roundoff or something in that ratio.
Mathematica has no such trouble, and I'm sure that's because it is
holding my hand through the numerical issues. Is there a way to
accomplish that in Matlab?

DOD <dcodea@gmail.com> wrote in message <7583ba5c-88b2-48a7-8f24-f341f096c335@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

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

DOD <dcodea@gmail.com> wrote in message <0963bec6-df81-4dd9-857e-6c255fc489d4@15g2000vbb.googlegroups.com>...
> 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)
> .........
- - - - - - - - -
  No, you are mistaken, Dennis. You have forgotten the necessary Jacobian in making the change of variable:

 u = x^A
 x = u^(1/A)
 du = A*x^(A-1)*dx

 int(exp(t*x)*A*x^(A-1)*dx,0,1) =
 int(ext(t*u^(1/A))*du,0^A,1^A) =
 int(ext(t*u^(1/A))*du,0,1) =

Roger Stafford

On Sep 22, 10:52 am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> DOD <dco...@gmail.com> wrote in message <0963bec6-df81-4dd9-857e-6c255fc48...@15g2000vbb.googlegroups.com>...
> > 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)
> > .........
>
> - - - - - - - - -
>   No, you are mistaken, Dennis.  You have forgotten the necessary Jacobian in making the change of variable:
>
>  u = x^A
>  x = u^(1/A)
>  du = A*x^(A-1)*dx
>
>  int(exp(t*x)*A*x^(A-1)*dx,0,1) =
>  int(ext(t*u^(1/A))*du,0^A,1^A) =
>  int(ext(t*u^(1/A))*du,0,1) =
>
> Roger Stafford

Ah, I see. I think you'd find int and quad still have trouble with
that. If anyone out there is interested, I have found a solution to
this problem, using exponential integrals. The MGF has the following
form, according to mathematica:

(-t)^-A Gamma[1 + A] - A ExpIntegralE[1 - A, -t]

This is implemented in matlab as
(-t)^(-A) * gamma(1 + A) - A * mfun('Ei',1-A,-t)

the derivative wrt to t is
-A (-t)^(-A -1) - A * mfun('Ei',-A,-t)

there is no problem here with accuracy, so when I evaluate
(-A (-t)^(-A -1) - A * mfun('Ei',-A,-t)) ./ ((-t)^(-A) * gamma(1 + A)
- A * mfun('Ei',1-A,-t))

I get accurate answers (in my opinion) up to t~ 1000, good enough for
me!

Thanks for your time and help with this

Dennis