I made some experiments and I conclude that if I change the interval (e.g. say c(-100,100)!!!) then the minimum value changes in respect to the simulated data.
So, probably the problem lies in the "optimal" values of the interval? On Sat, Jul 26, 2008 at 1:00 PM, Fotis Papailias <[EMAIL PROTECTED]>wrote: > Hi, > > thanks for your message. > > You mean to rewrite the function like that: > > lw <- function(d, x, im) > { > peri1 <- per(x) > len <- length(x) > m <- len/im > peri <- peri1[2:(m+1)] > z <- c(1:m) > freq <- ((2*pi)/len) * z > result <- log(sum(freq^(2*d-1)*peri))-(2*d)/m * sum(log(freq)) > } > > > and then to use: > > target <- function(d) lw(x, d, im=2) > k <- optimize(target, interval=c(0, 0.5)) > > or just, > > k <- optimize(lw, x, im) > > But I still get the same values > > $minimum > [1] 0.4999542 > > $objective > [1] 2.739509 > > after the optimization, no matter my "x" series. It seems that I am still > doing something wrong in the optimize() or in the result() function inside > the loop. But comparing it to the original S-PLUS code, it seems like it is > fine. > > Send me if you have any other ideas please. > > Thanks again. > > fotis > > > On Sat, Jul 26, 2008 at 12:48 PM, Duncan Murdoch <[EMAIL PROTECTED]>wrote: > >> On 26/07/2008 7:40 AM, Fotis Papailias wrote: >> >>> Dear R-users, >>> >>> I have sent another mail some hour ago about a matlab Code I was trying >>> to >>> translate in R. >>> >>> Actually I have found a simpler code originally written in S-PLUS for the >>> same function. >>> Author's page -> http://math.bu.edu/people/murad/methods/locwhitt/ >>> >>> ============================================================= >>> >>> rfunc_function(h, len, im, peri) >>> # h -- Starting H value for minimization. >>> # len -- Length of time series. >>> # im -- Use only len/im frequencies. >>> # peri -- Periodogram of data. >>> { >>> m <- len %/% im >>> peri <- peri[2:(m + 1)] >>> z <- c(1:m) >>> freq <- (2 * pi)/len * z >>> result <- log(sum(freq^(2 * h - 1) * peri)) - (2 * h)/m * >>> sum(log(freq) >>> ) # cat("H = ", h, "R = ", result, "\n") >>> drop(result) >>> } >>> >>> >>> locwhitt_function(data, h = 0.5, im = 2) >>> # data -- Time series. >>> # h -- Starting H value for minimization. >>> # im -- Use only N/im frequencies where N is length of series. >>> >>> { >>> peri <- per(data) >>> len <- length(data) >>> return(nlminb(start = h, obj = rfunc, len = len, im = im, peri = >>> peri)$ >>> parameters) >>> } >>> =============================================================== >>> >>> The author who has written the above S-PLUS code uses two functions (with >>> the locwhitt_function he lets the user to define the data and the >>> parameters >>> and with the rfunc_function he does the minimization.) >>> >>> Mine translation is in R is: >>> >>> where I use a joint function compared to the the above author >>> >>> >>> ================================================================ >>> >>> lw <- function(x, d, im) >>> { >>> peri1 <- per(x) >>> len <- length(x) >>> m <- len/im >>> peri <- peri1[2:(m+1)] >>> z <- c(1:m) >>> freq <- ((2*pi)/len) * z >>> result <- log(sum(freq^(2*d-1)*peri))-(2*d)/m * sum(log(freq)) >>> } >>> >>> ================================================================= >>> >>> which seems to run ok. >>> >>> But when I do >>> >>> k <- optimize(lw, x, im=2, interval=c(0, 0.5)) >>> >>> I always get the same result no matter the (simulated) data in x! >>> >>> The parameter of interest to be minimized is "d". Does anyone know how to >>> edit the function "optimize" so it can work properly? >>> >> >> optimize() is fine, but the way you're calling it is not. It optimizes a >> function over the first argument. So you could rewrite lw to put d first, >> or write a new function which calls it, e.g. >> >> target <- function(d) lw(x, d, im) >> >> and then >> >> optimize(target, interval=c(0, 0.5)) >> >> Because target is defined in the global environment, it will look there >> for x and im, and you don't need to pass them as arguments: unless x and im >> aren't defined there too! >> >> Duncan Murdoch >> > > > > -- > fp > -- fp [[alternative HTML version deleted]] ______________________________________________ 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.