[R] Area() artefacts??

[Prof Brian Ripley]

You seem to be be calculating a complex result and then taking the real 
part.  Do that inside fun(), and you can use integrate.

Did you think to look at the function you were trying to integrate 
numerically?

x <- seq(-pi, pi, length=1000)
plot(x, Re(fun(x, 65, 3, 0.5)), type = "l")

Not a pretty sight, and no wonder you need a small tolerance.
Numerical integration routines are designed for smooth functions, and the 
tolerance is related to how smooth.  I don't think the result is 'wrong', 
in fact a pretty good shot given that picture.

area() (sic) is a support function for the first (1994) edition of MASS, 
as the help page shows.  I think you failed to consult it as the posting 
guide asked you to, for this was covered there.


On Sat, 12 May 2007, Sergey Goriatchev wrote:

> Hello, everybody
>
> I run the following program, and depending on the size of eps I get
> different results.
> With eps=1e-05, the program calculates wrong values for x=65:67 and
> others. The program runs fine with eps=1e-07. Why is this so?
> Also, I am using area() instead of integrate() because I cannot make
> integrate to work, especially with imaginary numbers. Maybe someone
> can show me how to use integrate in this particular code? THanks in
> advance!
>
> #Computes the p.m.f. via (10.53) of the number of i.i.d. Ber(p) trials
> #required until m consecutive successes occur.
> #Requires MASS package
>
> #==========================================================#
> consecpmf <- function(xvec, m, p, eps=1e-05){
> library(MASS)
> f<-numeric()
> for(j in seq(xvec)){
> 	x <- xvec[j]
> 	f[j] <- area(fun, -pi, pi, limit=1000, eps=eps, x, m, p)
> }
> f<-Re(f)
> round(f,4)
> }
>
> fun <- function(t,x,m,p){
> I <- exp(-1i*t*x)*cf(t,m,p)/(2*pi)
> I
> }
>
> cf <- function(t,m,p){
> q <- 1-p
> if(m==1) {g <- p*exp(1i*t)/(1-q*exp(1i*t))} else
> {kk <- exp(1i*t)*cf(t, m-1, p); g <- (p*kk)/(1-q*kk)}
> g
> }
>
> #===================TESTING================================#
> consecpmf(seq(0,200),3,0.5)
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595