Thanks for your reply Michael, it seems I have a lot of things to learn yet but for sure, your response is being very helpful in this proccess. I will try to explore every point you said:
A doubt I have is, if I define "func <- function(x,y) x + y" how can I integrate it only in "x"? My solution for this would be to define "func <- function(x) x + y". Is not ok? Also, with respect to the helper functions I'd created, I am wondering if you can see a better organization for my code. It is so because this is the only way I can see. Particularly I do not like how I am using "results", but I can not think in another form. Thanks in advance. On Thu, Sep 1, 2011 at 2:44 PM, R. Michael Weylandt <michael.weyla...@gmail.com> wrote: > 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. >> >> > >> >> > >> > >> > > > ______________________________________________ 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.
