How can I chose the Threshold value, in change point detection for a set of data with unknown number of changepoints

A = [zeros(1,100) 1:100 99:-1:50 50*ones(1,250)] + 10*rand(1,500);

[TF,S1,S2] = ischange(A,'linear','Threshold',200);

segline = S1.*(1:500) + S2;

plot(1:500,A,1:500,segline)

legend('Data','Linear Regime')

this is a code given in MATTLAB Help center,Ducumentaion

https://uk.mathworks.com/help/matlab/ref/ischange.html?lang=en

Answer 1

[Kaiguang]

If anybody wants to try a Bayesian alternative for changepoint detection, one possbility is a Bayesian method called BEAST (Bayesian Estimator of Abrupt Change, Seasonality and Trend) available from this FileExchange entry or https://github.com/zhaokg/Rbeast. It can be instally installed by running eval(webread('http://b.link/beast',weboptions('cert',''))). Below is a quick illustration using the sample data provided above to illustrate its basic use:

% Quick installation of BEAST to a temporary path on your local drive
eval(webread('http://b.link/beast',weboptions('cert','')))  

% The simualted time series from the original question
y = [zeros(1,100) 1:100 99:-1:50  50*ones(1,250)] + 10*rand(1,500); 

% Apply beast to y. Here season='none' indicates that y has no periodic/seasonal component
o =  beast(y, 'season','none')    

printbeast(o)
plotbeast(o)

Below is a summary of the number and locations of the changepoints detected:

#####################################################################
#                      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.914|*********************************************** |
|Pr(ncp = 4 )=0.083|*****                                           |
|Pr(ncp = 5 )=0.002|*                                               |
|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.91: There is a 91.4% probability  |
|                  | that the trend component has  3 changepoint(s).|
|ncp_mean   = 3.09 | 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  = 3.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                 |199.000000                 |1.00000             |
|2                 |252.000000                 |0.92867             |
|3                 |96.000000                  |0.89042             |
|4                 |471.000000                 |0.01800             |
|5                 |413.000000                 |0.00733             |
|6                 |435.000000                 |0.00692             |
|7                 |483.000000                 |0.00679             |
|8                 |448.000000                 |0.00579             |
|9                 |343.000000                 |0.00204             |
|10                |63.000000                  |0.00154             |
.-------------------------------------------------------------------.

Here is the plot. The trend is fitted using a piecewise polynomial model. Again, as a Bayesian method, BEAST assumes the order of the polynomial as uknowns. The orders of the polynomial needed to adequately fit the trend are estimated over time, as depicted iin the tOrder subplot. The 1st and 4th segments are flat lines, so their estimated poly orders are close to zeros.

Answer 2

[Pranav Verma]

Hi Mohamad,

From your question I understand that you are trying change point detection for the data sample. You can try using the findchangepts function in MATLAB. In this function, there exists a Name,Value pair for setting the minimum threshold : MinThreshold.

In the documentation of the findchangepts, ity is explained how to decide upon the threshold value:

• The proportionality constant, denoted by β and specified in 'MinThreshold', corresponds to a fixed penalty added for each changepoint. findchangepts rejects adding additional changepoints if the decrease in residual error does not meet the threshold. Set 'MinThreshold' to zero to return all possible changes.
• If you do not know what threshold to use or have a rough idea of the number of changepoints in the signal, specify 'MaxNumChanges' instead. This option gradually increases the threshold until the function finds fewer changes than the specified value.

You can find the mathematical analysis of changepoint detection and the above mentioned points in the documentation.

Also you can refer to the below thread explaning findchangepts:

https://www.mathworks.com/matlabcentral/answers/621918-changepoint-analysis-findchangepts-how-does-it-work?s_tid=answers_rc1-2_p2_MLT

Thanks

[Mohamad Shukfa]

Thanks