So guys,
I am writing to ask for an opinion from you...
As I begin to write in the previous post, I developed a function
called myintegrate that is based in integrate but my function give
better results for the problem I am dealing. My development was
inspired in this post
http://r.789695.n4.nabble.com/integration-td843252.html, by Ravi
Varadhan.
To check if the results are good I compare with the ones from an
analytical version of the problem.
This is the new code,
dsad <- Vectorize(FUN function(y, frac, rate, sad, samp="Poisson",
trunc=0, ...){
f0 <- function(y,frac,n) {
f1 <- function(y,frac,n){
dpois(y,frac*n)
}
dcom <- paste("d",deparse(substitute(sad)),sep="")
dots <- c(as.name("n"),list(...))
f2 <- call(dcom,dots)
f <- function(n){
f1(y,frac,n)*f2(n,rate)
}
myintegrate <- function() {
r <- 0
r1 <- 1
x1 <- 0
dx <- 20
while(r1 > 10e-500){
r1 <- integrate(f,x1,x1+dx)$value
r <- r + r1
x1 <- x1 + dx
}
r + integrate(f,x1,Inf)$valu
}
myintegrate()
}
f0(y,frac,n)/(1-f0(trunc,frac,n))
},"y")
dsad(200, 0.1, 0.1, exp)
whose result is 6.223015e-61, correct.
And the below is the old one,
dsad <- Vectorize(FUN function(y, frac, rate, sad, samp="Poisson",
trunc=0, ...){
f0 <- function(y,frac,n) {
f1 <- function(y,frac,n){
dpois(y,frac*n)
}
dcom <- paste("d",deparse(substitute(sad)),sep="")
dots <- c(as.name("n"),list(...))
f2 <- call(dcom,dots)
f <- function(n){
f1(y,frac,n)*f2(n,rate)
}
myintegrate <- function() {
r <- 0
r1 <- 1
x1 <- 0
integrate(f,x1,Inf)$valu
}
myintegrate()
}
f0(y,frac,n)/(1-f0(trunc,frac,n))
},"y")
dsad(200, 0.1, 0.1, exp)
whose result is 6.554036e-80, incorrect.
So, what do you think of it?
I will really apprettiate any comments.
Thanks in advance.
On Mon, Aug 29, 2011 at 4:06 AM, . . <xkziloj at gmail.com>
wrote:> Hi all,
>
> I am utilizing integrate() to get the area below a function that shape
> like an exponential, getting very low for higher values of the
> horizontal coordinate y.
>
> In the code below, I present the tricky way I can improve the result
> obtained. I compare this result with the one I get from Mathematica
> and also with the exact solution I have for this particular case.
>
> /*----------
> func <- function(y, a, rate){
> ?x <- function(n,rate) {
> ? ?rate*exp(-n*rate)
> ?}
> ?boi <- function(y,n,a){
> ? ? ? ? ? ? ? ? ? ? ?w <- y*log(a*n)-lfactorial(y)-a*n
> ? ? ? ? ? ? ? ? ? ? ?exp(w)
> ? ? ? ? ? ? ? ? ? ?}
> ?f <- function(n){
> ? ?boi(y,n,a)*x(n,rate)
> ?}
> ?r <- 0
> ? ? ? ? ? ? ? ? ? ?r1 <- 1
> ? ? ? ? ? ? ? ? ? ?x1 <- 0
> ? ? ? ? ? ? ? ? ? ?dx <- 20
> ? ? ? ? ? ? ? ? ? ?while(r1 > 10e-1000){
> ? ? ? ? ? ? ? ? ? ? ?r1 <- integrate(f,x1,x1+dx)$value
> ? ? ? ? ? ? ? ? ? ? ?r <- r + r1
> ? ? ? ? ? ? ? ? ? ? ?x1 <- x1 + dx
> ? ? ? ? ? ? ? ? ? ?}
> ? ? ? ? ? ? ? ? ? ?r + integrate(f,x1,Inf)$valu
> }
> func(200,0.1,0.1)
> ----------*/
>
> Altought I get better results, the value of dx must be carefully
> selected. So I ask, there is another method that can give me better
> results?
>