A Real Application of Taylor Series

MTHSC108H
Daniel D. Warner
December 3, 1999

There is a function in statistics called the Student-t probability distribution. The most important factor is
, where stands for the number of degrees of freedom. It is not uncommon for to

be a very large number. If is less than 1 and is very large, then as we saw in class we may find that we get
an incorrect answer because we are not carrying enough digits.

Here's an example, showing how "catostrophic underflow" can occur.

> g1 := unapply((1+x^2/nu), x, nu);
h1 := unapply(g1(x,nu)^((nu+1)/2), x, nu);
g1(0.001, 10^7);
h1(0.001, 10^7);

With Maple we can ask for more digits and get the correct answer -- if we happen to notice that we had
a problem.

> Digits := 20;
g1(0.001, 10^7);
h1(0.001, 10^7);
evalf(%,10);

On other calculators and more numerically oriented computational programs asking for more digits is
not an option. However, Taylor series can come to the rescue . The way to calculate our function is
simply as . Now we can construct a series for the log, do the multiplication
and take the exponential. Start by letting .

> Digits := 10;
f := ln(1+y);
t7 := taylor(f,y,7);
p6 := convert(t7,polynom);

Since , we have

> g := simplify((x^2+y)/(2*y)*p6);
h := exp(g);
h2 := unapply(h,x,y);

Let's try a few values.

> xx := 0.0012345;
for k from 1 to 5 do
[h1(xx,10^k),h2(xx,(xx^2/10^k))]
od;

The function h2 is certainly changing more smoothly, but how accurate is it?

Taylor's Theorem with Remainder (for three terms) states

+ + -
The integral at the end is the remainder. There are four steps to finding a bound for the remainder. Let denote the k+1 derivative
of and let be a number for which for all in the interval from to . Also let denote the Taylor polynomial
of degree k. Then we have the following string of inequalities.






.

Now let's try to apply Taylor's Theorem with Remainder to the approximation we developed earlier, where .

> g[7] := diff(f,y$7);

Clearly gets smaller as gets larger. If we are only interested in positive values of , which happens to be
the case for the application under consideration, then we can let . Therefore

.

This means that the error is less than whenever y is less than

> evalf(root(7!*0.5*10^(-10)/720,7));

Since , we see that for and we have an error less than

> x := 2;
nu := 100;
evalf(subs(y=(x^2/nu),720*y^7/7!),10);

To get a little better idea of the accuracy of approximation we can do the calculations in Maple to many more digits
and plot the actual error.

> Digits := 20;
plot(h1(1,1/y)-h2(1,y),y=0..0.05);

>

Going back to the small table where and took on the values 10, 100, 1000, 10000, and 100000, we see that
the approximation produced the correctly rounded answer in every case and the straight forward formula did not. We can
confirm this by rerunning the calculation with more digits. Notice how much the answers for change in the first 10 digits.

> Digits := 15;
xx := 0.0012345;
for k from 1 to 5 do
[h1(xx,10^k),h2(xx,(xx^2/10^k))]
od;

Notice that in this last calculation, the answers do agree to the first 10 digits - but they don't agree to all the digits shown.
Which are the more corrent? Can you confirm your answer?