Matlab plot for exponential decay function

Question

I have empirical data of 9 sets of patients the data looks in this format

               input = [10      -1      1
                        20      17956   1
                        30      61096   1
                        40      31098   1
                        50      18446   1
                        60      12969   1
                        95      7932    1
                        120     6213    1
                        188     4414    1
                        240     3310    1
                        300     3329    1
                        610     2623    1
                        1200    1953    1
                        1800    1617    1
                        2490    1559    1
                        3000    1561    1
                        3635    1574    1
                        4205    1438    1
                        4788    1448    1
                      ];
                      calibrationfactor_wellcounter =1.841201569;

Here, the first column describes values of time and next one is concentration. As you can see, the concentration increases until a certain time and then decreases exponentially with increase in time.

If I plot the following characteristics, I obtain following curve enter image description here

I would like to create a script which represents the same behavior cited above. following is the script which i have formulated where concentration linearly increases till certain time period and aftermath it decays exponentially, but when i plot this function i am obtaining linear characteristics , kindly let me know if my logic is appropriate

function c_o = Sample_function(td,t_max,a1,a2,a3,b1,b2,b3)
   t =(0: 100 :5000); % time of the sample post injection in mins
   c =(0 : 2275.3 :113765);
   A_max= max(c);%Max value of Concentration (Peak of the curve)

   c_o = zeros(size(t));
   c_o(t>td & t<=t_max) = A_max*(t(t>td & t<=t_max)-td);

   c_o(t>t_max)=(a1*exp(-b1*(t(t>t_max)-t_max)))+(a2*exp(-b2*(t(t>t_max)-t_max)))+(a3*exp(-b3*(t(t>t_max)-t_max)));        

   fprintf('plotting Data ...\n');
   hold on;
   %figure ;
   plot(c_o,'erasemode','background');
   xlabel('time of the sample in minutes ');
   ylabel('Activity of the sample Ba/ml');
   title (' Input function: Activity sample VS time ');

   pause;
end

The figure i obtained is enter image description here In the above plot the decay is linear instead of exponential, let me know how to obtain 3rd order decay this is the line of code i have written to obtain 3rd order decay

c_o(t>t_max)=(a1*exp(-b1*(t(t>t_max)-t_max)))+(a2*exp(-b2*(t(t>t_max)-t_max)))+(a3*exp(-b3*(t(t>t_max)-t_max)));

possible duplicate of [Matlab :fitting 3rd order Exponential decay function](http://stackoverflow.com/questions/21795749/matlab-fitting-3rd-order-exponential-decay-function) — Wok, Feb 15 '14 at 16:37
Thank you, much better! But more questions: The green line in the first plot represents your empirical data, right? What is the blue line in that plot? – What values did you use for the parameters in your `Sample_function`? – How did you fit your function to your data? — A. Donda, Feb 15 '14 at 16:37
@Wok, the poster posted this question several times, but I think this is the most detailed and useful version of the question. I asked the moderator to mark the other ones as duplicates. — A. Donda, Feb 15 '14 at 16:38
@Devak You want to fit a sum of 3 exponentials (six parameters to learn) to less than 20 data points. Good luck finding relevant parameters, I guess we could find an infinity of them. — Wok, Feb 15 '14 at 16:40
Yes the Green line represents my empirical data,and the values are just mere sample values of a patient, the values might change if i take another clinical report, I have create a model function which accurately represents the behavior cited in Figure 1, The behavior will have some delay time, aftermath a steep increase in concentration till certain time, then an exponential decrease..The behavior describes 2-Compartment tissue model input function — Devak, Feb 15 '14 at 16:41
Okay, so you are not interested in fitting the model. Fine. Well, if your problem is the plot looks piece-wise linear, then increase the number of time points. Or revise your estimation procedure of the 2-compartment parameters. — Wok, Feb 15 '14 at 16:42
Blue line is reference i dont need to plot it , i just need to formulate the behavior described in green line. the Values are Td.. delay time, T_max,peak time,a1,a2,a3 are coefficients of exponential decay so are the values b1,b2,b3 — Devak, Feb 15 '14 at 16:44
@Devak, sure but what are the actual values of these parameters that you used for your second plot, and how did you get them? — A. Donda, Feb 15 '14 at 16:45
@wok, model may not represent the same characteristics, the whole point is behavior must be replicated, the behavior must increase linearly till certain time,I have called it as t_max and aftermath it must decrease exponentially its 3rd order decay as i need to use reference values from patients, 3rd order decay represents the characteristics more accurately so i used 2 tissue compartment model — Devak, Feb 15 '14 at 16:47
@Doonda, i just added them, if you observe the input empirical data, till first 10 minutes there is a delay in the function, delay sometimes varies may be sometimes 15 minutes, so Td value varies from 10-15, the concentration usually peaks by 60 minutes, sometimes 40 minutes so t_max range is 40-60 and i use some arbitrary values for a1,a2,a3 and b1,b2,b3 to obtain exponential decay charecteristics — Devak, Feb 15 '14 at 16:50
@Doonda My objective is to write a function which models the behavior, the actual values of Td-delay time,T_max,maximum time and coefficients a1,a2,a3 and b1,b2,b3 are user given ... if you read the code I have described the vector representing the model of concentration values as c =(0 : 2275.3 :113765); — Devak, Feb 15 '14 at 16:53
@Doonda anymore queries , please let me know, also kindly review if my logic of formulating the model is accurate — Devak, Feb 15 '14 at 17:08
Devak, I think I understand now, but a quick approach on my part didn't work out so far. I'll get back to you later. — A. Donda, Feb 15 '14 at 17:13
@Doonda the way I formulated the equation which represents the model is as follows : 1) If t>=0 & t<=td ,Concentration will be 0 2)if t>td & t<=t_max ,concentration will be A_max(t-td), where A_max is maximum value of concentration and 3)if t>t_max, the concentration value can be expressed as Summation of(a_n*exp(-b_n(t-t_max)) where n varies from 1 to 3. this is the representation of 3rd order exponential decay, now i have written a sample script for the same model, you can review the code which i have posted and plot it,u will obtain plot 2 — Devak, Feb 15 '14 at 17:22
@A.Donda Kindly post the script if you have worked on this question, any help is greatly appreciated — Devak, Feb 15 '14 at 19:39
@wok :Could you be more descriptive , how can i formulate this using my script, what do i have to adopt into my function — Devak, Feb 15 '14 at 19:41

A. Donda · Accepted Answer · 2014-02-16T17:05:40.423

I've come up with a solution using the functionality of Matlab's Curve Fitting Toolbox. The fitting result looks very good. However, I've found that it strongly depends on the right choice of starting values for the parameters, which therefore have to be carefully chosen manually.

Starting from you variable input, let's define the independent and dependent variables for the fit, time and concentration,

t = input(:, 1);
c = input(:, 2);

and plot them:

plot(t, c, 'x')
axis([-100 5000 -2000 80000])
xlabel time
ylabel concentration

data

These data are to be modeled with a function with three pieces: 1) constantly 0 up to a time td, 2) linearly increasing between td and tmax, 3) decreasing as a sum of three different exponentials after time tmax. In addition, the function is continuous, so that the three pieces have to fit together seamlessly. The implementation of this model as a Matlab function:

function c = model(t, a1, a2, a3, b1, b2, b3, td, tmax)

c = zeros(size(t));

ind = (t > td) & (t < tmax);
c(ind) = (t(ind) - td) ./ (tmax - td) * (a1 + a2 + a3);

ind = (t >= tmax);
c(ind) = a1 * exp(-b1 * (t(ind) - tmax)) ...
    + a2 * exp(-b2 * (t(ind) - tmax)) + a3 * exp(-b3 * (t(ind) - tmax));

Model parameters appear to be treated internally by the Curve Fitting Toolbox as a vector ordered alphabetically by the parameter names, so to avoid confusion I sorted the parameters alphabetically in the definition of this function, too. a1 to a3 and b1 to b3 are the amplitudes and inverse time constants of the three exponentials, respectively.

Let's fit the model to the data:

ft = fittype('model(t, a1, a2, a3, b1, b2, b3, td, tmax)', 'independent', 't');
fo = fit(t, c, ft, ...
    'StartPoint', [20000, 20000, 20000, 0.01, 0.01, 0.01, 10, 30], ...
    'Lower', [0, 0, 0, 0, 0, 0, 0, 0])

As mentioned before, the fitting works well only if the algorithm gets decent starting values. I here chose for the amplitudes a1 to a3 the number 20000, which is about one third of the maximum of the data, for b1 to b3 a value of 0.01 corresponding to a time constant of about 100, the time of the data maximum, 30, for tmax, and 10 as a rough estimate of the starting constant time td.

The output of fit:

fo = 

     General model:
     fo(t) = model(t, a1, a2, a3, b1, b2, b3, td, tmax)
     Coefficients (with 95% confidence bounds):
       a1 =        2510  (-2.48e+07, 2.481e+07)
       a2 =   1.044e+04  (-7.393e+09, 7.393e+09)
       a3 =   6.506e+04  (-4.01e+11, 4.01e+11)
       b1 =   0.0001465  (7.005e-05, 0.0002229)
       b2 =     0.01049  (0.006933, 0.01405)
       b3 =     0.09134  (0.08623, 0.09644)
       td =       17.97  (-3.396e+07, 3.396e+07)
       tmax =       26.78  (-6.748e+07, 6.748e+07)

I can't decide whether these values make sense physiologically. The estimates also don't appear to be too well defined, since many of the confidence intervals are huge and actually include 0. The documentation isn't clear on this, but I assume the confidence bounds are nonsimultaneous, which means it is possible that the huge intervals simply indicate strong correlations between the estimates of different parameters.

Plotting the data together with the fitted model

plot(t, c, 'x')
hold all
ts = 0 : 5000;
plot(ts, model(ts, fo.a1, fo.a2, fo.a3, fo.b1, fo.b2, fo.b3, fo.td, fo.tmax))
axis([-100 5000 -2000 80000])
xlabel time
ylabel concentration

shows that the fit is excellent:

fit

A close-up of the more interesting initial part:

close-up

Note that the estimated time and value of the true maximal concentration (27, 78000) depends only on the fit to the following decreasing part of the data, since the linear increase is characterized only by one data point, which does not pose a constraint.

The results indicate that the data are not sufficient to obtain precise estimates of the model parameters. You should consider either to increase the sampling rate of the data, particularly up to time 500, or to decrease the complexity of the model, e.g. by using a sum of two exponentials only; or both.

Donda ,Wow, thanks i will work out what you have stated here, and see if i can replicate them, please be in touch, thanks a lot — Devak, Feb 16 '14 at 17:37
I have one more query, i would like to perform Circular convolution of the above concentration equations with 2-tissue compartment model output equation.. for which i have to take the equation of Concentration from t>td & t<=tmax and multiply with 2-tissue compartment model and plot it, same for the Concentration values t> t_max ,i have to perform convolution with 2-Tissue compartmental models output equation and plot it..Do you have any ideas how i can approach it, i would like to write another function which takes the concentration values at t>td &t<=t_max and t>t_max, how to approach here — Devak, Feb 16 '14 at 20:12
Good answer. In my opinion, one exponential would allow a more realistic parameter estimation, which would be just as good in terms of error fit. — Wok, Feb 16 '14 at 21:19
@Wok Could you please let me know how i can possibly approach this problem which i have stated here http://stackoverflow.com/questions/21819133/matlab-continuous-convolution-and-plotting — Devak, Feb 17 '14 at 01:36
@Devak, why did you unaccept? My answer does answer your question, doesn't it? — A. Donda, Feb 17 '14 at 11:39
@A.Donda, i accepted your answer and upvoted it, more over i am interested in another query if you please review my question and let me know how i can approach it, i have stated the question here in this link http://stackoverflow.com/questions/21819133/matlab-continuous-convolution-and-plotting — Devak, Feb 17 '14 at 11:50
@Devak, it is not accepted any more. You can't accept more than one answer on a question, and I guess you accepted Wok's answer, therey unaccepting mine. You have to decide which one is the best answer, and no offence to Wok, but I think it is mine... — A. Donda, Feb 17 '14 at 11:53
@Wok, I tried it with one exponential, and it wasn't a good fit; the decay is too sharp at the beginning and too slow later. Two exponentials might work though. — A. Donda, Feb 17 '14 at 11:54
@A.Donda,I will un accept woks answer and will accept yours mean while could you please let me know how i can resolve the other query , i have provided the link for the same — Devak, Feb 17 '14 at 11:56
@Devak, I will have a look at you other question when I have the time, and maybe I have an idea. But, you need to understand that SO is not a forum to help people out with whatever they need at the moment as fast as possible – it is a place to ask good questions and get good answers. I think you are not patient enough. You can't repost the same question several times if there are no answers yet – people here help if they find the time in their lifes, which are not *mainly* about helping others. — A. Donda, Feb 17 '14 at 12:04
Moreover, when you post a question, take the time to really think it through, expose the problem clearly and without clutter, include all necessary information, etc. You expect people to be thorough and careful on their answers – they expect the same from your questions. Please take a moment to look at http://stackoverflow.com/tour, and at other questions and answers, to get an understanding how SO works. You can get a lot out of this site, it is extremely helpful – but only if you're willing to make an effort yourself! — A. Donda, Feb 17 '14 at 12:05
@A.Donda, I am just interested in plotting the convolution integral of the same function which was obtained by your earlier with another vector, i have performed cconv on 2 equations at respective time intervals and plotted them using plot(t,c_t),where c_t is convolution output function,but i am incurring assignment error about the number of elements — Devak, Feb 17 '14 at 12:21
@Devak The accepted answer should be the most useful answer. Was your problem a plotting problem or a parameter estimation problem? — Wok, Feb 17 '14 at 15:35

score 0 · Answer 2 · edited May 23 '17 at 12:07

0

Try this code from this question:

   x = input(:,1);
   c = input(:,2);
   c_0 = piecewiseFunction(x, max(c), td,t_max,a1,a2,a3,b1,b2,b3)

with:

function y = piecewiseFunction(x,A_max,td,t_max,a1,a2,a3,b1,b2,b3)
y = zeros(size(x));

for i = 1:length(x)
    if x(i) < td
        y(i) = 0;
    elseif(x(i) < t_max)
        y(i) = A_max*(x(i)-td);
    else
        y(i) = (a1*exp(-b1*(x(i)-t_max)))+(a2*exp(-b2*(x(i)- t_max)))+(a3*exp(-b3*(x(i)-t_max)))
    end
end

end

edited May 23 '17 at 12:07

Community

1
1

answered Feb 15 '14 at 19:46

Wok

4,956
7
42
64

T_d value changes accordingly with respective tracer, some tracers have higher delay time and how have lower, so value varies – Devak Feb 15 '14 at 19:54
All the values td,t_max,a1,a2,a3 and b1,b2.b3 change from patient to patient, so i would like to formulate the plot for different values of td,t_max,a1,a2 and b1... respectively – Devak Feb 15 '14 at 20:03
td :30 minutes ; t_max :100;a1=1231,a2=1442,a3=1125,b1=12,b2=14,b3=16...These are sample values of 1 of the patient – Devak Feb 15 '14 at 20:20
The code which you have mentioned for piece wise function is throwing errors could you be more descriptive do i have to describe the values x,c,and c_o outside the function? – Devak Feb 16 '14 at 13:42
Yes, you save the function as a .m file and you call it from Matlab. – Wok Feb 16 '14 at 15:49
MATLAB absolutely supports logical expressions. In fact, if you use `find`, the code checker in the editor will throw a warning telling you to use logical indexing instead because it is faster. – craigim Feb 16 '14 at 23:54
@craigim Correct. I have removed this part of my answer. – Wok Feb 17 '14 at 10:49

Matlab plot for exponential decay function

2 Answers2