Thread Subject:

Integration

Hi,

I have a simple question concerning numerical integration.

I need to numerical integrate a vector and use the trapz comment. However, I do not understand of how to do the non uniform spacing.

For example, when I run the following program to compute the CDF of a exponential distribution.

lambda(1)=1/8;
Num=10000;

for i=1:Num
      y(i,1) = i*10/Num;
      z(i,1) = exppdf(y(i,1),(1./lambda(1,1)));
end

trapz(y,z)

Why do I not get 1?

I appreciate any help!

Best

"steingre" wrote in message <ji1bhj$cvt$1@newscl01ah.mathworks.com>...
> lambda(1)=1/8;
> Num=10000;
> for i=1:Num
> y(i,1) = i*10/Num;
> z(i,1) = exppdf(y(i,1),(1./lambda(1,1)));
> end
> trapz(y,z)
>
> Why do I not get 1?
- - - - - - - - -
  You haven't gone far enough out in the range of y to approximate infinity. You have stopped where y is equal to 10 which makes exppdf equal to 1/8*exp(-10/8) = .0358 . This is not nearly far enough. Also you started at y = 0.001 instead of zero, another source of error. Finally, you should not expect exact answers because trapezoidal integration is only a numerical approximation to the ideal integral value.

  By the way, you *are* using uniform spacing here. Each interval in y is 0.001 .

Roger Stafford

"Roger Stafford" wrote in message <ji1q4g$on4$1@newscl01ah.mathworks.com>...
> You haven't gone far enough out in the range of y to approximate infinity. You have stopped where y is equal to 10 which makes exppdf equal to 1/8*exp(-10/8) = .0358 . This is not nearly far enough.
- - - - - - - -
  I made a mistake in calculation, the furthest value is 8*exp(-80) which is very small and should be far enough out. Whatever error you are getting is probably due to starting at .001. Besides correcting this, try taking smaller steps to see the effect. How large is your error at present?

Roger Stafford

"Roger Stafford" wrote in message <ji1qsr$qmr$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <ji1q4g$on4$1@newscl01ah.mathworks.com>...
> > You haven't gone far enough out in the range of y to approximate infinity. ...
> I made a mistake in calculation,....
- - - - - - - - -
  Blush! My first post was correct and my second one in error. You haven't gone far enough out after all. The last value of z is

 1/(1/lambda)*exp(-10/(1/lambda)) = 1/8*exp(-10/8) = .0358

That's too large.

Roger Stafford

I started at 0 and go until 10 with steps 0.01 and it works!

lambda(1)=1/10;
Num=1000;

x(:,1) = (exprnd((1./lambda(1,1)),Num,1));
for i=1:Num
      y(i,1) = (i-1)*100/Num;
      z(i,1) = exppdf(y(i,1),(1./lambda(1,1)));
end
trapz(y,z)


However, do you know a method of how I get the optimal spacing?

For example, I change lambda to 1/100 and I do not get 1 for the PDF only if I also change y.

"steingre" wrote in message <ji2p84$p2u$1@newscl01ah.mathworks.com>...
> I started at 0 and go until 10 with steps 0.01 and it works! ......
- - - - - - - -
  The reason it worked is that you went until y equals 100, not 10. If you had only gone to 10 that is only at the new mean value point, which is way short. A 100 is safely past there. When you change to lambda = 1/100, that sets the mean value to 100, so you would then need to go far past there.

  You would get a better feel for what you are doing if you made plots of z against y ahead of time. You would see why some settings are okay and others not so good.

 plot(y,z)

  You need to take into consideration the nature of your integrand as you alter the mean parameter. As you have it set up, 'lambda' is the reciprocal of the mean of this distribution. The smaller 'lambda', the broader the distribution and therefore the further out you need to integrate for approximating an infinite limit.

Roger Stafford

"Roger Stafford" wrote in message <ji3d3o$79e$1@newscl01ah.mathworks.com>...
> "steingre" wrote in message <ji2p84$p2u$1@newscl01ah.mathworks.com>...
> > I started at 0 and go until 10 with steps 0.01 and it works! ......
> - - - - - - - -
> The reason it worked is that you went until y equals 100, not 10. If you had only gone to 10 that is only at the new mean value point, which is way short. A 100 is safely past there. When you change to lambda = 1/100, that sets the mean value to 100, so you would then need to go far past there.
>
> You would get a better feel for what you are doing if you made plots of z against y ahead of time. You would see why some settings are okay and others not so good.
>
> plot(y,z)
>
> You need to take into consideration the nature of your integrand as you alter the mean parameter. As you have it set up, 'lambda' is the reciprocal of the mean of this distribution. The smaller 'lambda', the broader the distribution and therefore the further out you need to integrate for approximating an infinite limit.
>
> Roger Stafford

Thank you for the clarification!