Change Point Detection in the Frequency Domain with Statistical Reliability
Authors: Akifumi Yamada
Tomohiro Shiraishi
Shuichi Nishino
Teruyuki Katsuoka
Kouichi Taji
Ichiro Takeuchi
Venue: Accepted by TMLR ()

Abstract: Effective condition monitoring in complex systems requires identifying change points (CPs) in the frequency domain, as the structural changes often arise across multiple frequencies. This paper extends recent advancements in statistically significant CP detection, based on Selective Inference (SI), to the frequency domain. The proposed SI method quantifies the statistical significance of detected CPs in the frequency domain using $p$-values, ensuring that the detected changes reflect genuine structural shifts in the target system. We address two major technical challenges to achieve this. First, we extend the existing SI framework to the frequency domain by appropriately utilizing the properties of discrete Fourier transform (DFT). Second, we develop an SI method that provides valid $p$-values for CPs where changes occur across multiple frequencies. Experimental results demonstrate that the proposed method reliably identifies genuine CPs with strong statistical guarantees, enabling more accurate root-cause analysis in the frequency domain of complex systems.

Published in Transactions on Machine Learning Research (06/2025)
Change Point Detection in the Frequency Domain with Statistical Reliability
Akifumi Yamada
Nagoya University
Tomohiro Shiraishi
shiraishi.
Nagoya University
RIKEN
Shuichi Nishino
nishino
Nagoya University
RIKEN
Teruyuki Katsuoka
Nagoya University
Kouichi Taji
Laji@naggyya.u.pp
Nagoya University
Ichiro Takcuchi
takeachi, takeuchi. ragge wwgoyaau.na.............
Nagoya University
RIKEN
Reviewed on OpenReview: https: // openreview. net/ forum? id- FWWdaHH23gN
Abstract
Effective condition monitoring in complex systems requires identifying change points (CPs) in the frequency domain, as the structural changes often arise across multiple frequencies. This paper extends recent advancements in statistically significant. CP detection, based on Selective Inference (SI), the frequency domain. The proposed SI method quantifies the statistical significance of detected CPs in the frequency domain using p-values, ensuring that the detected changes reflect genuine structural shifts in the target system. We address two major technical challenges to achieve this. First, we extend the existing SI framework to the frequency domain by appropriately utilizing the properties of discrete Fourier transform (DFT). Second, we develop an SI method that provides valid p-values for CPs where changes occur across multiple frequencies. Experimental results demonstrate that the proposed method reliably identifies genuine CP's with strong statistical guarantees, enabling more accurate root-cause analysis in the frequency domain of complex systems.
1 Introduction
To detect failures in complex systems, it is crucial accurately identify the frequencies where the significant changes occur. By pinpointing these specific frequencies, the root causes of faults in the systems can be identified. In this paper, we address the change point (CP) detection problem in the frequency domain and propose a method that can identify slatistically significand changes in specific frequencies. Quantifying statistical significance CP detection ensures that the detected changes in the frequency domain reflect genuine structural alterations rather than random noise. Mcasures such as p-values provide a quantitative framework differentiate real changes from spurious ones, effectively mitigating the risk of incorrect decisions in data-driven systems.

---

Published in Transactions on Machine Learning Research (06/2025)
This study was motivated by recent developments in statistically significant CP detection based Selective Inference (ST) (Taylor & Tibshirani, 2015; Fithian et al., 2015; Lee & Taylor, 2014). Prior to ST, quantifying the significance of detected CPss was challenging due to the issue of double dipping, i.e. using the same data to both identify and test CPs inflates false positive findings. Traditional statistical approaches prior to SI primarily locused on lesting only whether a CP exists or not within a certain range using multiple testing framework based on asymptotic distribution (Page, 1954; Mika et al., 1999). Therefore, providing statistical significance measures, such as p-values, for specific points or frequencies was not feasible. SI is novel statistical inference framework designed for data-driven hypotheses. In the context of CP detection, it enables the evaluation for the statistical significance of detected CPs conditional on a specific CP detection algorithm, thereby resolving the aforementioned double dipping issue. The use of S1 for statistically significant CP detection in the time domain has been actively studied recently (Hyun et al., 2018; Duy et al., 2020). Our contribution in this study is to extend these methods and realize statistically significant CP detection in the froquency domain.
One of the diffieulties in frequeney-domain CP detection is that changes often appear across multiple fre- quencies due to shared underlying phenomena influencing broad signal characteristics. Figure 1 shows an example of froquency-domatt CP detection problems along with the statistically significant CPs identified using our proposed SI method. Panel (a) shows a time series transformed into signals across five frequencies. Panel (b) presents the CP detection results for each individual frequency d E (do, d1, d2, d3,d4) 1 Panel (c) integrates CPs across multiple frequencies, detecting three changes A, B, and C (at CP A, changes occur in frequencies d1 and d2; at CP B, d2 and d3; and at CP C, only in do). In table (d), the obtained p-values for changes A, B, and C using the proposed method are 0.008, 0.006, and 0.752 (denoted selective p-values); respectively, indicating that at a significance level of 0.05 = 0.0167 decided by Bonferroni correction, changes A and B are stalistically significant CPs2.
3
To perform statistically significant CP detection, as illustrated in Figure 1, two key technical challenges must be addressed. The first challenge is to extend the SI framework to the frequency domain. We tackle this by formulating the test statistic and the conditioning on the hypothesis selection that accurately account for the properties of discrete Fourier transform (DFT). The second challenge involves quantifying the statistical significance of CPs shared across multiple frequencies. Addressing this requires solving combinatorial optimization problem, which is computationally infeasible to solve globally optimally, necessitating heuristic approach to obtain the approximate solution. In this study, we employ method based on simulated annealing (SA) (Kirkpatrick et al., 1983; Černý 1985) and implement the SI framework to appropriately quantify the statistical significance of the approximate solution derived from the heuristic algorithm.
The proposed CP detection method consists of two stages. In the first stage, CP candidates in the frequency domain are selected. Since the selection of these CP candidates is formulated as a combinatorial optimization problem, a hcuristic algorithm is used to derive significance of each CP candidate selected in the first stage is quantified in the form of p-values using the SI framework. Among the CP candidates, only those with p-values below significance level (e.g., 0.05 or 0.01) are eventually detected as final CPs. The probability of the final detected CPs being false positives is theoretically guaranteed to be below the specified significance level.
approximate solution. In the second stage, the statistical
The rest of this paper is organizod as follows. Section formulates the problem. Section 3 explains the heuristic algorithm used to select. CP candidates in the first stage of the proposed method. Section 4 de- scribes the method for quantifying the statistical significance of CP candidates using the ST framework in the second stage. In Section 5, we demonstrate the effectiveness of our proposed method through compre- hensive numerical experiments using both synthetic and real-world data. For reproducibility, our implemen- tation is available at https://githu/ https////ggttbb...m////acchhhhht hwtpTTakeechh-Lata SI- COT-Group/stt sarffor___ff___//////e///eaee/e_ fwwefffeecy-domaatee change_ point_ ceatt_deeectiin. Finally, Section 6 concludes this paper.
Note that while the "amplitudo" spectra presented in the figures of this paper visualize the temporal variations of spectral sequences and the means for each sogment, the "complex" spectra are utilized in the actual CP detection and the hypothesis testing (sce Section 2 and the subsequent scctions for details).
In Figure 1, the values labeled as naive p-vuuues are inappropriate because they fail account for dopple-dipping, leading to excessively small values and an inflated rate of false positive findings.
2

---

Published in Transactions on Machine Learning Research (06/2025)
de
r r
d
d2
dP
Un Mww.mmmmmm
wm
Value
do
do
20
60
80
100
20
40
60
80
100
Time
Time Point
Time Point
(a) Time series signal.
(b) CP detection result for each frequency component.
(c) CP detection result after in- tegrating CPs across multiple froquencics.
(d) Statistical Lest for the CPs detected in (c).
Time Point
Detection
True
True
Wrong
prvalls
0.008
0.006
0.752
Naive p-value
0.000
0.000
0.000
Figure 1: Demonstration of the proposed method. Panel (a) shows the original time series signal. Panel (b)
illustrates the result of CP detection for each time variation of the five frequency components, and panel (c)
represents the CPs merged across multiple frequencies. In panel (c), the CPs for frequency dj at time point 28
and d2 at 30, and those for d2 at 68 and d3 at 74 merge into one CP at 28 (A) and 70 (B), respectively.
Actually, the CPs A and B are truely detected, while a CP C for do at 97 wrong detection. Table (d)
shows the results of statistical test for the CPs detected in (c). The proposed p-value is enough large for falsely detected CP C, while the naive metod causes false positive because the p-value is too small. Furthermore, the proposed method provides sufficiently small p-values for truely detected CPs A and B, which are statistically significant.
Our contribution. The main contribution of this study is the development of a statistical testing method, based on the framework of SI, to quantify the statistical significance of detected CPs in the frequency domain. While CP detection in the frequency domain is important problem in engineering and several algorithms have been proposed, to the best of our knowledge, no existing method properly quantifies the statistical significance of CPs in the form of p-values. By leveraging the recently advancing SI framework, we provide a solution to this unresolved issue. There are various CP detection algorithm in the frequency domain. In this study, as an example, we consider a CP detection algorithm that identifies changes common across multiple frequencics using dynamic programming (DP) and simulated anncaling. should be noted that our contribution does not lie in this CP detection algorithm itself, but in proposing a statistical test to quantify the significance of the detected CPs.
3

---

Published in Transactions on Machine Learning Research (06/2025)
Related work. The CP detection problem has long been studied with various applications in variety of fields, such as finance (Fryzlewiez & Subba. Rao, 2014; Pepelyshev & Polunchenko, 2017), bioinformat- ies (Chen & Wang, 2008; Muggeo & Adelfio, 2011; Pierre-Jean et al., 2015), climatology (Reeves et al., 2007; Beaulieu et al., 2012), and machine monitoring (Lu et al., 2017; 2018). In the statistics and machine learning communities, various methods have been proposed for identifying multiple CPs from univariate sequences. The most straightforward approach is repeatedly applying single CP detection algorithms, such as binary segmentation (Scott & Knott, 1974). Examples of such approaches include circular binary segmentation (Ol- shen et al., 2004) and wild binary segmentation (Fryzlewicz, 2014). Another line of research has proposed numerous approaches based on penalized likelihood, such as Segment. Neighbourhood (Auger & Lawrence, 1989), Optimal Partitioning (Jackson et al., 2005), PELT (Killick et al., 2012), and FPOP (Maidstone et al., 2017). In these studies, dynamic programming is employed to solve the penalized likelihood minimization problem. While most CP detection studies focus on the time domain, a few studies have targeted the froquency domain, such as Adak (1998), Last & Shumway (2008), and Preuss et al. (2015).
Most existing studies for statistical inference on detected CPs have relied on asymplotic theory under re- strictive assumptions, such weak dependency. For single CP detection problems, approaches such the CUSUM score (Page, 1954), Fisher discriminant score (Mika et al., 1999; Harchaoui et al., 2009), and MMD (Li et al., 2015) conduct statistical inference on the detected CPs based on asymptotic distribution of some discrepancy measures. For multiple CP detection problems, methods such as SMUCE (Frick et al., 2014) and the MOSUM procedure (Eichinger & Kirch, 2018) also employ asymptotic inference. However, these methods primarily focus on testing whether a CP exists within a certain range rather than directly quantifying the statistical significance of the location of the CP itself. Furthermore, asymptotic approaches often fail to control the false positive (type I crror) rate effectively, resulting in conservative testing with low statistical power (Hyun et al., 2018).
SI was initially introduced as a method of statistical inference for feature selection in linear models (Taylor & Tibshirani, 2015; Fithian et al., 2015) and later extended to various fcature selection algorithms, inchuding marginal screening (Lee Taylor, 2014), stepwise feature selection (Tibshirani el al., 2016), and Lasso (Lee et al., 2016). The core concept of SI is to derive the exact null distribution of the test statistic conditional the hypothosis selection event, thereby enabling valid statistical inference with controlled type I error rate. Recently, significant attention has been given to applying SI to more complex supervised learning algorithms, such as kernel models (Yamada et al., 2018), boosting (Rügamer & Greven, 2020), tree-structured models (Neufeld et al., 2022), and neural networks (Duy et al., 2022; Miwa et al., 2023; Shiraishi et al., 2024b). Furthermore, SI has proven valuable for unsupervised lcarning tasks, including clustering (Lee et al., 2015; Chen & Witten, 2023; Gao el al., 2024), outlier detection (Chen & Bien, 2020; Tsukurimichi et al., 2022), domain adaptation (Duy et al., 2024), and segmentation (Tanizaki et al., 2020; Duy et al., 2022). SI was first used for CP detection problem in Hyun et al. (2018), which focused on Fused Lasso algorithm. Since then, the framework has been explored in various CP detection algorithms, including the CUSUM- method (Umezu & Takeuchi, 2017), binary segmentation and its variants (Hyun et al., 2021), dynamic programming (Duy et al., 2020), and other related problems (Sugiyama et al., 2021; Jewell et al., 2022; Carrington & Fearnhead, 2024; Shiraishi el al., 2024a). Existing studies on CP detection have focused on the time domain, making this the first to offer valid statistical inferences for frequency- CPs using the SI framework.
2 Problem Setup
In this section, we first desrcribe the probabilistic model of time series on which the statistical inference is based, and then formulate the problem of CP detection in the frequency domain.
2.1
Probabilistic Model for Time Series Data
For statistical inference, we interpret that the observed time series data is realization of a random sequence following a certain probabilistic model. Let denote the univariate random sequence with length N by
X= = (X1, .... XN)T = =ste, o0 TNT),
(1)
4

---

Published in Transactions on Machine Learning Research (06/2025)
where E RRN is the unknown true signal vector, and E RN is the normally distributed noise vector with the covariance matrix INIx= 3. In other words, we assume that the observed time series data is sampled from the probabilistic model X ~ (s,,2xxx
Then, consider applying discrete short-time Fourier transform (STFT) to the random sequence X using rectangular window of width M without overlapping. In this case, the computations of DFT are required M i.e., N is a multiple of M for simplicity. We denote M-point
T N times, where we assume T
DFT matrix
WM W (0)
M ..... WM ( MMM
where why = (u ) and war e - for frequency d E {0,...,
(M-1)d M
............ M - 13 in which j is
M
the imaginary unit. Using whr we consider multiple spectral sequences across D due to the symmetry property of the spectrum. For frequency d, the sequence of spectra is written as
M + frequencies 2
F(4) = (F( F(1)) E CT, de (O...,,D- D-13.
where F((0) = (1ms @ (4)
T
X, and ss:e E IR is a vector whose elements from position s to e are set to 1,
and 0 otherwise for 1 S T.
2.2 Statistically Significant CP Detection in the Frequency Domain
As mentioned above, for statistical inference, we interpret that the observed time series is randomly sampled from the probabilistic model in (I), which is denoted by
x = (xI,...,,NN) ERR.
(2)
Similarly, the sequence of spectra for frequency d, obtained by applying the aforementioned STFT to the observed time series æ, is written as
y(0)= = [11").....f
E ECT, de DO.....D-11)
The goal of this study is to detect changes in the true signals of frequency spectral sequences {F(d) (ee based on the observed sequences ) (o.......
Among variety of changes, we focus in this paper on mean-shift of frequency spectral sequences. Let us denote the mean spectrum of frequency at time point by
m = E[ (d. t) E tO....,. D - 1} x [T].
where the expectation operator E[] is taken with respect to the probabilistic model in (1)4, and [T] [!, indicates the set of natural numbers up to T.
Considering segment from time points to e with 1 i e < T, we say that there is a mean-shift change in the frequency d at time point t E fs....,e - 1} if and only if
Σ V=8
( Mr c
Σ
(4) Pr
8
v=++1
As discussed in Section have prior knowledge that simultaneous changes occur across multiple distinct frequencies. In Section 3, we introduce a heuristic algorithm that generates CP candidates by incorporating this prior knowledge. Subsequently, in Section 4, we present an SI framework to quantify the statistical significance of each CP candidate in the form of p-values. Finally, we select the candidates with p-values smaller than the usor-specified significance level (c.g., 0.05 or 0.01) as our final CPs.
IIn the following discussion, the noise must be independent and identically distributed Gaussian vector with prede- termined covariance matrix. The robustness of our proposed method for unknown noise variance, non-Caussian noise, and correlated noise discussed in Appendix D.3.
"Since the frequency spectral sequences obtained through linear transformations of Ganssian random variables, the existence of their expertations is guaranteed.
5

---

Published in Transactions on Machine Learning Research (06/2025)
3 Heuristic Algorithm for CP Candidate Selection in Multiple Frequencies
As discussed in Section 1, the types of problems occurring in time series signals can be systematically determined by simultancously detecting CPs across multiple frequencies, such as the harmonics and sidebands of the characteristic frequencies. Thus, we formulate CP candidate selection as an optimization problem that aims to not only estimate the optimal number and location of CPs for each frequency, but also reduce the total number of change locations across all frequencies by aligning the positions.
3.1 Objective Function for CP Candidate Selection
Let ((4) be the number of selected CP candidates and 7 = (T(4) (d)
(d) } S [T-1] be the ordered set of R())
CP candidates (7(d) <... < 714) K(o and TO 0, K( = T) for frequency d E .O....,D D D--1 1}. Furthermore, we define the set of CP candidate locations as
D-1
U d=0
7 (4)
1 [T, (T.....,,xx Tx] S [T- -
(4)
where we denote the total number of CP candidate locations as K == |T|-
Let T = (7 ,..,, ((-11)) be the collection of all CP candidates across all the D frequencies. Then, given an observed time-series data æ in (2), the objective function of our CP candidates T written as
D-1 D-1K((41
E(T m) x)
=
Σ
1 (1 ) (
(
c
3(() + 7K.
(5)
d=0
k=1
The first term of the obejective function in (5) indicates the cost for quantifying the variability of segments between two adjacent CP candidates, where C (faa) indicates the cost function for a segment f!d (J!4) ((c) f. .! specifically defined as
c(110) =
(d) sym
(d)
f
8
t=
1'=8
where
1
7
(6)
sym
2 if if 0, 2
M
because the frequency spectra have complex conjugate symmetry. The second term indicates the penalty term for the number of CP candidates in each frequency, while the third term indicate the penalty for the number of total CP candidate locations, where (3(0) (..,,,,(DDD1)) E RD and y E R are hyper-parameters for controlling the balance between the three terms (the details of how to determine these hyper-parameten are presented in Appendix A). The third penalty term indicates that we take into account the trade-off between cost and penalty term for not only ((d) but also K in (5), hence CP candidates of multiple frequencies tend be detected the same time point.
3.2 Approximately Optimizing the Objective Function by Simulated Annealing
Unfortunately, since minimizing the objective function in (5) is challenging combinatorial optimization problem, we must rely on approximate solutions derived from heuristic algorithms. Following the approach in Lavielle (1998), employ simulated annealing (Kirkpatrick et al., 1983; Černý, 1985) to approximately solve the optimization problem. Simulated annealing is meta-heuristic algorithm widely applied to various practical problems, as it converges asymptotically to global solution with high probability under specific conditions. To approximately optimize the objective function in (5), we perform the following two steps:
Step Individually estimate the CP candidates for each frequency. Specifically, we minimize the objective function (5) with ? = 0. This can be optimally achieved by applying dynamic programming each sequence.

---

Published in Transactions on Machine Learning Research (06/2025)
Step 2: Refine the CP candidates for each frequency estimated in step using simulated annealing to estimate the changes shared across multiple frequencies.
The reason for obtaining approximate solutions using simulated annealing is that the global minimization of the objective function in (5) across multiple frequencies cannot be achieved by directly extending the dynamic programming approach used in step 1. This inability to obtain the global optimum stems from the fact that it is impossible to construct a Bellman equation for this optimization problem; specifically, additional information which is not preserved in the solutions of subproblems may be needed at later-stage decisions.
We note that approximately solving the combinatorial optimization problem in (5) is NOT our novel con- tribution (our key contribution, detailed in Scction 4, lics in providing theoretical guarantee for the false positive detection probability of the obtained approximate solution). For example, similar approach has been used for analogous problems in studies such as Lavielle (1998). Moreover, we do NOT claim that simulated annealing is the optimal approach for this problem; it is simply one of the reasonable choices, and other meta-heuristics could also be used as alternatives.
3.3 Step 1: Generating an Initial Solution by Dynamic Programming
We need to generate initial solution T (r init To obtain the initial solution, we first set 0 in (5) to detect CP candidates for each frequency and formulate an optimization problem as
...., 7 (iit(D-))) before applying simulated annealing.
KG0 +1
(iii((() = arg min (d)
Σ k=
C
(1))
+ g(0) +,000 k(d)
(7)
This optimization problem can be solved efficiently using dynamic programming algorithm which is called Optimal Partitioning (Jackson et al., 2005). This method employs the Bellman equation that recursively determines optimal solutions for simpler subproblems. Given Dinit the set of frequencies for which least one CP candidate is detected, it is sufficient to apply simulated annealing only to d E Dinit because introducing an additional CP to any frequency d requires a minimum penalty of 3(d) as shown in (5). Thus, we reduce the number of optimal solution candidates, and the computational cost. of simulated annealing can decreased without degrading the solution quality.
3.4 Step 2: Refining the Solution by Simulated Annealing
We perform multivariate CP candidate selection in the frequency domain using simulated anncaling that is applied for solving large combinatorial optimization problems for which finding global optima is difficult. In each step of simulated annealing, the Metropolis algorithm is used to accept transitions not only to improving solutions that decrease the objective function, i.e., AE(T', T. = E(T x) E(T, x) 0, where T and T' are current and new solutions, respectively, but also to deteriorating solutions that increase the objective function, i.e., (T', T, x) > 0, with the probability controlled by a temperature parameter e, which allows escape from the local solution. The pscudo code of the Metropolis algorithm is shown in Algorithm 1.
Algorithm metropolis_algeree
Input: Difference of objective function values AE, temperature 1: Uniformly sample & from 1),
2: if AE + clng < 0 then
3: status + Acceptance
4: else
status + Rejection
end if
Output: status

---

Published in Transactions on Machine Learning Research (06/2025)
Local search. In this paper, we consider four types of neighborhood operations applicd to the current solution, T.e., adding, removing, and moving a CP (Tavielle, 1998) for a uniformly selected frequency d, and merging two adjacent CP locations uniformly sampled from at a random position between them. Schematic illustrations of these four operations are provided in Figures 2 and 3. We set the number of searching iterations al a specific lemperature e equal to the size of neighborhoods, i.e., |Dinit| T (Aarts & Korst, 1989). Each operation is uniformly selected from adding, removing, and moving one CP. Subsequently, two adjacent CP locations are merged only once because this operation significantly fluctuates the objective function value. When the operation is not possible (e.g., removing operation for frequency with no CP), it is skipped. If no transition to neighborhoods occurs when these operations are repeated sufficiently at a certain temperature e, the search is terminated.
Hall y mnw
mummyn Mhn you herh
v M
30
68
30
68
30
68 Time Point
Time Point
Time Point
the you 1/
de
v. Mmm
5768 Time Point
68
30
73
Time Point
Time Point
(a) Adding at 57.
(b) Removing at 30.
(c) Moving from 68 to 73.
Figure 2: Illustrations of the three local search operations for frequency d. (a) Adding is to insert CP at randomly selected time point with no CP. (b) Removing means to delete CP at randomly selected time point with a CP. (c) Moving is to shift a CP that is randomly selected from ((d) to a random position between its adjacent CPs.
r N A. Myyn
d2
wmw M Mww.M
www. M Www.mmooao..
do
An M 97
28
68 74
97
28 Time Point
Time Point
Figure 3: Mustraaioo of the local search operation which merges two adjacent CP locations. In this figure, CP for frequency de at time point 68 and CP for d3 at 74 are merged at 70.
8

---

end if Published in Transactions on Machine Learning Research (06/2025)
Setting of initial temperature. We determine an initial temperature co by setting acceptance ratio
v(co) # accepted transitions
# proposed transitions
to desired value in preliminary experiment of simulated annealing. In practice, start by setting the temperature to a sufficiently small positive value, then multiply it with a constant factor X larger than 1, as follows:
-
where of represents the i-th temperature in the initial value setting, until the acceptance ratio exceeds the predefined criterion (Aarts & Korst, 1989). Therefore, the imitial temperature co is specified as the final value ct.
Decrement of temperature. We employ geometric cooling schedule, which is used as practical method of temperature control, altough it does not guarantee the convergence to global optimum. The decrement function of the i-th temperature & is given as
C++1
Aci,
where À is a constant factor smaller than but close to 1. The value typically lies between 0.8 and 0.99 (Aarts & Korst, 1989). That is because the decreasing rate must sufficiently slow to obtain a better solution for the optimization problem.
The overall procedure of the CP candidate selection using simulated annealing is shown in Algorithm 2.
Algorithm CP Candidate Selection by Simulated Annealing
Input:
Time sequence æ Obtain spectral sequences (()) for d E (O...., 1} in (3) by applying STFT to
Initialize the CP candidates T as T in (7) using dynamic programming
Determine the initial temperature in preliminary experiment
4: while true do
5: while The number of local search less than [Dinit T do
6:
Frequency d is uniformly sampled from Dinit
7: The operation is uniformly sampled from adding, removing, and moving a CP for Obtain T' by applying the operation status metropolis, iattroplls_llloriitttttoooeoeeei thmm (AA(AAFTT T)) 10: if status is Acceptance then
11:
T + T end if
12:
13:
end while
14:
Obtain T' by merging uniformly sampled two adjacent CP locations T
15:
status metropolis metropolis__lgorrttua AagorithmAAE(( e))
16: if status Acceptance then
17:
T+
18:
19: if no transition to neighborhoods occured at c then
20: break
21: end
22: Update the temperature using a constant factor À c x- Ac
23:
end while
Output: CP candidates T

---

Published in Transactions on Machine Learning Research (06/2025)
4 Selective Inference on CP Candidate Location
In this section, using the SI framework, we quantify the statistical significance of all locations selected as CP candidates by the algorithm in Section 3 in the form of p-values. By setting the significance level a (e.g. 0.05 or 0.01) and considering CP candidate locations with the p-values below a as the final CPs, it is theoretically guaranteed that the false positive detection probabilities (type error rates) of the final CPs are controlled below the specified significance level
To formalize the SI framework, let write the CP candidate selection algorithm in Section 3 as A: X H T 5. As detailed in Section 3, the output of this algorithm is collection of CP candidates represented T --((), ..... . (r())...,=((-1)) where each x(d) (T1 -...) for frequency d E {O..... - 1}. The final set of selected CP candidate locations is obtained by taking the union of these sets of CP candidates, i.c., T == U2_ 7(4) [TI...., TK} S [T 1]. Subsequently, we quantify the statistical significance of each CP candidate location ты E for k E [K], in the form of p-values within the SI framework. Thercfore, without loss of generality, formulate the hypothesis testing problem for the k-th CP candidate location 7he for simplicity.
(d) } S [ 1] is the ordered set of CP candidates
4.1 Statistical Test on CP Candidate Location
Hypotheses. To test the statistical significance of CP candidate location Tk, consider whether mean- shift occurs before and after this location тк at any frequency for which this location was selected as CP candidate. To formalize this, let Dx denote the set of frequencies where тı is selected as CP candidate (d) location, ie.. D2 {d {O..... D 1} T (()). For each d E Dk, let Tpre and (d) represent the CP candidate locations immediately before and after ты at frequency d, respectively. Note that although Tpre
and (d) depend The this dependence is omitted from the notation for simplicity, Given these notations, we formulate the null hypothesis Ho. and alternative hypothesis H1,k as follows:
Ho
(d)
H4
(d)
(d) Tsue
(d) vd DE.
(8)
TK
Tpre
(=mpre
+1
Tk
(=
v.8.
(d)
HI
(d) 14
(d) 3d E Dk-
(9)
(d)
(d)
Tk
Tpe
(=
+1
Tsuc
Tk
t=
Test statistic. We define the test statistic to test with the null hypothesis (8) and alternative hypothesis (9) as follows. First, for each frequency complex spectra in the two segments before and after the location ткe Then, we aggregate these differences across all frequencies d E D3 to define the test statistic. Therefore, the test statistic Ta (X) is formulated
Dx, we consider the difference between the average values of the
(d) (d) a len Gayan M
F(d) 2
Tx(X)= 0 -1
Σ D
F ((d)
(10)
Tpre +pr
+1:7k
Te +c+1:
where den = (T - TT) ((d)) (TTuc
Tpre) E is a scaling factor for correcting the segment lengths (d) before and after the location тı for each frequency DE, caym E R defined in (6) is required for considering
the complex conjugate symmetry of the spectrum, and ((d) E C is the average value of the complex spectra in the segment from time points s to e for frequency E Dr. ie.,
((d)
1 - 5 s
(Iere & w (1)
X.
(11)
c-5
Note that, to deterministic hehavior the algorithm A, the random seed is fixed to constant value at the heginning of the procedure.
10

---

Published in Transactions on Machine Learning Research (06/2025)
An important point when conducting the statistical test within the SI framework is that the test statistic is represented as the norm of projection of the data, i.e. (x P where Ps E RNNN is a projection matrix that depends on the observed time series x only through the selected CP candidates A(x) (see the next subsection for further details). We now define the orthogonal projection matrix Px as
(d) (d)
Ps
=
a Jen Caym
E NRN
(12)
M
dEDx
where
o (d)
(d) Tpre
ppre
+1:r
(c)
(d) Tsuc
1
11+1:
Tx+
(a) 2 w
E EC
(13)
Th
7
is vector that depends the selected CP candidates A(æ), and v(d)) is the complex conjugate of v(d)
Using (10) (13), the test statistic Ta (X) can be rewritten as
(d) (d) len Caym
(x) --
(ol) X
dEDx
(d) (d)
Hen Caym M
=0-
X
X
-
Ps (
(14)
4.2 Computing Selective p-value
Selective Inference (SI). To compute the p-value, need to identify the sampling distribution of the test statistice (x..). However, the projection matrix Pa within the test statistic ii Tx (X) depends the CP candidates A(X) (which, in turn, depends on the sequence X through the CP candidate selection algorithm A), the sampling distribution of the test statistic T₃ (X) is too complicated to characterize. Within the SI framework, to address this challenge, we consider the sampling distribution of the test statistic conditional on the event that the selected CP candidates A(X) for random sequence X the same as A(x) for the observed sequence æ, that is,
Ta (X) (A(X) A(z)).
(15)
To compute a selective p-value based on the conditional sampling distribution in (15), we introduce an additional condition on the sufficient statistic of the muisance parameter Q(X), which is defined as
Q(x) = (V(X), u(X)) E R. x RRN
(16)
with
OP,
V(X)
E ENN u(x) (IN- ER..
| X|
This additional conditioning on Q(X) is a standard approach for computational tractability in the SI liter- ature (Loftus & Taylor, 2015). Based on the additional conditioning on Q(X), the following theorem tells that the conditional distribution of the test statistic can be represented as a truncated x-distribution.
Theorem 1. Consider a random sequence X ~ o(x, and an observed sequence æ. Let A(X) and A(x) the detected CP candidates, by applying a CP candidate selection algorithm to X and x, respectively. Let Ps E R NxN be projection matrix depending on A(x), and consider a test statistic in the form of Tx(X) o (P,x Furthermore, define the nuisance parameler ((X) (V(X), U(X)) in (16).
Then, under the null hypothesis that the norm of the true signal vector projected by Ph is zero (i.e., ||Pas|] 0), the conditional distribution of the lest slatistic
Tx(X) (A(X) A(ar), 2(x) 2(x))
11

---

Published in Transactions on Machine Learning Research (06/2025)
is truncated X-distribution (C()) 2) with the degrees of freedom tr(Pk) and the truncation region 2. The truncation region 2 is defined as
2 {2 [ A + bæ) A(z)} U(x), V(x).
The proof of Theorem is deferred to Appendix B.1. By using the sampling distribution of the test statistic (X) conditional on A(X) A(x) and ( (2x) in Theorem 1, we can define the selective p-value as
Pk
seloctive = PHo. (T&(X) 2 (x) A(X) A(x), 0(x) Q(x)).
(17)
Theorem 2. The selective p-value defined (17) satisfies the property:
PHo. (eelectivee SalA(X) -(x))==())))))))))))) a, Va E (0, 1).
Then, the selective p- pavul also satisfies the following property of valid p-value
HH (prelective
<a) = =a Va E (0.
The proof of Theorem 2 is deferred to Appendix B.2. This theorem guarantees that the selective p-value is uniformly distributed under the null hypothesis Ha,k, and thus be used to conduct the valid statistical inference in (8) and (9). Figure 4 schematically illustrates the SI framework presented in this section.
Once the truncation region 2 identified, the selective p-value in (17) can be easily computed by Theorem 1. Thus, the remaining task is roduced to identifying the truncation region 2. this study, to identify the truncation region 2, we adopt the method based on parametric-proorammiii technique proposed by Duy & Takeuchi (2022). For further details, refer to Appendix C.
Test statistic (x))
(ax)
Parametrized line X by
Truncation region Z
C
A(x)
A(x²)
A(xx)
A(w)
A(x:)
Data space R
Figure 4: Schematic ralustration of the SI framework. point in the data space R corresponds to a sequence with length N. The darkly shaded regions the data space indicate that, if we input a point in these regions into the algorithm A, the CP candidates are the same as A(x) obtained from the observed sequence æ. By conditioning on these regions and Q(X)= = Q(z), the conditional sampling distribution of the test statistic (x( represented truncated x-distribution. Thus, selective p-values defined based the tail probability of such truncated x-distribution.
12

---

Published in Transactions on Machine Learning Research (06/2025)
5 Numerical Experiments
5.1 Methods for Comparison
prelective in (17) with the following In our experiments, we compared the proposed method (Proposed) using pk methods in terms of lype I error rale control and power.
OC: In this method, that is, a simple extension of ST literature to our setting, we consider p-values with additional conditioning (over-conditioning. described in Appendix C.1.
OptSeg-SI-oc, OptSeg-SI (Duy et al., 2020): These methods use p-values conditioned only on the dynamic programming algorithm, disregarding the conditioning on simulated annealing.
Naive: This method is a conventional statistical inference.
Bonferroni: This method applies Bonferroni correction for multiple testing correction.
The details of these comparison methods are provided in Appendix D.1.
5.2 Synthetic Data Experiments
Experimental setup. In all synthetic experiments, we set window size M E {512, 1024), the number of froquencies D + the length of sequence N = M T, where T specified for each experiment, the sampling rate Sa = 20480, and each element of mean vector as
Sn Σ
A sin (w) (I
dE(di ds)
1} were randomly selected without replacement for each simulation, A/4)) was defined for each experiment, and (()) {22(( L )d d= o...., [ We used BIC for the choice of penalty parameters 3 and y indicated in Appendix A, and set the parameters of simulated annealing ot 1000, X+ = 1.5, n)(o))= 0.5 and X = 0.8 in Section 3.4. After detecting CP candidates, a CP candidate location relt randomly seleeted from det was tested at the significance level a = 0.05.
where frequencies d1, d2, d3 E {O...., D D-1) 1}
In the experiments conducted to evaluate the type I error rate control, we generated 1000 null sequences, which did not contain true CPs in the frequency domain, = (z.....,..) ~ N (s. Is.o"IN), where A A (d) was randomly sampled from [0, 1) for d in each simulation, and for each T' (40, 60, 80, 100). Regarding the experiments to compare the power, generated sequences æ (IN) ~ N (8.) where
A(d)
f .
A/.0
A(d)
+
A
if M ( 1 M M. 1 if 1/ + +1<tET
AA4)
2A
with ((ds), E [0, 1) which were randomly sampled in each simulation, (A ((a)) .t ((db))
(18, 20, 22), (d))) (d2) (ds) (38, 40, 42), an intensity of the change A E (0.04, 0.08, 0.12, 0.16) and 1, for T' = 60. In each case, ran 1000 trials. Since tested only when CP candidate location correctly detected, the power was defined as follows:
Power (or Conditional Power) # correctly detected & rejected
# correctly detected
We considered the CP candidate location Tkk to be correctly detected it satisfied the following two conditions:
The set Dj of frequencies containing at least one CP candidate was subset of (di, da).
For DA satisfying the above condition, either min (d) de 1 S S maxt or min (d)
deD, deD, t2 S S maxt (d) deD,
(d)
held, that is, Tk was detected within the true CP locations for the frequencies in Dj.
13

---

Published in Transactions on Machine Learning Research (06/2025)
Experimental results. The results of experiments regarding the control of the type I error rate are shown in Figure 5. The Proposed, OC, and Bonferroni successfully controlled the type error rate below the significance level, whereas the OptSeg-SI, OptSeg-SI-oc, and Naive could not. That was because the OptSeg-SI and OptSeg-SI-on SI-oc used p-values conditioned only on the dynamic programming algorithm, excluding the conditioning on simulated annealing, and the Naive employed conventional p-values without conditioning. Since the OptSeg-SI, OptSeg-SI-oc, and Naive failed control the type error rate, we omitted the analysis of their power. The results of power experiments are shown in Figure Based on these results, the Proposed was the most powerful of all methods that controlled the type I error rate. The power of the OC was lower than that of the Proposed due to redundant conditions (see Appendix C for details), Furthermore, the Bonferroni method had the lowest power because it was a highly conservative approach that accounted for the huge number of all possible hypotheses. Additionally, we provide the compulational time the Proposed in both experiments and the information on the computer resources in Appendix D.2.
1.0-
1.0
Proposed
Propesed
OC
OC
0.9
Bonterroni
Bontemonni
Type Error Rate
OptSeg-SI
Type Error Rate
OptSeg-Si
0pt5eq-Si-oc
OptSeg-SS-oc
0.3
Nalve
0.3
Nalve
0.2
0.2
0.1
0.1
0.0
0.0
40
60
80
100
40
60
80
100
(a) M 512.
(b) M 1024.
Figure 5: Type I error rate.
1.0
Proposed
1.0
Proposed
OC
Bonferroni
Bonferreni
0.8
0.8
0.6
0.6
Power
Power
0.4
0.4
0.2
0.2
0.0
0.0
0.04
0.08
0.12
0.16
0.04
0.08
0.12
0.16
A
(a) M 512.
(b) M 1024.
Figure 6: Power.
14

---

Published in Transactions on Machine Learning Research (06/2025)
Additional experiments. We also conducted the following experiments to investigate the robustness of the Proposed in terms of type I error rate control.
Unknown noise variance: We considered the case where the variance x² was estimated from the same data.
Con-Gaassian noise: We also considered the case where the noise followed the five types of stan- dardized non-Gaussian distributions.
Correlated noise: Furthermore, we considered the sequence whose noise was correlated, i.e., the covariance matrix Σ * o2IN. In this case, although the test statistic did not theoretically follow x-distribution, we conducted the hypothesis testing using our proposed framework.
In addition, a sensitivity study for the penalty parameter y within the objective function in (5) was also performed. These details and results are shown in Appendix D.3 and D.4, respectively.
5.3 Real Data Experiments
To demonstrate the practical applicability of the Proposed, we applied the Proposed, DC, and Naive to a real-world dataset. We used the set No.2 of the IMS bearing datasel, which is provided by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati (Qiu et al., 2006) and is available from the Prognostic data repository of NASA (Lee et al., 2007). The experimental apparatus consisted of four identical bearings (bearings 1, 2, 3, and 4) installed on a common shaft, driven at a constant rotation speed by an AC motor under applied radial loading. In this dataset, the vibration signals were measured using accelerometers until the outer race of bearing failed at the end of the experiment, as shown in Figure 7, The signal with the sampling rate f. = 20480 was recorded for one second at intervals of 10 minutes over a period of about 7 days for each bearing. The analysis for the signal of bearing in the frequency domain, as conducted by Gousseau et al. (2016), had revealed that enhancements of the spectral intensity were detected in harmonics of the characteristic frequency (236 Hz) associated with the outer race failure (Ball Pass Frequency Outer race, BPFO) on 3-4 days. Based on this previous study, we conducted CP candidate selection for the sensor data, of bearing in the frequency domain for two periods: 0.25-2.25 days when significant changes in the spectra existed and 2.25 4.25 days when the BPFO harmonics exhibited the spectral amplification. For each period, we computed the DFT of M = 1024 consecutive points in the 20480 samples and repeated the procedure T = 2 days/10 min = 288 times. As a result of the DFT, we obtained components for 513 frequencies (0, 20, .... 10240 Hz) and conducted CP candidate selection using 131 frequencies within the range of 1400 4000 Hz, which included the 6th to 16th harmonics of the characteristic frequency. For other parameters related to the CP candidate selection algorithm, please refer to the experimental setup in Section 5.2. Subsequently, evaluate whether each of them was genuine CP location. The variance o² for testing was estimated from the data on 0-0.25 days that did not contain significant spectral changes and was not used in any of the experiments. The results for the signal of bearing 1 on 0.25 2.25 days and 2.25 4.25 days are shown in Figure 8. In panel (a), the time variation of a frequency spectrum (1920 Hz) where a CP candidate location was falsely detected is shown for the period of 0.25 2.25 days. It shows that p-values of the Proposed and DC are above the significance level 0.05, and therefore the result provides the validity of the inference, while p-value of the Naive too small. In panel (b), the time variations of the Sth and 15th harmonies of the BPFO where CP candidate locations were correctly detected are presented for the period of 2.25-4.25 days. In this case, p-values of the Proposed are below the significance level 0.05 0.025 decided by Bonferroni correction, thus it indicates that the inference is valid. In contrast, p-values of the OC are too large, due to the loss of power caused by the redundant conditions. In addition, since even the time sequences of the healthy bearings 2, 3, and 4 had been reported to indicate spectral amplification associated with the outer race fault in bearing 1 (Gousscau et al., 2016), we performed the same analysis for the three signals. The results are shown in Appendix D.5.
tested the detected CP candidate locations to
15

---

Published in Transactions on Machine Learning Research (06/2025)
Amelludde
Amplituds
Amgiludd
Amplitudo
Time (Day)
Time (Day)
Time (Day)
Time (Day)
(a) Bearing 1.
(b) Bearing 2.
(c) Bearing
(d) Bearing 4.
Figure 7: Vibration signals of the four bearings the No.2 of the IMS bearing dataset. Each signal was measured using an accelerometer for about 7 days until the outer race of bearing 1 failed. We used not only the signal of the damaged bearing 1 but also those of the other healthy bearings for the analysis because the spectral amplification associated with the failure of bearing had been reported to be observed in the signals of all four bearings (Gousscau et al., 2016).
14
Proposed: 0.000 OC: 0.548 Naive: 0.000
p-values
12-
Amplitude Spectrum
8-
14
Proposed:
14
0.005
OC:
OC: 0.202
12
Naive: 0.000
12
Naiiee: 0.000
Amplitude Spectrum
Amplitude Spectrum
10.
10
0.25 0.50 0.75
1.00 1.25
1.50
1.75
2.00
2.25
2.25 2.50 2.75 3.00 3.25 3,50 3.75 4.00 4.25 (Day)
Time (Day)
(a) Inference on a falsely detected CP candidate location for 1920 Hz (around the 8th harmonic) on 0.25 2.25 days.
(b) Inferences on truely detected CP candidato locations for 1880 Hz (the 8th harmonic, above) and 3540 Hz (the 15th harmonic, below) on 2.25- 4.25 days.
Figure 8: Results of the CP candidate selection for the signal of bearing in the frequency domain and the subsequent inference for the detected CP candidale locations. In panel (b), note that p-values for the first. CP candidate location were actually computed by considering not only CP candidate of the 15th harmonic but also a CP candidate of 1900 Hz (around the 8th harmonic).
16

---

Published in Transactions on Machine Learning Research (06/2025)
6 Conclusion, Limitation, and Future Work
In this paper, we developed statistical inference method to quantify the statistical significance of detected CP locations in the frequency domain. Our proposed framework contributes to various fields where time- frequency analysis is widely employed, such as condition monitoring of machine systems using vibrational, electrical, and acoustic signals, and medical diagnosis based on biosignals. We conducted comprchensive experiments on both synthetic and real-world datasets. The results theoretically confirmed that our method provided an unbiased evaluation based on SI framework and demonstrated its superior performance com- pared to existing methods. However, an important limitation is the assumption that the noise follows uncorrelated normal distribution, which is technically essential for deriving the conditional sampling dis- tribution of the test statistic as a truncated x-distribution. In future work, we will extend our met to the case of multi-dimensional sequences. For instance, analyzing different types of time series obtained from multiple sensors would reveal significant changes in the system unattainable through the univariate approach, and reliability guaranlee for the detections also provide valuable future contribution.
Acknowledgments
This work was partially supported by JST CREST (JPMJCR2ID3, JPMJCR22N2), JST Moonshot R&D (JPMJMS2033-06), and RIKEN Center for Advanced Intelligence Project.
References
Emile Aarts and Jan Korst. Simulltedd annealing and Boltzmann machines: a stochastic approach combi- natorial optimization and neural computing. John Wiley & Sons, Inc., 1989.
Sudeshna Adak. Time-dependent spectral analysis of nonstationary time series. Journal of the American Slatistical Associalion, 93(444):1448811 1401,1 1998.
Ivan E Auger and Charles Lawrence. Algorithms for the optimal identification of segment neighborhoods. Bulletin of mathematical biology, (1):):39 54, 1989.
Claudie Beaulieu, Jie Chen, and Jorge L Sarmiento. Change-poin point analysis as tool detect abrupt climate variations. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1962):1228 1228-11299 2012.
Rachel Carrington and Paul Fearnhead. Post-selection inference for quantifying uncertainty in changes in variance. arXiv preprint arxiv:2405.155 156700 2024.
Vladimír Corný. Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm. Journal of optimization theory and applications, 45:41 51, 1985.
Jie Chen and Yu-Ping Wang. A statistical change point model approach for the detection of dna copy number variations in array egh data. IEEE/ACM Transactions on Compulational Biology and Bioinformaties, 6 (4):529 541, 2008.
Shuxiao Chen and Jacob Bien. Valid inference corrected for outlier removal. Journal of Computational and Graphical Statistics, 29(2):323-333 2020.
Yiqun T Chen and Daniela M Witten. Selective inference for k-means clustering. Journal of Machine Learning Research, 24(152):1-41 2023.
Vo Nguyen Le Duy and Ichiro Takeuchi. More powerful conditional selective inference for generalized lasso by paramotric programming. Journal of Machine Learning Research, 23(300):1-37, 2022.
Vo Nguyen Le Duy, Hiroki Toda, Ryota Sugiyama, and Ichiro Takeuchi. Computing valid p-value for optimal changepoint by selective inference using dynamic programming. Advances in Neural Information Processing Systems, volume 33, pp. 11356 18866-13366 11367, 2020.
17

---

Published in Transactions on Machine Learning Research (06/2025)
Vo Nguyen Le Duy, Shogo Iwazaki, and Ichiro Takeuchi. Quantifying statistical significance of neural network- based image segmentation by selective inference. Advonces in Neural Informalion Processing Systems, 35: 31627 31639, 2022.
Vo Nguyen Le Duy, Hsuan- -Tee Lin, and Ichiro Takeuchi. Cad-da: Controllable anomaly detection after domain adaptation by statistical inference. In International Conference on Artificial Intelligence and Statistics, pp. 1828-1836. PMLR, 2024.
Birte Eichinger and Claudia Kirch. A MOSUM procedure for the cstimation of multiple random change points. Bernoulli, 24:526 564, 2018.
William Fithian, Jonathan Taylor, Robert Tibshirani, and Ryan Tibshirani. Selective sequential model selection. arXiv preprint orXic:1512.0256 2015.
Klaus Frick, Axel Munk, and Hannes Sieling. Multiscale change point inference. Journal of the Royal Statistical Society Series B: Statistical Methodology, (6(3):495-580) 2014.
Piotr Fryzlewicz. Wild binary segmentation for multiple change-point detection. The Annals of Statistics, 42(6):2243 2281, 2014.
Piotr Fryzlewicz and Suhasini Subba Rao. Multiple-haang--pinn detection for auto-regressive conditional heteroscedastic processes. Journal of the Royal Statistical Society Series B: Statistical Methodology, 76(5): 903 924, 2014.
Lucy L Gao, Jacob Bien, and Daniela Witten Selective inference for hierarchical clustering. Journal of the American Statistical Association, 119(545):332-3444 2024.
William Gousseau, Jérôme Antoni, François Girardin, and Julien Griffaton. Analysis of the rolling element bcaring data set of the center for intelligent maintenance systems of the university of cincinnati. In CM2016, Charenton, France, 2016.
Zaïd Harchaoui, Eric Moulines, and Francis R Bach. Kernel change-point analysis. In Advances in neural information processing systems, pp. 609-616, 616, 2009.
Sangwon Hyun, Max G'sell, and Ryan J Tibshirani. Exact post-selection inference for the generalized lasso path. Electronic Journal of Statistics, 12(1):1053 1009:1097999, 2018.
Sangwon Hyun, Kevin Lin, Max G'Sell, and Ryan J Tibshirani. Post-selection inference for changepoint detection algorithms with application to copy number variation data. Biometrics, 77(3):1037-1049, 2021.
Brad Jackson, Jefsrey D Scargle, David Barnes, Sundararajan Arabhi, Alina Alt, Peter Gioumousis, Elyus Gwin, Paungkaew Sangtrakuucharoen, Linda Tan, and Tun Tao Tsai, An algorithm for optimal partitioning of data on an interval. IEEE Signal Processing Letters, 12(2):105 108, 2005.
Sean Jewell, Paul Fearnhead, and Daniela Witten. Testing for a change in mean after changepoint detection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(4):1082-1104, 2022.
Rebecca Killick, Paul Fearnhead, and Idris Eckley. Optimal detection of changepoints with a lincar computational cost. Journal of the American Stalistical Association, 107(500):1590 1598, 2012.
Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi. Optimization by simulated anncaling. science, 220(4598):671 680, 1983.
Michael Last and Robert Shumway. Detecting abrupt changes in piecewise locally stationary time series.
Journal of multivariate analysis, 99(2)1191-214, 2008.
Mare Lavielle. Optimal segmentation of random processos. IEEE Transactions signal processing, 46(5): 1365 1373, 1998.
Jason D Lee and Jonathan Taylor. Exact post model selection inference for marginal screening. Advances in neural information processing systems, 27, 2014.
18

---

Published in Transactions on Machine Learning Research (06/2025)
Jason D Lee, Yuckai Sun, and Jonathan E Taylor. Evaluating the statistical significance of biclusters. Aduances in neural information processing systems, 28, 2015.
Jason D Lee, Dennis L Sun, Yuekai Sun, and Jonathan E Taylor Exact post-selection inference, with application to the lasso. The Annals of Statistics, 44(3):907 927, 2016.
Jay Lec, Hai Qiu, Gang Yu, Jing Lin, and Rexnord Technical Services. Bearing Data Set, NASA Ames Prog- nostics Data Repository (http://ti arc. nasa. patprjjccttppooooiiii-ddaaa--eee Ames Research Center, Moffett Field, CA, 2007.
NASA
Shuang Li, Yao Xie, Hanjun Dai, and Le Song. M-statistic for kernel change-point detection. In Advances in Neural Informalion Processing Systems, pp. 3366 3374, 2015.
Keli Liu, Jelena Markovic, and Robert Tibshirani. More powerful post-selection inference, with application to the lasso. arXiv preprint arXiv: 10X1:1100,,0077 2018.
Joshua R Loftus and Jonathan E. Taylor. Selective inference in regression models with groups of variables. arXiv preprint 1r11::1511.01440 2015.
Guoliang Lu, Yiqi Zhou, Changhou Lu, and Xueyong Li. A novel framework of change-point detection for machine monitoring. Mechanical Systems and Signal Processing, 83:533-544 2017.
Guoliang Lu, Jie Liu, and Peng Yan. Sraph-based structural change detection for rotating machinery monitoring. Mechanical Systems and Signal Processing, 99:73-82, 2018.
Robert Maidstone, Toby Hocking, Guillem Rigaill, and Paul Fearnhead. On optimal multiple changepoint algorithms for large data. Statistics and computing, 27:519 533, 2017.
Sebastian Mika, Gunnar Ratsch, Jason Weston, Bernhard Scholkopf, and Klaus-Robert Mullers. Fisher discrininant analysis with kernels. In Neural networks for signal processing LX: Proceedings of the 1999 IEEE signal processing sociely workshop (cal. no. 98th8468), pp. 41-48. Teee, 1999.
Daiki Miwa, Vo Nguyen Le Duy, and Ichiro Takeuchi. Valid p-value for deep learning driven salient region. Proceedings of the 11th International Conference on Learning Representation, 2023.
Vito M. R. Muggeo and Giada Adclfio. Efficient change point detection for genomic sequences of continuous measurements. Bioinformatics, 27(2):161 166, 2011.
Anna C Neufeld, Lucy L. Gao, and Daniela M Witten. Tree-values: selective inference for regression trees. Journal of Machine Learning Research, 23(305):1 43, 2022.
Adam B Olshen, E Seshan Venkatraman, Robert Lucito, and Michael Wigler. Circular binary segmentation for the analysis of array-based dna copy mumber data. Biostatistics, 5(4):557-572, 2004.
E. S. Page. Continuous inspection schemes. Biometrika, 41(1/2):100-113 1954.
Andrey Pepelyshev and Aleksey S Polunchenko. Real-time financial surveillance via quickest change- detection methods. Statistics and Its Interface, 10(1):93-106, 2017.
Morgane Picrro-Jcan, Guillem Rigaill, and Pierre Neuvial. Performance evaluation of dna copy number segmentation methods. Briefings in bioinformatics, 16(4):600 615, 2015.
Philip Preuss, Ruprecht Puchstein, and Holger Dette. Detection of multiple structural breaks in time series. Journal of the American Statistical Association, 110(510):654 668, 2015.
Hai Qiu, Jay Lee, Jing Lin, and Gang Yu. Wavelet filter-based weak signature detection method and its application on rolling element bearing prognostics. Journal of sound and vibration, 289(455):1066-100 2006.
19

---

Electronic Journal of Published in Transactions on Machine Learning Research (06/2025)
Jaxk Reeves, Jien Chen, Xiaolan L Wang. Robert Lund, and Qi Qi Lu. A review and comparison of changepoint detection techniques for climate data. Journal of opplied meteorology and climatology, 46(6): 900 005, 2007.
David Rügamer and Sonja Greven. Inference for 2- b-bossing. Statistics and computing, 30(2):279 289, 2020.
Andrew Jhon Scott and Martin Knott. A cluster analysis method for grouping means in the analysis of variance. Biometrics, pp. 507 512, 1974.
Tomohiro Shiraishi, Daiki Miwa, Vo Nguyen Le Duy, and Ichiro Takeuchi. Selective inference for change point detection by recurrent neural network. Newral Computation, pp. 1-33, 2024a.
Tomohiro Shiraishi, Daiki Miwa, Teruyuki Katsuoka, Vo Nguyen Le Duy, Kouichi Taji, and Ichiro Takeuchi. Statistical test for attention maps vision transformers. In Proceedings of the 41st International Confer- ence on Machine Learning, 2024b.
Ryota Sugiyama, Hiroki Toda, Vo Nguyen Le Duy, Yu Inatsu, and Ichiro Takeuchi. Valid and exact statistical inference for multi-dinensional multiple change-points by selective inference. orXiv preprint arXiv:2110.0800 2021.
Kosuke Tanizaki, Noriaki Hashimoto, Yu Inatsu, Hidekata Hontani, and lchiro Takeuchi. Computing valid p-values for image segmentation by selective inference. Proceedings of the IEEE/CVE conference on computer vision and pattern recognition, pp. 9553 9562, 2020.
Jonathan Taylor and Robert J Tibshirani. Statistical learning and selective inference. Proceedings of the National Academy of Sciences, 112(25):7629-7033 2015.
Ryan Tibshirani, Jonathan Taylor, Richard Lockhart, and Robort Tibshirani. Exact post-sclection inference for sequential regression procedures. Journal of the American Statistical Association, 111(514):600-B)) 6200
2016.
Toshiaki Tsukurimichi, Yu Inatsu, Vo Nguyen Le Duy, and Ichiro Takeuchi. Conditional selective inference for robust regression and outlier detection using piecewise-lineat homotopy continuation. Annals of the Institute of Statistical Mathematics, 74(6))1197--228,, 2022.
Yuta Umezu and Ichiro Takeuchi. Selective inference for change point detection in mulli-dimensional se- quences. arXiv preprint arXiv:1706.0051, 2017.
Chihiro Watanabe and Taiji Suzuki. Selective inference for latent block models. Statistics, 15(1)33137-3183, 2021.
Makoto Yamada, Yuta Umezu, Kenji Fukumizu, and Ichiro Takeuchi. Post selection inference with kernels. In International conference on artificial intelligence and statistics, pp. 152-160. PMLR, 2018.
A Determination of Penalty Parameters in the Optimization Problem
In this section, we derive the penalty term in (5) used for detecting CP candidates in the frequency domain. We begin by assuming that the true mean vector of f(d) for frequency d (O...., D-1) is piecewise constant as follows:
11 ~ CN ( (
(d) Tk-1= Tx .oof. te [ L,..,,,, ....,, Th (4)
(d)
ke [K(d] 1].
(18)
where Moeg [h-1 : T/ (4) represents the true mean of the k-th segment in frequency d, and 0f denotes the known variance such that o² Mo² by the property of DFT6 Then, we introduce BIC for deriving the
(d) (d)
6 Allhhough we impose the assumption of the piecewise constant mean structure for model selection based BIC, the theoretical validity of y-values obtained by our proposed methad guaranteed even when this assumption does not hold.
20

---

Published in Transactions on Machine Learning Research (06/2025)
values of the penalty parameters & in (5). Given the assumptions of (18), the unknown parameters 0 of the sequence f(() are expressed as
0
(T)
(d) Kia
(Tod) : ...... Msog ( (d)
.....
(T(()) [ x1))+++)
Therefore, the degrees of freedom can be computed as k(d) clyl (K( + 1). That is because the mean vector consists real numbers for frequency d 0, M and complex numbers for the other frequencies. Consoquently, the BIC of this model given by 2
BIC = -20gg + (K( +
((d) sym
(K + +1) log
where denotes the likelihood function for this model, and it can be rewritten by ignoring the terms that not contribute to its minimization as
K1d
BIC = C
Σ
((((d)
p(d)
11:(( ) (((5)) +1) Mo Mo
T) k(d).
Comparing the BIC with (7), 3(d) can be defined as
g = -clffn + +1) Mo logT.
(19)
While 3 can be determined based the BIC, to the best of our knowledge, there is theoretical method to define y. Therefore, referring to (19), we employ the value scaled by Mo log as 7 which is given by
? <Mo log
(20)
where k is a hyper-parameter determined by users based on problem settings. We heuristically set k = 0.5 in the synthetic data experiments, and и = 3 the real data experiment.
B Proofs
B.1 Proof of Theorem 1
Proof. According to the condition Q(X)= Q(x)), i.e., u(X) ((()) and V(X) V(xr), we have
u(X) = U(x)
* (IN P 1(xx)
X U(x)) V(X)=
X u(x) V(x)² (.-V(X) V(x))
# X a + bz,
where =((x), b V(x), and (xx oー P Then, we have
(X E N A(X) = A(x), 2(X) = ((x))
=(X E R A(X) = A(z), = a b2, R)
=(X + bz E R AA(a A(a+b=)) bz) A(x), R
=(X = a + E N | Z).
Therefore, by noting that | is zero, we obtain
T(X) (A(X) A(w), 2(x) 2(x)) TC((r(P)),
21

---

Published in Transactions on Machine Learning Research (06/2025)
B.2 Proof of Theorem 2
Proof. The sampling distribution of the test statistic conditional on A(X) A(æ) and @(x) Q(x) denoted by
Th (x) = A(x), 2(X) = 2(x))
is a truncated x-distribution with the degroes of freedom tr(Px) and the truncation region 2 defined in Theorem 1. Thus, by applying the probability integral transform, under the mull hypothesis,
prclective [A(X): A(x). Q(x) = Q(x) ~ Unif(0, 1),
which leads to
PHo.k (Pk (p S a A(X) = A(z), Q(x) = 2(x)) = a, Va E (0,
Next, for any
(0, 1), have
PH.k (DF (p S a A(X) A(x))
= Piook (Pk ( s a A(x), Q(X) = Q(x)) PHo,k (Q(X) Q(x) A(X) A(()))((x)
a / (2)x)(((X)) = Q(x) A(X) = A(z))dQ(x)
a.
Therefore, we obtain the result in Theorem as follows:
F k (p%
( S a)
Σ
PH (pk (peleective <a A(X) A(x)) (H.,)(A(X) = A(x))
A(w)
a
P = A(x))
((x)
a.
C Identification of Truncation Region
In general, it is difficult to identify the truncation region 2 in Theorem 1 directly because conditioning only on the result of the CP candidate selection, iee, A(a + bz) = A(æ), is intractable. In this case, we have to enumerate all patterns where A(x) appears as a result of simulated anncaling algorithm, which is computationally impractical. To address this issue, we first compute the region conditioned on the process of the algorithm A. That is to say, is guaranteed that the process remains identical within the region. This additional conditioning is often denoted as "over-condittoning" because it is redundant for valid inference. In the case of SI with over-conditee vee--codittonnng, the type error rate still be controlled at the significance level, while the power tends to be low (Lee et al., 2016; Liu et al., 2018; Duy & Takeuchi, 2022). Therefore, we apply an efficient line scarch method based on parametric programming to compute Z where the redundant conditioning is removed for the purpose of improving the power (Duy & Takeuchi, 2022). In the following, we first compute the over-conditioned region in Appendix C.1, and then identify the truncation region Z using parametric programming in Appendix C.2.
C.1 -ver-Conditioning
Since we only need consider onc-dimensional data space in R we define the over-conditioned region goc where the process of the algorithm A remains unchanged as
no aa + bz) = (r E IR (p)(r) SprossA))) S Spre-s(()), (2), = SSA (=))
(21)
where SpP, Spre-SA and SSA are the events characterized by the process of dynamic programming for gen- crating an initial solution, the preliminary experiment to determine the initial temperature, and simulated annealing for making the solution more sophisticated, respectively. This conditioning redundant because goe (a bz) is a subset of the minimum conditioned region Z.
22

---

Published in Transactions on Machine Learning Research (06/2025)
Over-conditioning dynamic programming. We compute fr R SDP(r) SDP(2)} by condi- tioning on all the operations based on the Bellman equation which is used in the dynamic programming algorithm to obtain the optimal solution in (7). Since the Bellman equation consists of the cost C(-) and the penalty B, the condition finally represented by a quadratic incquality, as in the following discussion of simulated annealing. The detail derivation is presenled in Duy al. (2020).
Over-condi simulated annealing containing the preliminary experiment. Considering the algorithm of simulated anncaling, we find that the procodure depends on X = a a++ be at only one specific point, that is, the Metropolis algorithm, which is given by
Acceptance if AE i1 Tu- + c; In(S,)) < 0,
status is Rejection
if AE 1. Ji In(G.)) 2 0,
where AE(T k.,, Tu.- a br) = E(T I.,, brbr S a, br), and T 1, T. respectively represent CP candidates before and after a transition to the neighborhood in the l-th iteration with parameter &,, for the i-th temperature Q. Thus, to compute the region {r R (SsA(r) simulated annealing is identical, need to condition all results of the Mctropolis algorithm at cach lemperature. This condition consists of multiple inequalities as follows:
Sz)(2)) where the process of
{r E R (SA = (z)(
1. I.
AETTT. Tu--1, a + br) + a Im(E,)) < 0 if status is Acceptance
un rER 0
(22)
A a. Tu-1, a + br) e In(G,) 2 0 if status Rejection
where 12 and Lg denote the number of temperature updates and operations at each lemperature for 2, respectively.
We subsequently consider solving this inequality for r. The cost C (F(C) s+1=e) used in the objective function E can be expressed as a quadratic form of X such that
C (p(th =
(d) sym
Σ
p(a))
p(a) 2 st1ie
1=s+1
1 & w (
T X-
( (1oticc @ w & ) T
(d) sym
X
e e
1=++1
T ((d))
x X,
where
((d) = (d)
u l u (d)-T
s4lle E NxN
1=#t1
(d)
1 & w
e-S (1ot1cc & ww E CN
sthet
(d)+ and u U is the complex conjugate of (d) be rewritten follows:
bee
Therefore, the objective function F(T, a br) (5) can
D-1K(t)
D-1
E(T, a + br) (a + br)
Σ
Σ
(dd) (2)
(a + br) + Σ
gld) k(d) +
k + +1:((1)
d=0
k=
d=0
= +eT +er+
(23)
23

---

Published in Transactions on Machine Learning Research (06/2025)
where
D-1 ((t)
62 b
(id)
(c)
d=0
k=
D-1 K +1
((d)
61 2a
R-1 +1:
d=0
K
((() +1
Σ d=0
Bld k(l) t.K.
cu a
+1:
(
a
Thus, the multiple inequalities in (22) can be easily solved after computing the coefficients in (23).
Note that the region fr E R SpreesA) Spr--Sa can be computed similarly to (22) because the preliminary experiment is also based on the Metropolis algorithm. Therefore, this condition is formulated as
fr E R SeSee--AAA)) = Spre-S(())
1x L -u rER =01=1
AE(T to Th a br) + ot im( < 0 if status is Acceptance |E(T T a br) + of ln(S)) 2 status is Rejection
a--1:
where Tit: -5 / , and L in the preliminary experiment correspond the respective parameters in (22). and Tú, T
Based on the above discussion, since the conditioning in (21) is represented as the interscction of multiple quadratic inequalities, the region 2o (a + bz) can be computed by solving them as
R U U
2 + bz)
where, for 2, Loc and Uhe tively, and Rg represents the number of the intervals.
2(r) denote the lower and upper bounds of the r-th over-conditioned region, respec-
C.2 Parametric Programming
Having derived (0 + bz) in the previous analysis, the region Z in Theorem conditioned the result of the algorithm A can be obtained using computational method called parametric programming as follows:
U
20 (a + bz).
(24)
A(a]A(a+bz))aaaa A()==(
The line search method based on (24) for obtaining the region Z required for the computation prelective is detailed in Algorithm 3. Furthermore, the overall procedure for computing selective p-values of the detected CP candidate locations is presented in Algorithm An overview the proposed scarch method is shown in Figure 9.
Watanabe & Suzuki (2021) proposed selective inference method for model selection using simulated an- nealing in latent block models; however, this approach was limited to the specific algorithm and computed an approximated truncation region. In contrast, our proposed method can be applied to not only CP detection in the frequency domain which the subject of this paper, but also a wide range of optimization problems solved using simulated anncaling. Even in such general case, we consider the over-conditioning basod on the process of algorithm as in (21) and can obtain the "exact" truncation region using the parametric programming approach in (24).
24

---

Published in Transactions on Machine Learning Research (06/2025)
Algorithm 3 compute_solution_pee
Input: Time seaquence æ, CP candidates T, CP candidate location T
1: gobs + (x) in (14)
2: Compute a and defined in Theorem
3: Obtain (a + bzhh)) by (21)
S 2 -- 2 bzob)
while S + do
Obtain (a + b2) for E S by (21) S U (a + ba)
8:
if A(a bz) = A(a bz) then
9:
2 Zu 1.UUZ"(A (a (a+ bz)
10: end
11: end while
Output: Truncation region Z
Algorithm SI for Detected CP Candidate Locations
Input: Time seaquence æ
T A+x(z)
Obtain t by (4)
3: for TA E do
4: + compute_sollutonn
T, T,Th) Tk)
5: Compute pk peleotive by (17)
end for
Output: Detected CP candidate locations and the corresponding selective p-values { (Tk, pk
peloctive) )
C
Parametrized line
X + bz
Tist
T
Dynamic programming
co(1)
(2)
(2)
Co
Co CO(3)
Preliminary experiment for initial temperature
Simulated annealing: Adding CP for da time point 57
Reject
Reject
Reject
Reject
Reject
Accept
Accept
Simulated anncaling:
Reject
Accept
Accept
Accept
Reject
Accept
Accept
Moving a CP for dy from time point 74 to 66
Reject
Simulated annealing: The last operation
A(xr)
A(xi
.A(z)
Truncation region: Z = (a IR A(a bz) = A(x))
Figure Schematic illustration of the proposed line search method for the identification of the trunca- tion region. We first compute the over-conditoneed region where the process of the algorithm remains unchanged. Then, identify the truncation region Z by removing the redundant conditioning using para- metric programming.
25

---

Published in Transactions on Machine Learning Research (06/2025)
D Details of Numerical Experiments
D.1 Detailed Descriptions of Comparison Methods
In our experiments, we compared the proposed method (Proposed) with the following methods.
OC: In this method, we consider p-values conditioned on the process of the algorithm A. The p-value is computed by obtaining the over-conditioned region goe in (21) for the observed test statistic (ar). This method is computationally efficient, however, its power is low due to over-conditioning,
OptSeg-SI-oc (Duy et al., 2020): This method uses a p-valuc conditioned only on the event SpP in (21), that is, we consider conditioning on the process of dynamic programming algorithm.
OptSeg-SI (Duy et al., 2020): This method removes over-conditioningg from OptSeg-ST-oc, that is, the p-value is conditioned only on the result of dynamic programming.
Naive: This method uses conventional p-value without conditioning, which is computed
naive = Fo. (Tx(X) 2 (x)))
Pk
Bonferroni:
This is a method to control the type I error rate by applying Bonferroni correction which is widely used as multiple testing correction. Since the number of all possible hypotheses is m = (2D () - 1), the bonferroni p-value is computed by pk ponferroon = min(1, m m.pp..ppeee plavve).
D.2 Computational Time and Computer Resources
We measured the computational time of our proposed method for the synthetic data experiments presented in Section 5.2 and computed the medians for each settings. The results for the type I error rate and power experiments are shown in Figure 10. Panels (a) and (b) indicate that the computational time increases exponentially with the sequence length. In addition, the computational time becomes shorter the signal intensity increases in panels (c) and (d). This may be because the results of hypothesis testing in the case of high intensity are more likely to be obvious, and the inference process be terminated early. All numerical experiments were conducted on a computer with 96-core 3.60GHz CPU and 512GB of memory.
Computational Time (s)
8000
8000
6000
6000
Computational
4000
4000
2000
2000
0
40
60
80
100
40
60
80
100
T
(a) Type error rate (M = 512).
(b) Type error rate (M = 1024).
Computational Time (s]
8000
8000
6000
6000
Computational
4000
4000
2000
2000
0.04
0.08
0.12
0.16
0.04
0.08
0.12
0.16
A
(c) Power (M = 512).
(d) Power (M = 1024).
Figure 10: Computational time in the type I error rate and the power experiments.
26

---

Published in Transactions on Machine Learning Research (06/2025)
D.3 Robustness of Type Error Rate Control
We evaluated the robustness of the Proposed in lerms of the lype I error rale control. For this purpose, we conducted experiments under three distinct noise conditions: (i) unknown noise variance, (ii) non- con-Gaussian noisc, and (iii) correlated noise. The dctails of cach experiment are given below.
Unknown noise variance. We generated 1000 mull sequences in the same manner as the type I error rate experiment with known variance presented in Section 5.2. To estimate the variance o² from the same data, we first applicd CP candidate selection algorithm to identify the segments, and then computed the empirical variance of each segment for all frequencies. Since the estimated variance tended to be smaller than the true value, we adopted the maximum value for each frequency. Given o the average of the values, o could be estimated from the property of DFT as ô = Ar The results for the significance levels a = 0.01, 0.05, 0.1 are shown in Figure 11 and the Proposed still could properly control the type I error ratc.
/M
Non- Gan-Gaussian noise. We considered the case where the noise followed the five non-Gaussian distributions:
skewnorm: Skew normal distribution family.
exponnorm: Exponentially modified normal distribution family.
gennormsteep: Generalized normal distribution lamily whose shape parameter B is limited to be steeper than the normal distribution, ie., 8 2.
gennormflat: Generalized normal distribution family whose shape parameter 3 is limited to be flatter than the normal distribution, i.e, , > 2.
t: Student's distribution family.
To generate sequences used in the experiment, first obtained noise distribution such that the 1- Wasserstein distance from the standard normal distribution (N(0,, 1) was {0.01, 0.02, 0.03, 0.04] in each afore- mentioned distribution family. Subsequently, we standardized the distribution to have mean of 0 and a variance of 1. Then, we generated 1000 null sequences x = (xI,...,NN) where the mean vector was speci- fied in the same manner as described in the type error rate experiment with known variance, and the noise followed the obtained distribution, for T 60. We applyed hypothesis testing using the test statistic with a a=1 for the detected CP candidate locations. The results for the significance levels a 0.05 are shown in Figure 12 and the Proposed could properly control the type error rate for all non-Gaussian distributions.
Correlated noise. We generated 10000 null sequences x (2I,...,NN) ~ M(ss, E), where mean vector was specified in the same manner described in the type error rate experiment with known variance, and covariance matrix Σ was defined as Σ = o² (p)==l) E RNNN with = 1 and & E (0.025, 0.05, 0.075, 0.1), for T 60. After CP candidate selection, we conducted the hypothesis testing using the Lest statistic with 0.01, 0.05, 0.1 are shown in Figure 13. For weak covariance, the Proposed could properly control the type I error rate. However, the type error rate could not be controlled with increasing noise correlation. This remains a challenge for future work.
= 1. The results for the significance levels
D.4 Sensitivity Study for the Penalty Parameter
This section presents a sensitivity study for the penalty parameter within the objective function in (5). As shown in Appendix A, the penalty parameter 3(d) for each frequency d can be theoretically derived using the BIC, leaving the penalty parameter ? to be determined heuristically by users. Thus, we focus on the sensitivity study for 7 in this section. Specifically, we conducted the experiments for the type I error rate (T = 60) and power (A 0.08) on synthetic data as described in Section 5.2, while varying the hyper- parameter in (20) across the values {1.0, 1.5, 2.0, 2.5). The results for the type I error rate and power experiments are shown in Figures 14 and 15, respectively. In Figure 14, the Proposed can control the type error rate at the significance levels = 0.01, 0.05, 0.1 for all values of k. Additionally, Figure 15 shows that the power of the Proposed is the highest of all methods, regardless of the settings of k. In particular, for the experiment with M 512, as the value of к increases, the power of the Proposed tends to be higher. This may be because larger penalty makes distinct spectral changes more likely to be the subject of the hypothesis testing.
27

---

Published in Transactions on Machine Learning Research (06/2025)
a = 0.1
a= 0.1
0.05
0.05
0.01
0.01
0.15
0.15
Type Error Rate
Type |Errrr
0.10
0.10
0.05
0.05
0.01
0.01
0.00
0.00
40
60
80
100
40
60
80
100
(a) M 512.
(b) M 1024.
Figure 11: Robustness of lype error control for estimated variance.
0.20
skewnorm
0.20
skewnorm
exponnorm
exponnorm
gennormsteep gennormflat
gennormsteep
gennormfiat
Type Error Rate
0.10
0.10
Type
0.05
0.05
0.00
0.00
0.02
0.04
0.02
Wasserstein Distance
Wasserstein Distance
(a) M 512.
(b) M = 1024.
Figure 12: Robustness of type T error control for non-Gaussian noise.
0.20
a 0.1
0.20
a= a=0.1
0.05
== 0.05
0.01
0.01
0.15
0.15
Type Error Rate
Type Error
0.10
0.10
0.05
0.05
0.01
0.01
0.00
0.00
0.025
0.05
0.075
0.1
0.025
0.05
0.075
0.1
p
(a) = 512.
(b) M = 1024.
Figure 13: Robustness of type I error control for correlation of noise.
28

---

Published in Transactions on Machine Learning Research (06/2025)
0.20
0.1
0.20
0.1
a 0.05
a 0.01
a= 0.01
0.15
0.15
Type Error Rate
Type Error
0.10
0.05
0.05
0.00
0.00
1.0
1.5
2.0
2.5
1.0
1.5
2.0
(a) M
512.
(b) M 1024.
Figure 14: Sensitivity study of type error rate control for the htyer-paaameee in the penalty parameter ?
1.0
Proposed
1.0
Proposed
OC
OC
Bonferroni
Bonferroni
0.8
0.6
0.6
Power
Power
0.4
0.4
0.2
0.2
0.0
0.0
1.0
1.5
2.0
2.5
1.0
1.5
2.0
2.5
K
(a)
512.
(b) M 1024.
Figure 15: Scnsitivity study of power for the hyper-parameter a in the penalty parameter 7,
D.5 More Results on Real Data Experiments
Additional results for the signals bearings 2, 3, and before and after the spectral intensity enhancements occurred in the BPFO harmonics of bearing are shown in Figures 16, 17, and 18. In panel (a) of each figure, the time variation of frequency spectra where CP candidate locations were falsely detected are shown for the period of indicate that the inferences using p-values of the Proposed and OC are valid. In panel (b), the time variations of the BPFO harmonics where CP candidate locations were correctly detected are presented for the period of 4-6 days (bearing 2), 4.75-6.75 days (bearing 3), and 4-6 days (bearing 4), respectively. In these cases, although the detection was delayed by several days relative to the occurrence of spectral amplification bearing 1, the outer race fault signatures were successfully identified in all bearings using the Proposed.
days when significant spectral changes in bearing 1 did not actually exist. The results
29

---

Published in Transactions on Machine Learning Research (06/2025)
12
Proposed:: 0.353 OC: 0.578
Proposse:: 0.000 Oc: 0.721
10
Amplltde Spectrum
Ampltude Spectum
2-
0.25 0.50 0.75
1.00
1.25 1.50 1.75
2.00 2.25
4,00 4.25 4.50
4.75 5.00 5.25 5.50
5.75 6.00
Time (Day)
(Day)
(a) Inference on a falsely detected CP can- didate location for 3260 Hz (around the 14th harmonic) on 0.25-2.22 days.
(b) Inference on a truely detected CP can- didate location for 1420 Hz (the 6th har- monic) on 4-6 days.
Figure 16: Results of bearing 2. In panel (b), p-values were actually computed by considering CP candidates of the 6th harmonic, 3240 Hz, and 3460 Hz (around the 14th and 15th harmonics).
18
p-values
18
p-values
16-
DC: 0.918 Naiver 1000
16
DC: 0.210
Naive: 0.000
Amplitude Spectrum
14:
Amplitude Spectrum
14
prealues Propessdi 0.000
12
Naive 0,000
10:
6
4-
2
0.75
1.00
1.25
1.50
2.00
4.75
5.00
5.25
5.50
5.75
6.50 6.75
Time (Day)
Time (Day)
(a) Inference a falsely detected CP can- didate location for 3380 Hz (around the 14th harmonic) on 0.25 2.25 days.
(b) Inferences on truely detected CP can- didate locations for 1420 Hz (the 6th har- monic) on 4.75-6.75 6.75 days.
Figure 17: Results bearing
proaluess Proposed:
16
paralues
14.
14.
0.029 Nalve: 0.000
Nalve: 0.000
Ampillude Spectrum
12
Ampitude Spectrum
12
6-
0.25 0.50 0.75
1.00
1.55 1.50 (Day)
1.75 2.00 2.25
4.00
4.25
4.50
4.75
5,00 5.25 5.50
5.75
6.00
(Day)
(a) Inference falscly detected CP can- didate location for 3540 Hz (the 15th har- monic) on 0.25-2.2222. 2.25 days.
(b) Inference on a trucly detected CP can- didate location for 1420 Hz (the 6th har- monic) on 4-6 days.
Figure 18: Results of bearing 4. In panel (b), p-values actually computed by considering not only CP candidate of the 6th harmonic but also CP candidate of 3440 Hz (around the 15th harmonic).
30


Reviews:
[{"review_id": "WgB74OvLBQ", "reviewer": "Reviewer_JRxG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n\n* The paper presents a novel approach by extending statistically significant CP detection using Selective Inference (SI) to the frequency domain, an area that is less explored than the time domain.\n* The proposed method quantifies the statistical significance of detected CPs using $p$-values, which helps in distinguishing real changes from noise.\n* The experimental results on synthetic and real-world data demonstrate the reliability of the proposed method in identifying true CPs while controlling the type I error rate below the specified significance level.\n* The paper examines the robustness of the proposed method to unknown noise variance and non-Gaussian noise, showing promising results under these conditions.\n\n**Weakness:**\n\n* The selection of CP candidates in the first stage relies on a heuristic algorithm (simulated annealing) due to the complexity of the combinatorial optimization problem. While simulated annealing is a practical method for finding approximate solutions, it does not guarantee obtaining the globally optimal set of CPs.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper tackles the problem of change point (CP) detection in the frequency domain to effectively identify anomalies in complex systems.\nThe core contribution of this work is the application of the Selective Inference (SI) framework to the detection of CP in the frequency domain.\nThe proposed SI method quantifies the statistical significance of detected CPs in the frequency domain using $p$-values, ensuring that the detected changes reflect genuine structural shifts in the target system rather than random noise.\n\nThe paper addresses two major technical challenges to achieve this. First, it extends the existing SI framework to the frequency domain by appropriately utilizing the properties of the discrete Fourier transform (DFT). \nSecond, it develops an SI method that provides valid p-values for CPs where changes occur across multiple frequencies. \nThe proposed method consists of two stages: a heuristic selection of CP candidates in the frequency domain, followed by a quantification of the statistical significance of each candidate using the SI framework.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I don't really have changes to request, the points below should rather be seen as potential improvements that the authors could (or not) integrate in a future revision of the paper:\n\n* The authors could discuss the limitations of using a heuristic algorithm (simulated annealing) for CP candidate selection in more detail and potentially explore or suggest alternative or improved heuristic approaches. They could also explain in a few words why the dynamic programming approach that is used for each frequency cannot be extended to the cross-frequency problem.\n* It would be helpful to provide more insight into the impact of hyperparameters $\\beta(d)$ and $\\gamma$ on the results. In this regard, a sensitivity study included in the paper's main body could have been a plus.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "4sttNdEV62", "reviewer": "Reviewer_nMys", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "While the techniques used are known, the proposed multi-resolution CP detection method described in the first half of the paper appears solid and useful.\n\nUnlike the first half, which presents the CP detection approach itself, Section 4, which focuses on the SI approach, is somewhat vague and difficult to follow.\n\nThe procedure seems to be based on the assumption that the observations are i.i.d. with a constant variance $\\sigma^2$ in Eq. (1), the validity of which should be verified. The underlying intuition appears to be that the quadratic statistic in Eq. (11) should follow a Chi-squared-type distribution. However, no clear procedure is provided especially for determining the linear transformation defined in Eq. (17) in a data-driven fashion.\n\nAlthough the authors present a real-world data analysis example in Section 5.3, the details of the analysis are not clearly described, which may leave many readers confused.\n\nOverall, while the first half of the paper presents a useful framework for multi-resolution CP detection, the second half on SI does not seem to add significant value.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a multi-resolution change point (CP) detection approach, which incorporates three main ideas:\n\n1. Performing CP detection simultaneously at multiple frequencies.\n1. Utilizing the dynamic programming (DP) approach for individual CP detection.\n1. Consolidating the individual CPs using a simulated annealing approach.\n\nAdditionally, the authors introduce a method for computing the confidence of the identified CPs using selective inference (SI).", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See the weakness.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "FhKrCBSYLt", "reviewer": "Reviewer_Fed7", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- The title and abstract of the paper are misleading, as the paper is not about anomaly detection in time series, but rather about the detection of a specific type of change points. \n- Without any substantial experience in the field, it is not clear what the contribution of the paper is. The authors point to section 4 multiple times for the novelty, but most results are \"standard approaches\" or results from other papers according to the paper.\n- It is not obvious how the works discussed in the related work section relate to the proposed methodology.\n- Figure 2 is essentially contained in Figure 1.\n- The writing is inconsistent and the notation is overly complicated, making the paper very difficult to follow.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a statistical test for the significance of mean-shift change points in the frequency domain.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Significant changes:\n- The title and text should be changed to avoid confusion with time-series anomaly detection. If anomaly detection is an intended application, their relation to change points should be clearly discussed.\n- Figure 2 should probably be removed.\n- The authors should clearly point out and discuss the contributions of the paper, and improve the comparison to related work.\n\nGeneral Writing improvements:\n- Some equations that should be numbered are not and some that are should not be, e.g. Why does eq. (4) have a dedicated equation? It seems to only be used for the definition of $K$.\n- In general there is too much empty space caused by dedicated equations, which should probably be inline, e.g. (18) and (19).\n- I am not familiar with the notation used in the definition of $c_{sym}^{(d)}$ and $A_n^{(d)}$. Usually, a case is followed by a statement that evaluates to true or false.\n- The algorithms should be streamlined, e.g.\n    - The nested conditions in Algorithm 1 can be combined without loosing expressiveness.\n    - Algorithm 2 uses a mix of precise and imprecise notation, e.g. What is the \"desired value\"? Is \"d\" also uniformly sampled from $D^{init}$? Are the operations uniformly sampled? etc.\n- The indices of $\\tau$'s are hard to read. It seems, for every $\\tau$ the authors only need its predecessor, its successor, and its associated frequency. If this is indeed the case, the equations could be simplified, e.g. by using $\\tau_d^{pred}$, $\\tau_d$, $\\tau_d^{suc}$.\n- The definitions of the two means in eq. (11) should be generalized for an arbitrary range $\\bar{F}^{(d)}_{i:j}$ to remove redundancy.\n- There are missing citations, e.g. \"This additional conditioning on Q(X) is a standard approach for computational tractability in the SI literature.\"\n- The cumulative distribution function shares a name with the Fourier Transform, which can be easily overlooked, since they both also come with multiple indices each.\n\nFurther points which would benefit from more explanation:\n- What is the motivation for $\\sigma^{-1}$ and $a_{k^{(d)}}$?\n- It is not obvious, why eq. (12) holds. \n- The step \"[...] $Z = T_k(X)$ [...] that satisfy $X = a + bZ$ [...]\" is not clear.\n- Theorem 2.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This work presents a novel statistical testing method for change-point detection in the frequency domain, which is supported by experiments on toy and real data. The authors have revised the paper, especially for readability, according to the reviewers' feedback and majority are positive about the paper. The revised paper clearly meets the criteria for acceptance. The authors can consider further improving readability for a more general machine learning audience.
Decision: Accept (poster)