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)
>
> ______________________________________________
> [email protected] 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, [EMAIL PROTECTED]
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
______________________________________________
[email protected] 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.
