[R] Additional args to fun in integrate() not found?

[Dimitris Rizopoulos]

you need 'p = p' in the integrate() call; try this version:

gedCDF <- function (yvec, p = 2, mu = 0, scale = 1, numint = 0) {
    k <- sqrt(2)
    # k <- 1
    scale <- scale*k
    zvec <- (yvec - mu) / scale
    cdf <- numeric(length(zvec))
    for (i in seq(zvec)) {
        z <- zvec[i]
        cdf[i] <- if (numint == 0) {
            if (z <= 0) {
                0.5 * (1 - 1/gamma(1/p) * Igamma((1/p), (-z)^p,
                    lower = TRUE))
            } else {
                1 - (0.5 * (1 - 1/gamma(1/p) * Igamma((1/p), (z)^p,
                    lower = TRUE)))
            }
        } else {
            integrate(geddenstandard, lower = -35, upper = z,
                rel.tol = 100*.Machine$double.eps, p = p)$value
        }
    }
    cdf
}

geddenstandard <- function (z, p) {
    p / (2 * gamma(1/p)) * exp(-abs(z)^p)
}

###########

gedCDF(c(1,2,3,4,5), numint=1)


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Sergey Goriatchev" <sergeyg at gmail.com>
To: <r-help at stat.math.ethz.ch>
Sent: Tuesday, February 27, 2007 10:49 AM
Subject: [R] Additional args to fun in integrate() not found?


> Hello, fellow Rdicts,
>
> I have  the code for the program below. I need to integrate a 
> function
> of "x" and "p". I use integrate to integrate over "x" and pass "p" 
> as
> an additional argument. "p" is specified and given default value in
> the argument list. Still, integrate() cannot read "p", unless I
> explicitly insert a numeric value in the integrate() argument list.
> And when I do that, I get the right result, but still some warnings.
>
> Please, help me with these problems:
> 1) why is "p" not recognized?
> 2) what are these warning messages?
>
> PROGRAM CODE:
> ---------------------------
>
> #THIS LIBRARY IS NEEDED FOR THE INCOMPLETE GAMMA FUNCTION
>
> library(zipfR)
>
> #------------------------------------------------------------------------
>
> gedCDF = function(yvec, p=2, mu=0, scale=1, numint=0)
> {
>
> #-------------------------------------------------------------------------
> #Setting k to sqrt(2) and the GED with p=2 coincides with standard 
> normal.
> #Set k=1 and GED with p=1 coincides with Laplace.
>
> k<-sqrt(2)
> #k<-1
> scale<-scale*k
>
> zvec<-(yvec-mu)/scale
> cdf<-matrix(0, length(zvec),1)
>
> for(i in 1:length(zvec))
> {
> z<-zvec[i]
>
> if(numint==0)
> {
> if(z<=0)
> {
> t<-0.5*(1-1/gamma(1/p)*Igamma((1/p),(-z)^p,lower=TRUE))
> }
> else
> {
> t<-1-(0.5*(1-1/gamma(1/p)*Igamma((1/p),(z)^p,lower=TRUE )))
> }
> }
> else
> {
> t<-integrate(geddenstandard, -35, z, subdivisions=1000,
> rel.tol=100*.Machine$double.eps, abs.tol=rel.tol,
> stop.on.error=TRUE, keep.xy=FALSE, aux=NULL,p)
> }
> cdf[i]<-t
> }
> cdf
> }
>
> #-------------------------------------------------------------------------
> geddenstandard = function(z,p)
> {
> f<-p/(2*gamma(1/p))*exp(-abs(z)^p)
> }
> -----------------------------------------------------------------------------------------------------------------------
>
> If I run with this definition I get the following error message and 
> abort:
>> gedCDF(c(1,2,3,4,5), numint=1)
> Error in eval(expr, envir, enclos) : ..1 used in an incorrect 
> context,
> no ... to look in
>
> If I replace "p" in integrate() with 2, I get correct answers, but
> still some warning messages:
>
>> gedCDF(c(1,2,3,4,5), numint=1)
> [[1]]
> [1] 0.8413447
>
> [[2]]
> [1] 0.9772499
>
> [[3]]
> [1] 0.9986501
>
> [[4]]
> [1] 0.9999683
>
> [[5]]
> [1] 0.9999997
>
> Warning messages:
> 1: number of items to replace is not a multiple of replacement 
> length
> 2: number of items to replace is not a multiple of replacement 
> length
> 3: number of items to replace is not a multiple of replacement 
> length
> 4: number of items to replace is not a multiple of replacement 
> length
> 5: number of items to replace is not a multiple of replacement 
> length
>
> -------------------------------
> Do I get these warnings because I define cdf as a matrix and the
> output-cdf is a list?
>
> Please, help me with these!
> Email to my gmail account, please: sergeyg at gmail.com
>
> THanks in advance
> Sergey
>
> ______________________________________________
> 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.
> 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm