Fitting a monotonically increasing spline function

I want to fit a monotonously increasing smooth spline function for a dataset

x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]

x = 1×17

0 0.7500 1.8000 2.2500 3.7500 4.5000 6.4500 6.7500 7.5000 8.3250 10.8750 11.2500 12.5250 12.7500 15.0000 20.8500 21.0000

y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,
    3.56250726, 3.6209845,  3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

The current fit using interp1d looks like the above. I would like to know how to fit a monotonously increasing spline function.

Answer 1

[Bruno Luong]

One way is to use my https://uk.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation

x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,3.56250726, 3.6209845,  3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]

nknots=10;
opt=struct('shape',struct('p',1,'lo',zeros(1,nknots),'up',inf(1,nknots)));
pp=BSFK(x,y,4,nknots,[],opt); %FEX

xi=linspace(min(x),max(x),100);
yi=ppval(pp,xi);
plot(xi,yi,'-',x,y,'or')

[Deepa Maheshvare]

@Bruno Luong Thanks so much. This is great.

Could you please explain a bit about how nknots=10 is set?

Also, may I know if there is a python implementation of this function.

[Bruno Luong]

I select 10 knots because you have 13 points, do 10 should be enough. I believe at tte end BSFK remove the knots and keep relevent 4 or 5.

No there is no equivalent python implementation of my function.

[Deepa Maheshvare]

Thank you for the clarification.

[Deepa Maheshvare]

@Bruno Luong, Could we get the polynomial function of the fitted curve?

[Bruno Luong]

You could use the same technique I gave here :

https://uk.mathworks.com/matlabcentral/answers/1794985-how-to-constrain-the-resulting-equation-from-a-polynomial-surface-fit-to-a-positive-range?s_tid=srchtitle

instead of imposing constraints of f(xi)>=0, you imposes df/dx(xi) >= 0 for some xi "sufficiently" dense.

[Deepa Maheshvare]

@Bruno Luong, Thank you. I am not sure if I understand how to get k, l and C.

% cofficients of 2D polynomial 3d order
k = [0 0 1 0 1 2 0 1 2 3];
l = [0 1 0 2 1 0 3 2 1 0];
C = [xn.^k.*yn.^l]; % please no comment about my use of bracket here
% d = z;

Also could you please help with including

"instead of imposing constraints of f(xi)>=0, you imposes df/dx(xi) >= 0 for some xi "sufficiently" dense."

in the below code?

% ref: https://in.mathworks.com/matlabcentral/answers/1794985-how-to-constrain-the-resulting-equation-from-a-polynomial-surface-fit-to-a-positive-range#answer_1042065
x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,3.56250726, 3.6209845,  3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]
% Stuff needed to normalize the data for better inversion
[xmin, xmax] = bounds(x);
[ymin, ymax] = bounds(y);
xynfun = @(x,y)deal((x(:)-xmin)/(xmax-xmin),(y(:)-ymin)/(ymax-ymin));
[xn,yn]=xynfun(x,y);


% cofficients of 2D polynomial 3d order
k = [0 0 1 0 1 2 0 1 2 3];
l = [0 1 0 2 1 0 3 2 1 0];
C = [xn.^k.*yn.^l]; % please no comment about my use of bracket here
% d = z;

% Constraint positive of 3 x 3 points in the recatagular domain to be positive,
% it should be enough
[XNC,YNC] = meshgrid(linspace(0,1,3),linspace(0,1,3));
A = -[XNC(:).^k.*YNC(:).^l]; % please no comment ...
b = 0+zeros(size(A,1),1); % A*P<=0 means polynomial at (xnc,ync)>=0
P = lsqlin(C,d,A,b);

[Bruno Luong]

Monotonic polynomial

x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.];
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,3.56250726, 3.6209845,  3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498];

% Stuff needed to normalize the data for better inversion
[xmin, xmax] = bounds(x);
xnfun = @(x)(x(:)-xmin)/(xmax-xmin);

xn=xnfun(x);
% cofficients of 2D polynomial 3d order
k = 0:7;
C = [xn.^k]; % please no comment about my use of bracket here
d = y;
% Constraint positive of 3 x 3 points in the recatagular domain to be positive,
% it should be enough
XNC = linspace(0,1,41);
A = -[k.*XNC(:).^(k-1)]; % please no comment ...
A(:,k==0)=0;
b = 0+zeros(size(A,1),1); % A*P<=0 means polynomial at (xnc,ync)>=0
P = lsqlin(C,d,A,b);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

% Graphical check
% Create a grided model surface
xi=linspace(xmin,xmax,201);
Xin=xnfun(xi);
Yi=[Xin.^k]*P;  % please no comment about my use of bracket here

close all
plot(xi,Yi);
hold on
plot(x,y,'or')
xlabel('x')
ylabel('y')