Fit a piecewise regression in matlab and find change point

Question

In matlab, I want to fit a piecewise regression and find where on the x-axis the first change-point occurs. For example, for the following data, the output might be changepoint=20 (I don't actually want to plot it, just want the change point).

data = [1 4 4 3 4 0 0 4 5 4 5 2 5 10 5 1 4 15 4 9 11 16 23 25 24 17 31 42 35 45 49 54 74 69 63 46 35 31 27 15 10 5 10 4 2 4 2 2 3 5 2 2];
x = 1:52;
plot(x,data,'.')

Answer 1

If you have the Signal Processing Toolbox, you can directly use the findchangepts function (see https://www.mathworks.com/help/signal/ref/findchangepts.html for documentation):

data = [1 4 4 3 4 0 0 4 5 4 5 2 5 10 5 1 4 15 4 9 11 16 23 25 24 17 31 42 35 45 49 54 74 69 63 46 35 31 27 15 10 5 10 4 2 4 2 2 3 5 2 2];
x = 1:52;
ipt = findchangepts(data);
x_cp = x(ipt);
data_cp = data(ipt);

plot(x,data,'.',x_cp,data_cp,'o')

The index of the change point in this case is 22.

Plot of data and its change point circled in red:

Answer 2

I know this is an old question but just want to provide some extra thoughts. In Maltab, an alternative implemented by me is a Bayesian changepoint detection algorithm that estimates not just the number and locations of the changepoints but also reports the occurrence probability of changepoints. In its current implementation, it deals with only time-series-like data (aka, 1D sequential data). More info about the tool is available at this FileExchange entry (https://www.mathworks.com/matlabcentral/fileexchange/72515-bayesian-changepoint-detection-time-series-decomposition).

Here is its quick application to your sample data:

% Automatically install the Rbeast or BEAST library to local drive
eval(webread('http://b.link/rbeast')) %

data = [1 4 4 3 4 0 0 4 5 4 5 2 5 10 5 1 4 15 4 9 11 16 23 25 24 17 31 42 35 45 49 54 74 69 63 46 35 31 27 15 10 5 10 4 2 4 2 2 3 5 2 2];

out = beast(data, 'season','none') % season='none': there is no seasonal/periodic variation in the data
printbeast(out)
plotbeast(out)

Below is a summary of the changepoint, given by printbeast():

#####################################################################
#                      Trend  Changepoints                          #
#####################################################################
.-------------------------------------------------------------------.
| Ascii plot of probability distribution for number of chgpts (ncp) |
.-------------------------------------------------------------------.
|Pr(ncp = 0 )=0.000|*                                               |
|Pr(ncp = 1 )=0.000|*                                               |
|Pr(ncp = 2 )=0.000|*                                               |
|Pr(ncp = 3 )=0.859|*********************************************** |
|Pr(ncp = 4 )=0.133|********                                        |
|Pr(ncp = 5 )=0.008|*                                               |
|Pr(ncp = 6 )=0.000|*                                               |
|Pr(ncp = 7 )=0.000|*                                               |
|Pr(ncp = 8 )=0.000|*                                               |
|Pr(ncp = 9 )=0.000|*                                               |
|Pr(ncp = 10)=0.000|*                                               |
.-------------------------------------------------------------------.
|    Summary for number of Trend ChangePoints (tcp)                 |
.-------------------------------------------------------------------.
|ncp_max    = 10   | MaxTrendKnotNum: A parameter you set           |
|ncp_mode   = 3    | Pr(ncp= 3)=0.86: There is a 85.9% probability  |
|                  | that the trend component has  3 changepoint(s).|
|ncp_mean   = 3.15 | Sum{ncp*Pr(ncp)} for ncp = 0,...,10            |
|ncp_pct10  = 3.00 | 10% percentile for number of changepoints      |
|ncp_median = 3.00 | 50% percentile: Median number of changepoints  |
|ncp_pct90  = 4.00 | 90% percentile for number of changepoints      |
.-------------------------------------------------------------------.
| List of probable trend changepoints ranked by probability of      |
| occurrence: Please combine the ncp reported above to determine    |
| which changepoints below are  practically meaningful              |
'-------------------------------------------------------------------'
|tcp#              |time (cp)                  |prob(cpPr)          |
|------------------|---------------------------|--------------------|
|1                 |33.000000                  |1.00000             |
|2                 |42.000000                  |0.98271             |
|3                 |19.000000                  |0.69183             |
|4                 |26.000000                  |0.03950             |
|5                 |11.000000                  |0.02292             |
.-------------------------------------------------------------------.

Here is the graphic output. Three major changepoints are detected (i.e., vertical lines). As a Bayesian method, it tells not just whether there are any changepoints but also the probability of changepoint occurrence as in the Pr(tcp) subplot where the peaks correspond to the detected changepoints.