Using cumulative areas as an error threshold

Hello!

I have a function r(t), which in this case is a logistic function, but for all intents and purposes could be any function r(t).

My goal is to begin at the first value of r and compare it to the second value of r. If the area formed by the triangle under the curve is above a certain threshold then I keep the second value of r. If it is below, then I move to the third value of r and compare the area of the triangle now formed between the first and third r values to the threshold. I repeat until I find the value of r that is above the threshold. For the parameters below, this r value is the 75th index.

I then begin the process again between the 75th and 76th index and so on and so on. My final result is therefore a logical array that is true when r is above the threshold.

I have written a looping procedure that carries this out successfully, but my instinct is telling me there should be a more efficient solution, perhaps involving some loop-less cumulative function. Or a clever use of a step function...

Would appreciate any suggestions / ideas!

r0 = 3; % start r
r1 = 6; % end r
t0 = 5; % halfway point
k = 2; % exponential growth
r_ = @(t) 1./(1+exp(-k*(t-t0))) * (r1-r0) + r0;
t = linspace(0,15,301);
r = r_(t);

dt = mean(diff(t));
dr = abs(diff(r));

threshold = 5e-3;

A_ = 0;
lgc = false;
for idr = 1:numel(dr)
  A = dt*(dr(idr))/2;
  A_ = A_ + A;
  if A_ >= threshold
    lgc(idr) = true;
    A_ = 0;
  else
    lgc(idr) = false;
  end
end
lgc = [true lgc];

hold on
plot(t,r)
plot(t(lgc),r(lgc),'x')

Answer 1

[Voss]

You can use cumsum to perform the cumulative sum over all dr at once. Then loop over sections of that cumulative sum, finding the first index where it exceeds threshold each time and subtracting off that element. You still have a loop, but it has fewer iterations.

r0 = 3; % start r
r1 = 6; % end r
t0 = 5; % halfway point
k = 2; % exponential growth
r_ = @(t) 1./(1+exp(-k*(t-t0))) * (r1-r0) + r0;
t = linspace(0,15,301);
r = r_(t);

dt = mean(diff(t));
dr = abs(diff(r));

threshold = 5e-3;

lgc = [true false(1,numel(dr))];
A = dt*cumsum(dr)/2;
while true
    idx = find(A >= threshold, 1);
    if isempty(idx)
        break
    end
    lgc(idx+1) = true;
    A = A-A(idx);
end

hold on
plot(t,r)
plot(t(lgc),r(lgc),'x')

[Jack]

Thank you. Your code is definitely cleaner!

I've spotted in both mine and your code that the area of each consecutive triangle is used for the threshold, whereas I think I am actually looking for the area from the last index above the threshold.

For example,

If |r(2)-r(1)|*dt/2 < Threshold, then the next check should be |r(3)-r(1)|*(2*dt)/2, and not ( |r(2)-r(1)|*dt/2 + |r(3)-r(2)|*dt/2 ). The 2nd aproach would always underestimate the area of the triangle.

I hope that is somewhat clear.

Do you think this is possible using cumsum?

[Voss]

Not with cumsum, but you can do that using the expression for A below.

r-r.' is the difference between each element of r and every other element of r (including itself), giving a matrix of those differences.

toeplitz(0:numel(r)-1) gives you a matrix of integers with zeros along the main diagonal, ones along the 1 and -1 diagonals (one away from main), twos along the +/-2 diagonals, etc. This matrix is essentially the lengths of the base of each triangle divided by dt, i.e., it's j-i when you want to calculate |r(j)-r(i)|*dt*(j-i)/2 .

Then triu is used to keep only the upper triangular section of the matrix so that j>=i, i.e., only considering differences that look forward not backward.

Once the matrix of areas A is calculated, you can search through it like in my method before, except now you must jump rows basically. Starting with the first row, find the index where the area exceeds the threshold. Then move to that row to find the next one, etc.

r0 = 3; % start r
r1 = 6; % end r
t0 = 5; % halfway point
k = 2; % exponential growth
r_ = @(t) 1./(1+exp(-k*(t-t0))) * (r1-r0) + r0;
t = linspace(0,15,301);
r = r_(t);

dt = mean(diff(t));
threshold = 5e-3;

A = triu(dt*abs(r-r.')/2.*toeplitz(0:numel(r)-1));

lgc = [true false(1,numel(r)-1)];
row = 1;
while true
    idx = find(A(row,:) >= threshold, 1);
    if isempty(idx)
        break
    end
    lgc(idx+1) = true;
    row = idx;
end

hold on
plot(t,r)
plot(t(lgc),r(lgc),'x')

[Jack]

Thank you, this is exactly what I was looking for!

[Voss]

You're welcome!

[Jack]

I have one more question actually.

Is this adaptable to changing from the area of the triangle to the area under the curve (within the triangle)? I.e. the hypotenuse of the triangle from before becomes the curve between the two points, whereas the opposite and adjacent sides of the triangle stay the same and its now that area you calculate.

Perhaps some way of using cumtrapz, trapz or integral...

[Voss]

Maybe something like this. See if it seems right.

r0 = 3; % start r
r1 = 6; % end r
t0 = 5; % halfway point
k = 2; % exponential growth
r_ = @(t) 1./(1+exp(-k*(t-t0))) * (r1-r0) + r0;
t = linspace(0,15,301);
r = r_(t);

dt = mean(diff(t));

cA = cumtrapz(t,r-r(1));
dA = cA-cA.';

threshold = 5e-3;

lgc = [true false(1,numel(r)-1)];
row = 1;
while true
    idx = find(dA(row,:) >= threshold, 1);
    if isempty(idx)
        break
    end
    dA = dA-dA(row,idx);
    lgc(idx+1) = true;
    row = idx;
end

hold on
plot(t,r)
plot(t(lgc),r(lgc),'x')