## I see two types of difficulties, one type is actual errors in
## usage, and the other are what I consider less than optimal
## programming style.
##
## Errors
##
## a. The primary errors are indeed in scoping, as you guessed.
## Functions look in the global space for names that they do not find
## inside the current function. Thus the variable references to j, k,
## d inside the integrand function do not look inside the test
## function for a possible location. Therefore integrand does not see
## the values in your j and k loops. The way around that is to give
## integrand a second argument and then use it. See the revision
## below which fixes all(?) the scoping errors. Similarly, the test
## function may not see the value of "n".
##
##
## There are five stylistic issues that contribute to the difficulty in spotting
## the problem.
##
## b. I recommend always using " <- " for assignment and never using
## "_". For comparison, beginning with R 1.9.1, R no longer honors the
## underscore for assignment.
##
## I think the definition of test with "_" and of integrand with " <- "
## interfered with recognizing that this is a nested set of function
definitions.
##
## c. Each level of nesting is not indented within the previous level.
## This makes it hard to see the nesting pattern.
##
## d. I recommend spaces after all commas. This makes it easier for people
## to read the code.
##
## e. The integrand function is identical throughout the double loop.
## Therefore I took it out of the loop. This has two effects. One is
## that the program will be much more efficient since it only has to
## define the function once. The other is that it is easier to read,
## hence debug, since the constants are not inside the loop.
##
## f. There is potential danger in using the name "c" for a variable.
## "c" is the name of the catenate function which is used inside every
## S-Plus function. You run the risk of masking the function and generating
## very untraceable errors. I changed its name in your function to "cc".
## An additional note, which is probably a consequence of the simplification,
## is that the result of the integration is never used or returned.
## With all this, there is something still wrong. I get the message
## Problem in .Fortran("divset",: subroutine divset:
## Argument 2 has zero length
## I believe the problem is detected inside nlminb, but I can't find it.
## Rich
n <- 5
result <- diag(0,n)
test <- function(b, n) {
on.exit(browser())
cc <- b[1]
d <- b[2]
integrand <- function(x, j, k, d) {
sin(x/2)^2/(2*pi*j-x)/(2*pi*k+x)*abs(x/(2*pi*j))^(-2*d)
}
for(j in 1:n) {
for (k in 1:n) {
integrate(integrand,
subdivisions=1000,
lower=-Inf, upper=Inf,
j=j, k=k, d=d
)$integral
}
}
return(cc-exp(d))
}
testout <- nlminb(start=c(0.1,1.0),
obj=test,
lower=c(0.05,0.5),
upper=c(0.45,5.0),
n=n)$parameters
testout
|