KCB February 2016

Fitting boxcar (rectangle) function in R

I am trying to fit a boxcar / rectangle function to a dataset in R. I am using nls and a custom function describing a unit pulse of varying width. Here is what I have so far:

# Defines the unit pulse function
# The left side of the pulse is at x0, the right at x1
pulse <- function(x, x0, x1) {
  if (x >= x0 & x <= x1) {
    return (1)
  } else {
    return (0)
  }
}

xdata <- 1:30
ydata <- c(sample(-1:1, 10, replace = TRUE), sample(10:15, 10, replace = TRUE), sample(-1:1, 10, replace = TRUE))

plot(xdata, ydata)

df <- data.frame(xdata, ydata)

fitfit <- nls(ydata ~ I(A * pulse(xdata, L, R) + B), df, start = list(L = 0, R = 1, B = 0, A = 10))

I am having trouble understanding the error I get:

Error in qr(.swts * attr(rhs, "gradient")) : dims [product 4] do not match the length of object [30] In addition: Warning messages: 1: In if (x >= x0 & x <= x1) { : the condition has length > 1 and only the first element will be used 2: In if (x >= x0 & x <= x1) { : the condition has length > 1 and only the first element will be used 3: In if (x >= x0 & x <= x1) { : the condition has length > 1 and only the first element will be used 4: In if (x >= x0 & x <= x1) { : the condition has length > 1 and only the first element will be used 5: In if (x >= x0 & x <= x1) { : the condition has length > 1 and only the first element will be used 6: In if (x >= x0 & x <= x1) { : the condition has length > 1 and only the first element will be used 7: In .swts * attr(rhs, "gradient") : longer object length is not a multiple of shorter object length

Answers


HubertL February 2016

your pulse function doesn't return a vector, try this one:

pulse <- function(x, x0, x1) {
     ifelse (x >= x0 & x <= x1,1,0)
}


G. Grothendieck February 2016

Non-smooth functions can cause problems. Try brute force over the two boundaries of the pulse and optimization over the linear parameters. Note that below we have added a set.seed to make ydata reproducible.

library(nls2)

set.seed(123)
xdata <- 1:30
ydata <- c(sample(-1:1, 10, replace = TRUE), sample(10:15, 10, replace = TRUE), 
  sample(-1:1, 10, replace = TRUE))
df <- data.frame(xdata, ydata)

pulse <- function(x, x0, x1) (x >= x0 & x <= x1) + 0
st <- subset(expand.grid(L = xdata, R = xdata), L < R)

nls2(ydata ~ cbind(pulse(xdata, L, R), 1), df, start = st, alg = "plinear-brute")

giving some error messages which you can ignore and finally this output where .lin1 and .lin2 correspond to A and B in the question:

Nonlinear regression model
  model: ydata ~ cbind(pulse(xdata, L, R), 1)
   data: df
    L     R .lin1 .lin2 
 11.0  20.0  12.4   0.2 
 residual sum-of-squares: 49.6

Number of iterations to convergence: 435 
Achieved convergence tolerance: NA

Post Status

Asked in February 2016
Viewed 3,355 times
Voted 12
Answered 2 times

Search




Leave an answer