Thread Subject:

Numerical Integration of Frechet: Precision problem

Hi,

I try to numerical integrate a Frechet distribution and have the following problem of precision. In particular when choosing a very low show parameter close to 1. If I chose k>2 then the precision problem disappears.

I appreciate any help!

Here is an example of calculating the mean:

% Define number of points;
Num=10000;
for i=1:Num
      y(i,1) = (i)*1000/Num;
end

% Shape parameter of the revenue distribution
k=1.19;

% Scale parameter
s=1/10;

% Random variable
x=y;

% PDF of Frechet
pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
trapz(x,pdf)

% Empircial mean of Frechet
z=x.*pdf;
trapz(x,z)

% Actual mean of Frechet
mean=1./s.*gamma(1-1./k)

On 22 Feb., 16:12, "steingre " <wal...@steingress.com> wrote:
> Hi,
>
> I try to numerical integrate a Frechet distribution and have the following problem of precision. In particular when choosing a very low show parameter close to 1. If I chose k>2 then the precision problem disappears.
>
> I appreciate any help!
>
> Here is an example of calculating the mean:
>
> % Define number of points;
> Num=10000;
> for i=1:Num
>       y(i,1) = (i)*1000/Num;
> end
>
> % Shape parameter of the revenue distribution
> k=1.19;
>
> % Scale parameter
> s=1/10;
>
> % Random variable
> x=y;
>
> % PDF of Frechet
> pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
> trapz(x,pdf)
>
> % Empircial mean of Frechet
> z=x.*pdf;
> trapz(x,z)
>
> % Actual mean of Frechet
> mean=1./s.*gamma(1-1./k)

If - as in your case - the function you want to integrate is
explicitely known,
don't use trapz, but MATLAB's quadgk. It distributes the x-mesh
adaptively
such that the error of integration is below a specified tolerance.

Best wishes
Torsten.

Torsten <Torsten.Hennig@umsicht.fraunhofer.de> wrote in message <2cb6e221-6635-494a-baa1-11d271aeffdb@o6g2000vbz.googlegroups.com>...
> On 22 Feb., 16:12, "steingre " <wal...@steingress.com> wrote:
> > Hi,
> >
> > I try to numerical integrate a Frechet distribution and have the following problem of precision. In particular when choosing a very low show parameter close to 1. If I chose k>2 then the precision problem disappears.
> >
> > I appreciate any help!
> >
> > Here is an example of calculating the mean:
> >
> > % Define number of points;
> > Num=10000;
> > for i=1:Num
> >       y(i,1) = (i)*1000/Num;
> > end
> >
> > % Shape parameter of the revenue distribution
> > k=1.19;
> >
> > % Scale parameter
> > s=1/10;
> >
> > % Random variable
> > x=y;
> >
> > % PDF of Frechet
> > pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
> > trapz(x,pdf)
> >
> > % Empircial mean of Frechet
> > z=x.*pdf;
> > trapz(x,z)
> >
> > % Actual mean of Frechet
> > mean=1./s.*gamma(1-1./k)
>
> If - as in your case - the function you want to integrate is
> explicitely known,
> don't use trapz, but MATLAB's quadgk. It distributes the x-mesh
> adaptively
> such that the error of integration is below a specified tolerance.
>
> Best wishes
> Torsten.

Hi,

Thank you for your post! The reason why I use a known pdf is that I want to learn more about trapz. I have a data vector where I do not know the distribution. I estimate the cdf with ecdf and then want to integrate over that cdf. First, I test for a known CDF, here the CDF of Frechet:

clear all;

% Define number of points;

Num=100000;

% Example of numerical integration that works

% Getting the mean form a random variable knowing the pdf

x(:,1) = (exprnd((1./lambda(1,1)),Num,1));
x=sort(x);

% **************************************************************

% Doing the same for Frechet

% Shape parameter
k=2;

% Scale parameter (s=1/s in Wikipedia)
s=1;

x=(x.^k);

% Known PDF of Frechet
pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));

% Check whether pfd integrates to 1
trapz(x,pdf);

% Known CDF of Frechet
cdf=exp(-(x.*s).^(-(k)));

% Empirical mean
mean=trapz(x,(1-cdf))

% Actual mean of Frechet
mean=1./s.*gamma(1-1./k)


Then I want to compare this with a CDF I estimate from data:

[F z] = ecdf(x,'function','cdf');
mean=trapz(z,(1-F))

But this does not work. I should get the same mean!

"steingre" wrote in message <ji30kb$lri$1@newscl01ah.mathworks.com>...
> Hi,
>
> I try to numerical integrate a Frechet distribution and have the following problem of precision. In particular when choosing a very low show parameter close to 1. If I chose k>2 then the precision problem disappears.
>
> I appreciate any help!
>
> Here is an example of calculating the mean:
>
> % Define number of points;
> Num=10000;
> for i=1:Num
> y(i,1) = (i)*1000/Num;
> end
>
> % Shape parameter of the revenue distribution
> k=1.19;
>
> % Scale parameter
> s=1/10;
>
> % Random variable
> x=y;
>
> % PDF of Frechet
> pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
> trapz(x,pdf)
>
> % Empircial mean of Frechet
> z=x.*pdf;
> trapz(x,z)
>
> % Actual mean of Frechet
> mean=1./s.*gamma(1-1./k)
- - - - - - - - -
  You picked a tough function to try 'trapz' on. The integrand which you called z = x.*pdf is a very difficult one to do numerical integration on. Near x = 0 it rises and drops rapidly and therefore would need very close spacing of x points in that vicinity, and yet it needs a rather large value for x before it becomes small enough to stop the integration even though the points could be spaced far apart there. In other words it very much needs non-uniform spacing.

  Instead of trying to find an optimum kind of non-uniform spacing I would recommend appropriate changes of variable. First define

 u = (s*x)^(-k)

Then manipulating the differentials properly we get

 pdf(x)*dx = -exp(-u)*du

and this integral from infinity to 0 is obviously 1 from elementary calculus (which demonstrates that Frechet had the normalizing factor right.)

  To find the mean value I would suggest the further change of variable

 v = u^(1-1/k)

This leads to the following:

 x*pdf(x)*dx = -k/(k-1)/s*exp(-v^(k/(k-1)))*dv

The right side here is something 'trapz' can handle accurately. To see that just run the following:

 clear
 format long
 k=1.19;
 s=1/10;
 v = linspace(0,2,8192); % <-- 2 is far enough
 m = k/(k-1)/s*exp(-v.^(k/(k-1)));
 trapz(v,m)
 1/s*gamma(1-1/k)

On my machine these both display as: 58.24171302007969 .

  To see how much nicer this latter form is for 'trapz', do a plot of v against m and compare that with a plot of x against x*pdf(x).

Roger Stafford

"Roger Stafford" wrote in message <ji6qu8$m3s$1@newscl01ah.mathworks.com>...
> "steingre" wrote in message <ji30kb$lri$1@newscl01ah.mathworks.com>...
> > Hi,
> >
> > I try to numerical integrate a Frechet distribution and have the following problem of precision. In particular when choosing a very low show parameter close to 1. If I chose k>2 then the precision problem disappears.
> >
> > I appreciate any help!
> >
> > Here is an example of calculating the mean:
> >
> > % Define number of points;
> > Num=10000;
> > for i=1:Num
> > y(i,1) = (i)*1000/Num;
> > end
> >
> > % Shape parameter of the revenue distribution
> > k=1.19;
> >
> > % Scale parameter
> > s=1/10;
> >
> > % Random variable
> > x=y;
> >
> > % PDF of Frechet
> > pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
> > trapz(x,pdf)
> >
> > % Empircial mean of Frechet
> > z=x.*pdf;
> > trapz(x,z)
> >
> > % Actual mean of Frechet
> > mean=1./s.*gamma(1-1./k)
> - - - - - - - - -
> You picked a tough function to try 'trapz' on. The integrand which you called z = x.*pdf is a very difficult one to do numerical integration on. Near x = 0 it rises and drops rapidly and therefore would need very close spacing of x points in that vicinity, and yet it needs a rather large value for x before it becomes small enough to stop the integration even though the points could be spaced far apart there. In other words it very much needs non-uniform spacing.
>
> Instead of trying to find an optimum kind of non-uniform spacing I would recommend appropriate changes of variable. First define
>
> u = (s*x)^(-k)
>
> Then manipulating the differentials properly we get
>
> pdf(x)*dx = -exp(-u)*du
>
> and this integral from infinity to 0 is obviously 1 from elementary calculus (which demonstrates that Frechet had the normalizing factor right.)
>
> To find the mean value I would suggest the further change of variable
>
> v = u^(1-1/k)
>
> This leads to the following:
>
> x*pdf(x)*dx = -k/(k-1)/s*exp(-v^(k/(k-1)))*dv
>
> The right side here is something 'trapz' can handle accurately. To see that just run the following:
>
> clear
> format long
> k=1.19;
> s=1/10;
> v = linspace(0,2,8192); % <-- 2 is far enough
> m = k/(k-1)/s*exp(-v.^(k/(k-1)));
> trapz(v,m)
> 1/s*gamma(1-1/k)
>
> On my machine these both display as: 58.24171302007969 .
>
> To see how much nicer this latter form is for 'trapz', do a plot of v against m and compare that with a plot of x against x*pdf(x).
>
> Roger Stafford

Roger, great solution! Thank you!

"steingre" wrote in message <ji8n6s$ag$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <ji6qu8$m3s$1@newscl01ah.mathworks.com>...
> > "steingre" wrote in message <ji30kb$lri$1@newscl01ah.mathworks.com>...
> > > Hi,
> > >
> > > I try to numerical integrate a Frechet distribution and have the following problem of precision. In particular when choosing a very low show parameter close to 1. If I chose k>2 then the precision problem disappears.
> > >
> > > I appreciate any help!
> > >
> > > Here is an example of calculating the mean:
> > >
> > > % Define number of points;
> > > Num=10000;
> > > for i=1:Num
> > > y(i,1) = (i)*1000/Num;
> > > end
> > >
> > > % Shape parameter of the revenue distribution
> > > k=1.19;
> > >
> > > % Scale parameter
> > > s=1/10;
> > >
> > > % Random variable
> > > x=y;
> > >
> > > % PDF of Frechet
> > > pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
> > > trapz(x,pdf)
> > >
> > > % Empircial mean of Frechet
> > > z=x.*pdf;
> > > trapz(x,z)
> > >
> > > % Actual mean of Frechet
> > > mean=1./s.*gamma(1-1./k)
> > - - - - - - - - -
> > You picked a tough function to try 'trapz' on. The integrand which you called z = x.*pdf is a very difficult one to do numerical integration on. Near x = 0 it rises and drops rapidly and therefore would need very close spacing of x points in that vicinity, and yet it needs a rather large value for x before it becomes small enough to stop the integration even though the points could be spaced far apart there. In other words it very much needs non-uniform spacing.
> >
> > Instead of trying to find an optimum kind of non-uniform spacing I would recommend appropriate changes of variable. First define
> >
> > u = (s*x)^(-k)
> >
> > Then manipulating the differentials properly we get
> >
> > pdf(x)*dx = -exp(-u)*du
> >
> > and this integral from infinity to 0 is obviously 1 from elementary calculus (which demonstrates that Frechet had the normalizing factor right.)
> >
> > To find the mean value I would suggest the further change of variable
> >
> > v = u^(1-1/k)
> >
> > This leads to the following:
> >
> > x*pdf(x)*dx = -k/(k-1)/s*exp(-v^(k/(k-1)))*dv
> >
> > The right side here is something 'trapz' can handle accurately. To see that just run the following:
> >
> > clear
> > format long
> > k=1.19;
> > s=1/10;
> > v = linspace(0,2,8192); % <-- 2 is far enough
> > m = k/(k-1)/s*exp(-v.^(k/(k-1)));
> > trapz(v,m)
> > 1/s*gamma(1-1/k)
> >
> > On my machine these both display as: 58.24171302007969 .
> >
> > To see how much nicer this latter form is for 'trapz', do a plot of v against m and compare that with a plot of x against x*pdf(x).
> >
> > Roger Stafford
>
> Roger, great solution! Thank you!

Hi,

In fact the method works great for the PDF. Now, I want to integrate over the cdf and the same change of variables does not give me a good approximation.

I guess the magic question is: How do I chose a good change of variables? Is there any rule?

Best,
steingre

Here is the original program

% Define number of points;
 Num=100000;
 for i=1:Num
       y(i,1) = (i)*1000/Num;
 end
 
 % Shape parameter of the revenue distribution
 k=1.19;
  
 % Scale parameter
 s=1/1;
 
 % Random variable
 x=y;
 
 % PDF of Frechet
 pdf=s.*k.*(x.*s).^(-(1+k)).*exp(-(x.*s).^(-(k)));
 trapz(x,pdf)
 
 % Empircial mean of Frechet
 z=x.*pdf;
 trapz(x,z)

 % Alternative empirical mean from CDF of Frechet
    
 cdf = exp(-(x./s).^(-(k)));
 
 trapz(x,1-cdf)

"steingre" wrote in message <ji8oe8$4f5$1@newscl01ah.mathworks.com>...
> In fact the method works great for the PDF. Now, I want to integrate over the cdf and the same change of variables does not give me a good approximation.
>
> I guess the magic question is: How do I chose a good change of variables? Is there any rule?
> ........
> k=1.19;
> s=1/1;
> % Alternative empirical mean from CDF of Frechet
> cdf = exp(-(x./s).^(-(k)));
> trapz(x,1-cdf)
- - - - - - - - -
  To integrate the cdf I think we have to first resort to integration by parts to help 'trapz' along. First a scale change: u = x/s gives us

 int('s*(1-exp(-u^(-k))),'u',0,inf)

  Integrating this by parts gives s*u*(1-exp(-u^(-k))) evaluated at u = inf minus its value at u = 0, together with

 -int('s*u*d(1-exp(-u^(-k)))/du,'u',0,inf)

Good old l'Hospital's rule shows us that the first part is zero at both u = inf and u = 0, so we are left with just this last integral. Doing the differentiation gives:

 int('s*k*u^(-k)*exp(-u^(-k))','u',0,inf)

Finally we change variables: v = u^(-(k-1)) and obtain

 int('s*k/(k-1)*exp(-v^(k/(k-1))),'v',0,inf)

Again this is something our delicate 'trapz' can handle. In fact it is the same function as before except for the different value of s. We have simply arrived at it by a different path.

 clear
 format long
 s = 1;
 k = 1.19;
 v = linspace(0,2,8192);
 f = s*k/(k-1)*exp(-z.^(k/(k-1)));
 trapz(v,f) = 5.824171302007956

  You may consider it cheating to use integration by parts to assist 'trapz' in this manner but otherwise it would have had severe trouble handling something like the Frechet cdf. At least I couldn't think of a change of variable that would enable accurate integration by itself.

Roger Stafford

"Roger Stafford" wrote in message <ji95qb$jeh$1@newscl01ah.mathworks.com>...
> "steingre" wrote in message <ji8oe8$4f5$1@newscl01ah.mathworks.com>...
> > In fact the method works great for the PDF. Now, I want to integrate over the cdf and the same change of variables does not give me a good approximation.
> >
> > I guess the magic question is: How do I chose a good change of variables? Is there any rule?
> > ........
> > k=1.19;
> > s=1/1;
> > % Alternative empirical mean from CDF of Frechet
> > cdf = exp(-(x./s).^(-(k)));
> > trapz(x,1-cdf)
> - - - - - - - - -
> To integrate the cdf I think we have to first resort to integration by parts to help 'trapz' along. First a scale change: u = x/s gives us
>
> int('s*(1-exp(-u^(-k))),'u',0,inf)
>
> Integrating this by parts gives s*u*(1-exp(-u^(-k))) evaluated at u = inf minus its value at u = 0, together with
>
> -int('s*u*d(1-exp(-u^(-k)))/du,'u',0,inf)
>
> Good old l'Hospital's rule shows us that the first part is zero at both u = inf and u = 0, so we are left with just this last integral. Doing the differentiation gives:
>
> int('s*k*u^(-k)*exp(-u^(-k))','u',0,inf)
>
> Finally we change variables: v = u^(-(k-1)) and obtain
>
> int('s*k/(k-1)*exp(-v^(k/(k-1))),'v',0,inf)
>
> Again this is something our delicate 'trapz' can handle. In fact it is the same function as before except for the different value of s. We have simply arrived at it by a different path.
>
> clear
> format long
> s = 1;
> k = 1.19;
> v = linspace(0,2,8192);
> f = s*k/(k-1)*exp(-z.^(k/(k-1)));
> trapz(v,f) = 5.824171302007956
>
> You may consider it cheating to use integration by parts to assist 'trapz' in this manner but otherwise it would have had severe trouble handling something like the Frechet cdf. At least I couldn't think of a change of variable that would enable accurate integration by itself.
>
> Roger Stafford

Thank you Roger! Indeed, to find a change of variable to integrate over the CDF doesn't seem to be an easy task. Do you know any other method that allows me to accurately integrated over the Frechet CDF?

Again, thank you for precious help!

clear all;

Num=100000;
x(:,1) = (exprnd((1./lambda(1,1)),Num,1));

% Shape parameter of the revenue distribution
k=1.19;

% Scale parameter (s=1/s in Wikipedia)
s=lambda

% The random variable distributed as a Frechet with parameter k and s

t=(x.^(-1/k));
tt=sort(t);

% Integrate over the sample CDF:

format long

% Traditional uniform spacing
xx = 0:0.001:max(tt);

% Log S spacing
a=quantile(tt,10^(-3))
b=quantile(tt,0.999)

a=log10(a);
b=log10(b);

xl = logspace(a,b,Num);


F = exp(-(xx./s).^(-(k)));
Fl = exp(-(xl./s).^(-(k)));

trapz(xx,F)
trapz(xl,Fl)