✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
Smoothing splines are standard methods of nonparametric regression for obtaining smooth functions from noisy observations. But as splines are twice differentiable by construction, they cannot capture potential discontinuities in the underlying signal. The smoothing spline model can be augmented such that discontinuities at a priori unknown locations are incorporated. The augmented model results in a minimization problem with respect to discontinuity locations. The minimizing solution is a cubic smoothing spline with discontinuities (CSSD) which may serve as function estimator for discontinuous signals, as a changepoint detector, and as a tool for exploratory data analysis. However, these applications are hardly accessible at the moment because there is no efficient algorithm for computing a CSSD. In this paper, we propose an efficient algorithm that computes a global minimizer of the underlying problem. Its worst case complexity is quadratic in the number of data points. If the number of detected discontinuities scales linearly with the signal length, we observe linear growth of the runtime. By the proposed algorithm, a CSSD can be computed in reasonable time on standard hardware. Furthermore, we implement a strategy for automatic selection of the hyperparameters. Numerical examples demonstrate the applicability of a CSSD for the tasks mentioned above.
⛄ 部分代码
function output = cssd_cv(x, y, folds_cell, delta, startingPoint, varargin)%CSSD_CV K-Fold cross validation for cubic smoothing spline with discontinuities %% cssd_cv(x, y, folds_cell, delt, xx, delta) automatically determines% the model parameters p and gamma of a CSSD model based on minimization of % a K-fold cross validation score. % % Note: The minimization process uses standard derivative free% optimizers (simulated annealing and Nelder-Mead simplex downhill). So it% is not guaranteed that the result is a global minimum. % Trying different starting points often helps if the results with the default starting point % are unsatisfactory.%% Input% x: vector of data sites%% y: vector of same lenght as x or matrix where y(:,i) is a data vector at site x(i)%% folds_cell (optional): a cell array of index vectors corresponding to the% K folds (if [] the 5 folds are chosen randomly)% % startingPoint (optional): Starting point startingPoint = [p; gamma] for the optimizer% (default is [0.5; 1])%% delta: (optional) weights of the data sites. delta may be thought of as the% standard deviation of the at site x_i. Should have the same size as x.% - Note: The Matlab built in spline function csaps uses a different weight% convention (w). delta is related to Matlab's w by w = 1./delta.^2% - Note for vector-valued data: Weights are assumed to be identical over% vector-components. (Componentwise weights might be supported in a future version.)%% Output% output = cssd_cv(...)% output.p: % p: parameter between 0 and 1 that weights the rougness penalty% (high values result in smoother curves)% output.gamma: parameter between 0 and Infinity that weights the discontiuity% penalty (high values result in less discontinuities, gamma = Inf corresponds to a classical smoothing spline% output.cv_score: best cv_score achieved by the optimizer% output.cv_fun: a function handle to the CV scoring function%% See also CSSDif nargin < 3 folds_cell = {};endif nargin < 4 delta = [];endif nargin < 5 startingPoint = [];endp = inputParser;addOptional(p,'verbose', 0);addOptional(p,'maxTime', Inf);parse(p,varargin{:});verbose = p.Results.verbose;maxTime = p.Results.maxTime;if isempty(delta) delta = ones(size(x));end% checks arguments and creates column vectors (chckxywp is Matlab built in)% Matlab uses the parameter w which is related to delta of De Boor's book by w = 1./delta.^2w = 1./delta.^2;[xi,yi,~,wi] = chckxywp(x,y,2,w,0);deltai = sqrt(1./wi);[N,D] = size(yi);% create folds if not givenif isempty(folds_cell) folds_cell = kfoldcv_split(N, 5);endgamma_tf = @(b) atan(b) * 2/pi;gamma_itf = @(a) (tan(a * pi/2));% generate cv score functioncv_fun = @(p, gamma) cssd_cvscore(xi, yi, p, gamma, deltai, folds_cell);cv_fun_vec = @(z) cv_fun(z(1), gamma_itf(z(2))); % vectorised and transformed version for optimization% perform optimizationsaoptions = {@simulannealbnd,'MaxTime', maxTime};options = {};if verbose saoptions = {saoptions{:}, 'Display','iter','PlotFcns', {@saplotbestx,@saplotbestf,@saplotx,@saplotf}}; options = {options{:}, 'Display','iter', 'PlotFcns', {@optimplotx, @optimplotfunccount, @optimplotfval}};end% set starting point of optimizersif isempty(startingPoint) startingPoint = [0.5; 1]; % starting values: p =0.5, gamma = 1end% transform starting point to lie in the unit square [0,1]^2startingPoint_tf = [startingPoint(1); gamma_tf(startingPoint(2))];% invoke chain of standard derivative-free optimizers on [0,1]^2: simulated% annealing followed by Nelder-Mead simplex downhill[improvedPoint_tf, cv_score] = simulannealbnd(cv_fun_vec, startingPoint_tf, [0;0], [1;1], optimoptions(saoptions{:})); % simulated annealing[improvedPoint_tf, cv_score] = fminsearch(cv_fun_vec, improvedPoint_tf, optimset(options{:})); % refine result of simulated annealing by Nelder Mead downhill% improved p and gammap = improvedPoint_tf(1);gamma = gamma_itf(improvedPoint_tf(2)); % transform gamma back to [0, Inf)% set ouptputoutput.p = p;output.gamma = gamma;output.cv_score = cv_score;output.cv_fun = cv_fun;end
⛄ 运行结果
⛄ 参考文献
M. Storath, A. Weinmann, "Smoothing splines for discontinuous signals", arXiv:2211.12785, 2022
