# Developers Planet

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

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
``````