# How to find the best piecewise-smooth fit

I have data points $(\theta,y)$ that when plotted in a x-y (with $\theta$ being the x) graph can be fit by the curve

$$y=-\alpha-\frac{\beta}{\cos(\theta_{0}-\theta)}$$

i.e. I can find $\alpha, \beta,\theta_{0}$ that make the above curve fit my data via least squares.

There are really three ranges of $\theta$ in each of which I can do a fit with the above equation (each range will have different $\alpha,\beta,\theta_{0}$ parameters).

Is there a formalism that provides the boundaries of these ranges? I.e. that gives me the angle $\theta$ that is the boundary of the ranges?

Perhaps there is a way of minimizing a cost function that combines the residuals of all 3 ranges? How to define such a cost function?

• @Glen_b I was thinking that perhaps there is a formalism that would stitch together the 3 fits guaranteeing some smoothness among the 3 fits rather than just discontinuities. Jan 7, 2017 at 18:35
• Such as a piecewise-smooth but everywhere continuous fit, say, which would somewhat akin to a piecewise linear model? Jan 7, 2017 at 20:28
• Yes I think you are describing what I have in mind. I manually know that there are roughly 3 regions in each of which I can fit the curve with the model above. But doing it manually is somewhat imprecise/discontinuous. It is possible that a piecewise spline regression would be easier but what I would like is to fit the data with the model given above in a piecewise smooth way if possible. Jan 7, 2017 at 22:02

It looks like what you are looking for is a threshold model, also called segmented regression, or change-point model (also the latter refers more to time series settings). In a threshold model, parameters change depending on the threshold variable (I guess your $\theta$) are within some ranges.

### Estimation

To simplify, let us denote by $x$ the covariate that is also threshold variable, by $x_0$ the transition/threshold point. To estimate your mode, you can do a conditional/concentrated estimation: the idea is that, conditional on a given $x_0$, you can use the standard procedure to estimate your parameters $\alpha_0(x_0)$, $\beta_0(x_0)$, and obtain a cost function (sum of squares, likelihood...), let's write it $SSR(x_0)$. Now you simply look for the $x_0$ that minimizes that cost function. Given that the regime indicator is a discontinuous function, you can not use standard algorithms. A simple solution is to run a grid search over possible values. As you are looking for multiple threshold, you will have a 2, 3-dimensional grid. To avoid the dimensionality problem, you could run a conditional grid search, looking first for only one threshold, there for a second given the first one, eventually iterating, etc.

In steps:

1. Write your cost function, given regression parameters $\alpha^1$, $\beta^1$,$\theta_0^1$, $\alpha^2$, $\beta^2$,$\theta_0^2$ and threshold value $x_0$.

2. Write a function optimizing step 1, with a given threshold value. I.e. find $\hat{\alpha}^1(x_0)$, $\hat{\beta}^1(x_0)$,$\hat{\theta}_0^1(x_0)$, $\hat{\alpha}^2(x_0)$, $\hat{\beta}^2(x_0)$,$\hat{\theta}_0^2(x_0)$, given $x_0$.

3. Select a range of candidates for $x_0$. Typically, one sort x, take unique values, and keep values from the 10% to 90% quantiles. You can use less than the full values obviously. use however existing ones (i.e., if only 2.1 and 2.2 are observed, using 2.01, 2.02 etc as candidates will give same result).

4. Run the function of step 2 for every possible candidates

5. Pick the x value that gives smallest cost function in step 4. This is your estimate of $\hat{x}_0$

6. Use the function in 2, given now the optimal $\hat{x}_0$. This leads the optimal parameters. $\hat{\alpha}^1(\hat{x}_0)$, $\hat{\beta}^1(\hat{x}_0)$, etc.

Find below a small code as example, in R:

n <- 500
x <- sort(rnorm(n, mean=5, sd=2))
fo <- function(x, alpha, beta, theta) alpha -beta / cos(theta-x)
df <- data.frame(x=x, y=ifelse(x<5, fo(x, 1, 3, 5), fo(x, 2, 4,     5))+rnorm(n, sd=10))
plot(y~x, data=df, type="l")

## return sum of squares: simple and with regimes (2 here)
nls_ssr <- function(par, x, y) sum( (y-fo(x, par, par, par))^2 )

nls_ssr_2reg <- function(par, x, y, thresh) {
ssr_reg1 <- nls_ssr(par[1:3],x[x<thresh], y[x<thresh])
ssr_reg2 <- nls_ssr(par[4:6],x[x>=thresh], y[x>=thresh])
ssr_reg1+ssr_reg2
}

## grid search for all regimes
x_search <- x[ (n*0.1):(n*0.9)]
ssrs <- sapply(x_search, function(thresh) optim(c(2,3,5, 2, 3, 5),     nls_ssr_2reg, y=df$y, x=df$x, thresh=thresh)$value) #3 select best, re-estimate x_search[which.min(ssrs)] optim(c(2,3,5, 2, 3, 5), nls_ssr_2reg, y=df$y, x=df\$x, thresh=x_search[which.min(ssrs)])


This gives regression parameters values of 1.719521 3.147032 5.000069 1.136988 4.053732 4.999990, and threshold value 6.5, not too far from the true values 1,3,5,2,4,5 and 5

• How does your solution account for the continuity constraints? (to be "piecewise smooth") Jan 8, 2017 at 19:21