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?