Path: news.mathworks.com!not-for-mail
From: "leo nidas" <bleonidas25@yahoo.gr>
Newsgroups: comp.soft-sys.matlab
Subject: Re: trapz, numerical integration
Date: Sun, 7 Jun 2009 16:53:01 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 124
Message-ID: <h0gr9d$6db$1@fred.mathworks.com>
References: <h048vl$2v$1@fred.mathworks.com> <h04bpo$2vo$1@fred.mathworks.com> <h0fudh$pfo$1@fred.mathworks.com> <h0gade$kfk$1@fred.mathworks.com>
Reply-To: "leo nidas" <bleonidas25@yahoo.gr>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1244393581 6571 172.30.248.37 (7 Jun 2009 16:53:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 7 Jun 2009 16:53:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1194208
Xref: news.mathworks.com comp.soft-sys.matlab:545349



> I had originally suggested the use of a Gauss-Laguerre
> numerical integration. If your function is this:
> 
> fun = @(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m.*exp(-t./m)
> 
> then it can be broken down into two pieces.
> 
>     exp(b0+b1.*t)./(1+exp(b0+b1.*t)).*1/m
> 
> and 
> 
>     exp(-t./m)
> 
> I'll look at the part which is not exponential to see if a
> Gauss-Laguerre will work well.
> 
> K =@(t) exp(b0+b1.*t)./(1+exp(b0+b1.*t))./m;
> ezplot(K)
> 
> It is not really that polynomial, although it is well behaved,
> with a sigmoidal shape.
> 
> I then tried a Gauss-Laguerre quadrature on this function. It
> is not polynomial enough that a low order Gauss-Laguerre
> works, and roots fails above 14 nodes to generate any higher
> order Gauss-Laguerre rules than that. (I should have used a
> better routine to solve for the Gauss rule. No matter, since
> there are better choices to be found.) Look instead at a
> numerical tool like quadgk.
> 
> quadgk(fun,0,30,'reltol',1e-12)
> ans =
>          0.584687911313577
> 
> quadgk(fun,0,100,'reltol',1e-12)
> ans =
>          0.584733311243336
> 
> quadgk(fun,0,200,'reltol',1e-12)
> ans =
>          0.584733311243339
> 
> Or, try quadl. Note that quadl used more than 1200
> function evaluations to get this result, but I had set
> the tolerance rather small.
> 
> [I,fev] = quadl(fun,0,200,1e-12)
> I =
>          0.584733311243339
> fev =
>         1248
>  
> The default tolerance for quadl gave me a reasonable
> result at a cost of only 168 function evaluations. In fact
> this is far better accuracy than I would get from trapz
> anyway.
> 
> [I,fev] = quadl(fun,0,200)
> I =
>          0.584733310572222
> fev =
>    168
> 
> How well does trapz do? At a cost of 3000 function
> evaluations, you only get about 4 digits of accuracy,
> whereas the tools like quadl or quadgk will yield
> twice as many correct digits for vastly less work.
> 
> t = 0:0.01:30;
> trapz(t,fun(t))
> ans =
>          0.584687862068707
> 
> So I'm not at all sure why you have chosen to use
> trapz here, perhaps only that you did not know these
> other tools exist in MATLAB. Trapezoidal integration
> is rarely a good choice.
> 
> Can we do better yet? Certainly so, if I knew that your
> function would always have the same essential character
> of a sigmoid multiplied by a negative exponential. 
> 
> I'll stop here for now, unless you wanted me to explore
> this last avenue more deeply. Ask again if so.
> 
> John



Thanx again John,

It's me again.

You were right that I didn't know these other tools in MATLAB. I tried them though and with trial and error I think that quadgk is the right choise. Although it didn't work for the next:

setb0=-4.5313
setb1=2.7481
mhat=2.6470
pc=quadgk(@(gr) exp(setb0+setb1.*gr)./(1+exp(setb0+setb1.*gr)).*1/mhat.*exp(-gr./mhat),0,Inf,'reltol',1e-12)

For small changes of the values setb0 setb1 mhat it might work. As you can suppose setb0, setb1 mhat are random variables and I get a value of them in each iteration of the code. Then I try to integrate the above function which is the (sigmoidal*exponential distribution) . 

To tell you the whole problem I am interested in solving the following integrals for whatever values of setb0 setb1, and  mhat>0.

Integrals for the exponential distribution:

pc=  quadgk(@(gr) exp(setb0+setb1*gr)/(1+exp(setb0+setb1*gr))*1/mhat*exp(-gr/mhat),0,Inf);

muc1= quadgk(@(gr) exp(setb0+setb1*gr)/((1+exp(setb0+setb1*gr))^2)*1/mhat*exp(-gr/mhat),0,Inf); 

muc2= quadgk(@(gr) gr*exp(setb0+setb1*gr)/((1+exp(setb0+setb1*gr))^2)*1/mhat*exp(-gr/mhat),0,Inf);

As you can see, in all the above there is a product (something * exponential distribution).

The integrals for another distribution (    specifically the extended generalized gamma with p.d.f.=(1/gamma(g)).*(g.^g)*(l.^(a.*g)).*(x.^(a.*g-1)).*exp(-g.*((l.*x).^a))   )    where  -Inf<a<Inf , l>0, g>0, x>0.

are the corresponding ones that come if one substitute the pdf of exponential with the pdf of the GG. 

Any help for the above will be more than welcome...

Maybe if you can guide me through the example I gave you on top I might find me way through to the rest...

Thanx again for your time John!