Leaving aside some other issues that this whole email chain has opened up, I'd guess that your most immediate problem is that you are trying to numerically integrate the PMF of a discrete distribution but you are treating it as a continuous distribution. If you took the time to properly debug (as you were instructed yesterday) you'd probably find that whenever you call dpois(x, lambda) for x not an integer you get a warning message.
Specifically, check this out > integrate(dpois,0,Inf,1) 9.429158e-13 with absolute error < 1.7e-12 > n = 0:1000; sum(dpois(n,1)) 1 I could be entirely off base here, but I'm guessing that many of your problems derive from this. On another basis, please, please read this: http://google-styleguide.googlecode.com/svn/trunk/google-r-style.html or this http://had.co.nz/stat405/resources/r-style-guide.html And, perhaps most importantly, don't rely on the black magic of values moving in and out of functions (lexical scoping). Seriously, just don't do it. If you have helper functions that need values, actively pass them: you will save yourself hours of trouble when (not if) you debug your functions. I'm looking, for example, at g() in the first big block of code you provided. Call it g(a,n) and spend the extra 4 keystrokes to pass the values. It makes everyone happier. Michael On Thu, Sep 1, 2011 at 12:37 PM, . . <xkzi...@gmail.com> wrote: > So, please excuse me Michael, you are completely sure. I will try > describe I am trying to do, please let me know if I can provide more > info. > > The idea is provide to "func" two probability density functions(PDFs) > and obtain another PDF that is a compound of them. In a final analysis > this characterize an abundance distribution for me. The two PDFs are > provided through "f" and "g" and there is some manipulation here > because I need flexibility to easily change this two funcions. > > In the code provided, "f" is the Exponential distribution and "g" is > the Poisson distribution. For this case, I have the analytical > solution, below. This way I can check the result. But I am also > considering other combinations of "f" and "g" that have difficult, or > even does not have analitical solution. This is the reason why I am > trying to develop "func". > > func2 <- function(y, frac, rate, trunc=0, log=FALSE) { > is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) > abs(x - round(x)) < tol > if(FALSE %in% sapply(y,is.wholenumber)) > print("y must be integer because dpoix is a discrete PDF.") > else { > f <- function(y){ > b <- y*log(frac) > m <- log(rate) > n <- (y+1)*log(rate+frac) > if(log)b+m-n else exp(b+m-n) > } > f(y)/(1-f(trunc)) > } > } > > func2(200,0.05,0.001) > [1] 0.000381062 > > In theory, the interval of integration is 0 to Inf, but for some tests > I did, go up to 2000 may still provide reasonable results. > > Also, as it seems, I am still writing my first functions in R and > suggestions are welcome, please. > > Again, appologies for my previous mistake. It was not my intention to > blame about "integrate". > > On Thu, Sep 1, 2011 at 11:49 AM, R. Michael Weylandt > <michael.weyla...@gmail.com> wrote: > > I'm going to try to put this nicely: > > > > What you provided is not a problem with integrate. Instead, you provided > a > > rather unintelligible and badly-written piece of code that (miraculously) > > seems to work, though it's not well documented so I have no idea if > 1.3e-21 > > is what you want to get. > > > > Let's try this again: per your original request, what is the problem with > > integrate? > > > > If instead you feel there's something wrong with your code, might I > suggest > > you just say that and ask for help, rather than passing the blame onto a > > perfectly useful base function. > > > > Oh, and since you asked that I propose something: comment your code. > > > > Michael > > > > On Thu, Sep 1, 2011 at 10:33 AM, . . <xkzi...@gmail.com> wrote: > >> > >> Hi Michael, > >> > >> This is the problem: > >> > >> func <- Vectorize(function(x, a, sad, samp="pois", trunc=0, ...) { > >> result <- function(x) { > >> f1 <- function(n) { > >> f <- function() { > >> dcom <- paste("d", sad, sep="") > >> dots <- c(as.name("n"), list(...)) > >> do.call(dcom, dots) > >> } > >> g <- function() { > >> dcom <- paste("d", samp, sep="") > >> lambda <- a * n > >> dots <- c(as.name("x"), as.name("lambda")) > >> do.call(dcom, dots) > >> } > >> f() * g() > >> } > >> integrate(f1,0,2000)$value > >> # adaptIntegrate(f1,0,2000)$integral > >> > >> # n <- 0:2000 > >> # trapz(n,f1(n)) > >> > >> # area(f1, 0, 2000, limit=10000, eps=1e-100) > >> } > >> return(result(x) / (1 - result(trunc))) > >> }, "x") > >> func(200, 0.05, "exp", rate=0.001) > >> > >> If you could propose something I will be gratefull. > >> > >> Thanks in advance. > >> > >> On Thu, Sep 1, 2011 at 10:55 AM, R. Michael Weylandt > >> <michael.weyla...@gmail.com> wrote: > >> > Mr ". .", > >> > > >> > MASS::area comes to mind but it may be more helpful if you could say > >> > what > >> > you are looking for / why integrate is not appropriate it is for > >> > whatever > >> > you are doing. > >> > > >> > Strictly speaking, I suppose there are all sorts of "alternatives" to > >> > integrate() if you are willing to be really creative and build > something > >> > from scratch: diff(), cumsum(), lm(), hist(), t(), c(), .... > >> > > >> > Michael Weylandt > >> > > >> > On Thu, Sep 1, 2011 at 9:53 AM, B77S <bps0...@auburn.edu> wrote: > >> >> > >> >> package "caTools" > >> >> see ?trapz > >> >> > >> >> > >> >> . wrote: > >> >> > > >> >> > Hi all, > >> >> > > >> >> > is there any alternative to the function integrate? > >> >> > > >> >> > Any comments are welcome. > >> >> > > >> >> > Thanks in advance. > >> >> > > >> >> > ______________________________________________ > >> >> > R-help@r-project.org 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. > >> >> > > >> >> > >> >> -- > >> >> View this message in context: > >> >> > >> >> > http://r.789695.n4.nabble.com/Alternatives-to-integrate-tp3783624p3783645.html > >> >> Sent from the R help mailing list archive at Nabble.com. > >> >> > >> >> ______________________________________________ > >> >> R-help@r-project.org 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. > >> > > >> > > > > > > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org 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.
