[S] About uniroot(f,...) when f is an integration [long & cranky]

Bill Venables (wvenable@attunga.stats.adelaide.edu.au)
Sat, 30 May 1998 12:11:32 +0930


[An S-news reader] writes:
> ------------------------------------------------------------------
> mx<-300; my<-50; a0<-30
> b<-1.5;a<-gamma(1+1/b)/mx;
> d<-my/(gamma(1+1/b)*a^b*b*a0)
> my.fucn <- function(c)
> {
> as.numeric(c)
> weibx<-function(x) {x^(b-1)*exp(-(a*x)^b+c*x)}
> d-integrate(weibx,0,Inf)[1]
> }
> uniroot(my.func,c(1,100))[1]
>
> Error in c * x: Non-numeric first operand
> Dumped
> ----------------------------------------------------------------

I really do not want to single this (all too non-unique) user out
for harsh treatment, but it seems to me this is an excellent if
negative example on which to make a couple of points about how to
approach such a problem and how to approach S-news.

For the given values for parameters a, b and c the integral does
not converge. This would have been obvious if the user had
simply evaluated the integrand, weibx, for given a and b, c
anywhere between 1 and 100 and for x = 1:200, say. So that early
check has not been done.

Supposing that the integrand problem could be fixed, integrate()
may have worked at the command-line level but it would certainly
not have worked inside my.func() since the value of the parameter
`c' is not visible to the integration routine from the frame of
the function my.func(). (This is the notorious S scoping rule
anomaly, and where S and R conspicuously differ.) It follows
that my.func() cannot have been checked, either.

Before calling uniroot() it is elementary that you should check
that the function has opposite signs at the end of the range.
This cannot have been done.

It is pretty clear what has been done in this case. The user has
encountered a simple enough looking mathematical problem, coded
up the entire thing in S with a mental come-back for anything
that didn't work (and with very little whitespace, note...)
pressed the button and got some obscure message about a
non-numeric argument. The first attempt to fix this has been to
put a desparate but futile `as.numeric(c)' statement into the
integrand. The next was to clip the transcript, edit it a bit
(making at least one typo in the process, my.func --> mu.fucn)
and lay it at the feet of S-news. (I may well be too harsh, but
this is what it appears to be, and it is a common methodology.)

May I suggest a better way to approach such problems is to check
the code step by step as you put it together. And by `check' I
mean *both* look carefully at the code and run it on the computer
on several test cases to assure yourself that it does do what you
think it should do.

In particular, in this case I think the user should have

1. Checked the integrand function, weibx, to see that the
integration looked feasible, first.

2. Checked the integrate step to see that the integral did work
in selected cases at least.

3. Checked that my.func() worked with several values of c.

4. Checked that my.func(c) gave values with opposite signs at the
ends of the range in which the zero is supposed to occur.

5. If the step by step checks all pass but the problem still does
not work and there is nothing more he could do, constructed a
full example transcript in a form others will find easy to
read and include it verbatim in an enquiry to be sent to S-news.

6. When the error message generated a "Dumped" notice, observed
the golden rule: Put in a call to traceback() and *looked* at
it. Include that in the transcript also if it does not
suggest further leads (which it very often will, of course).

7. Finally, contacted S-news at the first instance of
inexplicable behaviour, not the last.

--

Of course you might well decide to go for broke, string it all together first and try your luck. If it works you may be OK (but don't be too sure) but if it fails there is no real alternative to checking it, incrementally, *on the computer*. I have known many users whose debugging tactic seems to be to stare at the code, get angrier and angrier and finally become totally frustrated because the mother-so&so computer doesn't argue back. The "I tell you that code should work!" tactic, doesn't.

Uncle Bill.

[BTW I have already replied directly to this user in a much more constructive and friendly manner and I do sincerely hope his problem is solved.] ----------------------------------------------------------------------- This message was distributed by s-news@wubios.wustl.edu. To unsubscribe send e-mail to s-news-request@wubios.wustl.edu with the BODY of the message: unsubscribe s-news