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 [EMAIL PROTECTED]
To: r-help@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.683
[[5]]
[1] 0.997
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: [EMAIL PROTECTED]
THanks in advance
Sergey
__
R-help@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
__
R-help@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.
