Incorporating local regression and smoothing splines to smooth data, can this be achieved by one of them?

Hi, we are trying to build a curve-fitting model according to some constraints, and we need a bit of advice here, to see whether we are on the right track or not. The model uses two different methods of smoothing. The first one is local regression using the function

smoothdata(y,"loess",20);

The results of this function will smoothen the curve for us and would allow us to achieve certain constraints, mainly nonlinear shifting of some data points according to certain conditions, but it will not generate a completely smooth curve, there will be some parts of the generated curve that needs further smoothing. For that, we used the second method which is smoothing splines using the curve fitting toolbox.

Now, we are afraid that using two smoothing methods might be unnecessary, and we are not aware if what we want can be achieved using one method, so the question is, would there be a way to use one of these two methods to fulfill the constraints that we have and get a full smoothed curve at the end, or both methods has to be incorporated.

this is the code for the two methods and a sample of the data is attached

            [myb,ib] = max(yb2);
            mxb = xb2(ib);

            xs = xb2;              
            ys = smoothdata(yb2,"loess",20);

            if max(ys) > ys(ib)
                ys(ib) = max(ys);
            end
            % the goal is to make the peak of the fitted curve as close as
            % it can be to the the two data points  hmf2_RO and NmF2_RO

            % non linear y  shift, trying to create a y shift based on x distance from peak point (exp
            % distribution)
            dx = abs(xs-hmf2_RO);
            dy = abs(ys(ib)-NmF2_RO);
            y_shift = dy*exp(-dx.^2/3000);
            yss = ys + y_shift;

            [fitrslt, gof1] = createFit2(xs, yss);
             h1=plot(fitrslt,xs,yss);
            curveHandle = findobj(h1,'DisplayName', 'fitted curve'); 
            X = curveHandle.XData; 
            Y = curveHandle.YData;
            
            function [fitresult, gof] = createFit2(xs, yss)
            [xData, yData] = prepareCurveData( xs, yss );
            
            % Set up fittype and options.
            ft = fittype( 'smoothingspline' );
            
            % Fit model to data.
            [fitresult, gof] = fit( xData, yData, ft );
            
            % Plot fit with data.
            figure( 'Name', 'untitled fit 1' );
            h = plot( fitresult, xData, yData );
            legend( h, 'yss vs. xs', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
            % Label axes
            xlabel( 'xs', 'Interpreter', 'none' );
            ylabel( 'yss', 'Interpreter', 'none' );
            grid on
            
            end

[Mathieu NOE]

hello

I played with your data (and had to use some available alternatives to spline smoothing as I don't have the required toolbox) but e final goal and what you want to do is still unclear to me

maybe you should tell us a bit more about what is the end goal and then we can come up with a solution.

clearvars, close all;

load('data.mat')
x = xs;
y = yb2;
clear xs yb2

id = x<700;
x=x(id);
y=y(id);

% smoothdata realisation
ys = smoothdata(y,"loess",20);

% B spline smoothing 
m = 8;% B-Spline Degree
s = f_Bspline([x';y'], m, 0.5);% B-Spline
xspl = s(1,:);
yspl = s(2,:);

% curve fit using fminsearch
f = @(a,b,c,x) a.*exp(-(x-b).^2 / c.^2);
obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);

sol = fminsearch(obj_fun, [max(y),max(x)/2,max(x)/6]);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
xx = linspace(min(x),max(x),300);
y_fit = f(a_sol, b_sol,c_sol, xx);
yy = interp1(x,y, xx);
Rsquared = my_Rsquared_coeff(yy,y_fit); % correlation coefficient

plot(x,y, 'r',x,ys, 'k',xspl,yspl, 'g',xx, y_fit, '-' , 'Linewidth', 2.5)

title(['Gaussian Fit / R² = ' num2str(Rsquared) ], 'FontSize', 5)
ylabel('Intensity (arb. unit)', 'FontSize', 14)
xlabel('x(nm)', 'FontSize', 14)
legend('raw data','smoothdata','B spline smooth','gaussian fit')

%%%%%%%%%%%%%%%%%%%%%%%%%
function s = f_Bspline(C, m, step)
% Calculates a B-Spline of degree m using the control points C.
% 
% C: 2-dimensional control points (x_0, ... , x_n; y_0, ... , y_n)
% m: B-Spline degree
% step: Step size of t.
%
% s: B-Spline. s(1,:) -> x component, s(2,:) -> y component

%% Parameters

% Control point's x and y components
x = C(1,:);
y = C(2,:);

% Number of control point - 1
n = size(x,2) - 1;

% Knot vector
T = f_BSpline_KnotVector(m,n);

% B-Spline intervall
t = 0:step:(n-m+1);

%% Calculate B-Spline

for z=1:1:size(t,2)
    ti = t(z);
    s(1,z) = 0;
    s(2,z) = 0;
    for i=0:1:n
        % Base B-Spline
        B = f_BSpline_BaseSpline(i, m, ti, T);
        % x component
        s(1,z) = s(1,z) + x(i+1) * B;
        % y component
        s(2,z) = s(2,z) + y(i+1) * B;
    end
end

end

function T = f_BSpline_KnotVector(m,n)
% Calculate knot sequence.
% m: Degree of B-Spline
% n: Number of control points - 1
%
% T = [t0, ... , t_n+m+1] : Knot sequence / vector

T = zeros(1, (n+m+2));
for j=0:1:(n+m+1)
    if j <= m
        Ti = 0;
    elseif j >= (m+1) && j <= n
        Ti = j - m;
    elseif j > n
        Ti = n - m + 1;
    end
    T(j+1) = Ti;
end

end

function B = f_BSpline_BaseSpline(i, k, t, T)
% Calculate Base B-Spline B_i,k
% k: B-Spline degree
% T: Knot sequence
% t: Current t parameter
%
% B: Base B-Spline at t.

% Index shift
j = i + 1;

if k == 0
    % Corrected end of recursion
    if (t >= T(j) && t < T(j+1))
        B = 1;
    % Exception: the last basis function is equal to 1 also at the last knot.
    % There might be a more elegant way of writing it; I am no Matlab expert.
    elseif ( T(j+1) == T(end) && t == T(end))
        B = 1;
    else
        B = 0;
    end
else
    % Check dividing by zero
    if T(j+k) ~= T(j)
        A = (t -T(j))/(T(j+k) - T(j));
    else
        A = 0;
    end
    if T(j+k+1) ~= T(j+1)
        B = (T(j+k+1) - t) / (T(j+k+1) - T(j+1));
    else
        B = 0;
    end
    % Calculate base B-Spline
    B1 = f_BSpline_BaseSpline(i,   k-1, t, T);
    B2 = f_BSpline_BaseSpline(i+1, k-1, t, T);
    B = A * B1 + B * B2;
end


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

[Malak]

Thank you for your replay. I can see that Gaussian fit is givving you a very smooth and nice curve. However, our ultimate goal is not just to fit the curve that we have.The curve that we have includes some coruppted data and we are aiming to correct this data so for each curve we have a point (nmf2, hmf2) the ones that I am attaching with the data, and w want to make the peak of our original profile as colse as it can be to this point or even be exactly on. Another issue is that, we do not know where exactly the peak of our original curve is, because again the data there are corrupted and we can't assume that the heightest value is the peak, the peak should be some how in the midile of the profile, (just as you are getting with the Gaussian fit). SO for that purpuse, I used smoothdat first to be able to some how see where the peak can be, then I applied some non liear shifftng to shift the peak point to match the point (nmf2, hmf2). but the resultant profile was not smooth, so I used smoothing splines again to obtain a smoother curve. Now can this be done in less steps and using one approach maybe and not two or is there an easier way to do this?

This is what we get when we try both approaches together

[Mathieu NOE]

hello again

a modified gaussian fit that goes exactly through the peak point

not other smoothing required ;

does it answer your need ?

load('data.mat')
x = xs;
y = yb2;
clear xs yb2
id = x<700;
x=x(id);
y=y(id);
% curve fit using fminsearch
f = @(a,c,x) a+ (NmF2_RO-a).*exp(-(x-hmf2_RO).^2 / c.^2);
obj_fun = @(params) norm(f(params(1),params(2),x)-y);
sol = fminsearch(obj_fun, [y(end) max(x)/10]);
a_sol = sol(1);
c_sol = sol(2);
xx = linspace(min(x),max(x),300);
y_fit = f(a_sol,c_sol, xx);
yy = interp1(x,y, xx);
Rsquared = my_Rsquared_coeff(yy,y_fit); % correlation coefficient

plot(hmf2_RO,NmF2_RO,'dk',x,y, 'r',xx, y_fit, 'b-' , 'Linewidth', 2.5)
title(['Gaussian Fit / R² = ' num2str(Rsquared) ], 'FontSize', 15)
ylabel('Intensity (arb. unit)', 'FontSize', 14)
xlabel('x(nm)', 'FontSize', 14)
legend('peak data','raw data','gaussian fit')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end