Trying to adjust weighting on PART of a sigmoid fit

I included a picture of the current fit that I am getting, what is in black is the actual data. So far the fit is pretty good but my PI would like it to fit a bit better and he suggested weighting the first part of the fit less than the later parts since we really only care from about 500 seconds onwards. Any suggestions on how to do this? I can include the actual matlab file for the fit if that would be helpful for anyone

Thanks!

Answer 1

[Image Analyst]

Then simply don't include that data in your training data set. Try

firstIndex = find(x > 500);
partialYData = y(firstIndex : end); % Extract out only y data where x is more than 500.
% Continue with however you did the fit, like with fitnlm or whatever (Demo attached).

Then just do whatever you did with partialYData instead of the full y data.

Answer 2

[Sam Chak]

Hi @ROBERTA

The data points do not exhibit a symmetrical sigmoidal trend, suggesting that the commonly used logistic function may not be suitable. You could consider fitting the data with asymmetrical sigmoidal models such as the Generalized logistic function, the Gompertz function, or the cumulative Gumbel distribution function.

Upon initial observation (based on experience), it appears that a 3rd-order dynamically-damped system with time delay could be a reasonably good fit. This type of trend also bears resemblance to the popular First-Order-Plus-Dead-Time (FOPDT) Model.

w   = 0.01;     % <-- tune this value only
amp = 6.5;      % amplitude (measurable)
Td  = 300;      % time delay (measurable)
G   = amp*tf(w^3, [1 3*w 3*w^2 w^3], 'InputDelay', Td)

G =
 
                             6.5e-06
  exp(-300*s) * ---------------------------------
                s^3 + 0.03 s^2 + 0.0003 s + 1e-06
 
Continuous-time transfer function.

plotoptions      = timeoptions('cstprefs');
plotoptions.Grid = 'on';
stepplot(G, 4000, plotoptions)

[Image Analyst]

She said she was fine with ignoring all data below time=500, so why did you include it?

[Sam Chak]

Because I can 'see' what she probably cannot, I might have a reasonably good guess of what could be a suitable fit. It's important to note that the measured data points likely originate from a dynamic system, which may lack a true analytical solution. In certain industrial processes, introducing time delays is logical when liquids or gases have to traverse long distances from the upstream to downstream areas.

From the image, as the response approaches the steady-state (around 6.5), the measurements exhibit some fluctuation accompanied by noise. This phenomenon is likely due to the placement of the sensor device near the steady-state level. It further strengthens my belief that the measurements are derived from a dynamic system rather than discrete events.

[ROBERTA]

Thank you so much for the detailed response! I sort of did a combination of filtering out a bit of data as well as using the Gompertz fit

This is probably a super dumb math question, but it looks like the fit extends past 0, I've tried to end the fit at 0 and then it messes up the whole thing. Is there a way to bound it without compromising the rest of the fit?

[Image Analyst]

I'm not familiar with the plot functions he used. But will you actually try to use the estimated formula to estimate values down in the region where the values are all essentially zero?

[Sam Chak]

Hi @ROBERTA

The coefficient 'd' in the Gompertz function represents the offset distance that shifts the function along the y-axis. The standard Gompertz function provided in the library models already includes this offset. If you wish to exclude the offset effect from the Gompertz function, you can either fix 'd = 0' at both lower and upper bounds, or manually remove the 'd' term from the equation. See Example below:

If you would like us to assist you in finding a better fit, please attach the data sheet (by clicking the paperclip icon). To save the variables 'xData' and 'yData' from the workspace to a mat-file, use the following syntax:

save myData.mat

Example:

%% data generation
a   = 0.9911;
b   = 0.0058;
c   = 426.7144;
d   =-0.1998;
f   = @(x) (d + (a - d)*exp(-exp(-b*(x - c))));
x0  = fzero(f, 320);
x   = x0:0.1:3500;
y   = (d + (a - d)*exp(-exp(-b*(x - c))));              % Gompertz fcn with offset d

%% fitting preparation using Gompertz function
[xData, yData] = prepareCurveData(x, y);
fo  = fitoptions('Method', 'NonlinearLeastSquares', ...
                 'Lower', [0.9, 0.005, 450], ...        % lower bounds of {a, b, c}
                 'Upper', [1.0, 0.009, 500]);           % upper bounds of {a, b, c}

ft  = fittype('a*exp(-exp(-b*(x - c)))', 'options', fo);    % Gompertz fcn with d = 0

%% fit the data 
a0          = 0.95;      % initial guess
b0          = 0.007;
c0          = 475;
[yfit, gof] = fit(xData, yData, ft, 'StartPoint', [a0, b0, c0])

yfit = 
     General model:
     yfit(x) = a*exp(-exp(-b*(x - c)))
     Coefficients (with 95% confidence bounds):
       a =      0.9897  (0.9896, 0.9898)
       b =    0.006504  (0.006497, 0.00651)
       c =       483.3  (483.2, 483.4)

gof = struct with fields:
           sse: 1.5262
       rsquare: 0.9988
           dfe: 31729
    adjrsquare: 0.9988
          rmse: 0.0069

%% Plot result
plot(x(:,1:80:end),  y(:,1:80:end), '.', 'markersize', 4), hold on
plot(yfit), grid on
grid on, xlabel('t'), ylabel('y')
legend('Data', 'Fitted Gompertz')

[S0852306]

@ROBERTA,

Most data measured from dynamics system do not have analytic solutions. so simple explicit model (such as Gompertz or logistic function) probably won't work.

There are 2 approaches,

1. System identification: instead of finding the explicit function , (θ are the unknown parameters)find x implicitly by fitting the state function f, i.e. , θ could be found by dynamic programing (system identification toolbox is able to do this). once f is known you can evaluate x using ode solvers, this approach perform well when the "structure" of the dynamics system is well known, e.g. single pendulum: (find a and b).
2. Universal approximator: using some complex black box model to fit the data; if your purpose is just to get a good fit this method usually works. however, the fitted black box function has no physical meaning. here's some results I got:

The underlying model is a neural net which is suitable to perform highly complex function approximation. the maximu absolute error (inf-norm) is about 0.02.

code: FEX

%% Generate data for fitting, follow Sam Chak's code.
clear; clc;
a   = 0.9911;
b   = 0.0058;
c   = 426.7144;
d   =-0.1998;
f   = @(x) (d + (a - d)*exp(-exp(-b*(x - c))));
x0  = fzero(f, 320);
x   = x0:0.1:3500;
y   = (d + (a - d)*exp(-exp(-b*(x - c))));  

x0 = linspace(0, 326, 10000);
y0 = zeros(1, 10000);

x1 = [x0 x];
y1 = [y0 y];
y1 = smooth(y1, 100);
y1 = y1.';
%% Initialize the model.
LayerStruct = [1, 8, 8, 8, 16, 1];
NN.InputAutoScaling = 'on';
NN.LabelAutoScaling = 'on';
NN.ActivationFunction = 'Sigmoid';
NN=Initialization(LayerStruct,NN);
%% Perform least squares.
option.MaxIteration = 300;
NN = OptimizationSolver(x1, y1, NN, option);
%% Plot the results.
Prediction = NN.Evaluate(x1);
scatter(x1, y1, 3)
hold on
plot(x1, Prediction, 'LineWidth', 2)
legend('Fitting', 'Data')
%% Stats and error analysis
% R = FittingReport(x1, y1, NN);