How do I fit two function in a same data at a time depending on the range ?

I have a data set consists of 1000 points . It follows y=t^2 upto first few data points and after that it follows y=-t. How do I find out the point of deviation (tc) from y=t^2 to y=t just by fitting the function

y=x^2 if t<tc

and y= tc^2+(tc-t) if t>tc ?

[Walter Roberson]

Is there a particular reason to do that, instead of looking at the numeric derivative of the function first to estimate where the inflection point is? gradient(t,y) would be positive until right near tc, then turn negative so you can look for the change in sign of the gradient to bracket where tc must be.

Answer 1

[Torsten]

t = load("t.mat");
t = t.t;
y = load("y.mat");
y = y.y;

value_min = Inf;
for i = 1:numel(t)
   fun1 = @(time)time.^2;
   fun2 = @(time)t(i)^2+(t(i)-time);
   value1 = sum((fun1(t(1:i))-y(1:i)).^2);
   value2 = sum((fun2(t(i+1:end))-y(i+1:end)).^2);
   value = value1 + value2;
   if value < value_min
     value_min = value;
     i_min = i;
   end
end
t(i_min)

ans = 0.4000

hold on
plot(t,y)
plot(t(i_min),y(i_min),'o')
hold off

[Sayanta Goswami]

Thanks for the respons. I want to use cfit tool and visualize the fitting .i mean, I need to define the two functions and find out tc from there where tc will be a fitting parameter. How to do that ?

Answer 2

[the cyclist]

If the data are truly smooth, as what you posted (and does not have statistical fluctuation), then you just need to find the first point where y is smaller than the previous value. Here is one way.

load("t.mat","t")
load("y.mat","y")

tc = t(find(diff(y)<0,1))

tc = 0.4000

Even, easier, you could find the peak location using max:

[~,ymax_ind] = max(y);

tc = t(ymax_ind)

tc = 0.4000

If you truly need to fit statistical data, you could try fitnlm, but I'm not sure how well it will do on the discontinuity. (I'd have to try to it, but it's not worth doing if the above is good enough for you.)

[Walter Roberson]

Not exactly -- tc might happen to be in-between two of the input t values, so once you figured out where the dividing point is, you would want to fit the exact tc value from a few samples after the boundary.

[the cyclist]

Agreed. I neglected to mention the nuances of needing the point just before, or the point just after, or an estimate of the "true" peak. Depends on how accurate one needs to be.

Answer 3

[Walter Roberson]

Provided that tc is not at or beyond the end of the data, you can calculate

tc = sqrt(4*t(end) + 4*y(end) + 1)/2 - 1/2;

This would be subject to noise in the y or t measurements, so you might want to calculate over several points and take the mean()

Answer 4

[the cyclist]

Here is a solution in which I fit the non-linear function (of the form you prescribe) to the data.

(For fun, I added some statistical noise to y. You can take that out.)

Note that the model fit finds the value of TC quite accurately. (It would be perfect, if I had not added the noise.)

Interestingly, the model fit fails completely if I set the starting guess to be greater than 0.4. I did not do a deep dive into why, but be aware that non-linear fits are often sensitive to the initial parameter guesses.

% Load the data
load("t.mat","t")
load("y.mat","y")

% Convert to column vectors
t = t';
y = y';

% Add a touch of noise to y
y = y + 0.005*randn(length(y),1);

% Tabulate the data
tbl = table(t,y);

% %  Fit the model

% Define function that will be used to fit data
% (F is a vector of fitting parameters)
f = @(F,t) (t < F(1)).*t.^2 + (t > F(1)).*(F(1).^2 - (t-F(1)));
beta0 = 0;
mdl = fitnlm(tbl,f,beta0)

mdl = 
Nonlinear regression model:
    y ~ (t<F1)*t^2 + (t>F1)*(F1^2 - (t - F1))

Estimated Coefficients:
          Estimate        SE        tStat     pValue
          ________    __________    ______    ______

    F1    0.39993     0.00011082    3608.7      0   


Number of observations: 1001, Error degrees of freedom: 1000
Root Mean Squared Error: 0.00489
R-Squared: 0.999,  Adjusted R-Squared 0.999
F-statistic vs. zero model: 1.33e+06, p-value = 0

% Calculate the model values at the empirical t
y_predicted = predict(mdl,t);

% Plot the data and fit
figure
plot(t,y,'.',t,y_predicted,'r');
legend('data','fit')