[R] S-PLUS code in R

[Duncan Murdoch]

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