如果您有信号处理工具箱，则可以直接使用findchangepts函数（有关文档，请参阅https://www.mathworks.com/help/signal/ref/findchangepts.html）：

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')

字符串
在这种情况下，变化点的索引是22。
数据图及其变化点以红色圈出：

的数据

我知道这是一个老问题，但我只是想提供一些额外的想法。在Maltab中，我实现的一个替代方案是贝叶斯变点检测算法，该算法不仅估计变点的数量和位置，还报告变点的发生概率。在当前的实现中，它只处理类似时间序列的数据（也称为1D顺序数据）。有关该工具的更多信息，请参阅此FileExchange条目（https://www.mathworks.com/matlabcentral/fileexchange/72515-bayesian-changepoint-detection-time-series-decomposition）。
下面是它在示例数据中的快速应用：

% 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)

字符串
下面是由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             |
.-------------------------------------------------------------------.

型
以下是图形输出。检测到三个主要的变化点：

的数据

你可以使用sgolayfilt函数，也就是多项式拟合数据，或者重现OLS方法：http://www.utdallas.edu/~herve/Abdi-LeastSquares06-pretty.pdf（有 a+bx 符号而不是 ax+b）
对于 ax+B 的线性拟合：
如果用长度为2n+1的常数向量替换x：[-n，…0…n]在每一步中，您将获得以下用于滑动回归系数的代码：

for i=1+n:length(y)-n
    yi = y(i-n : i+n); 
    sum_xy = sum(yi.*x);

    a(i) = sum_xy/sum_x2;
    b(i) = sum(yi)/n;
end

字符串
请注意，在此代码中，B表示数据的滑动平均值，a是最小二乘斜率估计值（一阶导数）。