parameter estimation and minimization

I have to estimate parameters(k1 ,k2) of the following equation

dx/dt=k1(3-x)-k2*ca0*x^2

I have data available of x vs t at different inital ca0.I have minimization function as follows.

where N=total no. of runs at different ca0(=6) and N'= total no. of runs at same ca0(=12). kindly help me to build a matlab code for the same bcoz i have done only with one summation error minimization.

[Devyani]

my code looks like this

xt_data=[15 0.1
       12 0.3];
Ca0=[6 9 12 18];
sum=0;
for j=1:6
  [x t]=ode45(@paraesti,[0 0]);
  for i=1:12
      sum=sum+(xt_data(j,i)-x(j,i))^2;
  end

[Star Strider]

You only have 4 data points and 4 values of Ca0 and you want to fit 2 parameters?

I don’t understand the 6 and 12 iteration loops.

[Devyani]

i am sorry actually i was just trying with four random data points.. i posted tht code only,change the loop to 4 thn.and the function paraesti looks like this

function dxdt=paraesti(x,t,theta)
global theta;
k1=theta(1);
k2=theta(2);
dxdt=k1*(3-x)-k2*ca0*x^2;

[Devyani]

and i have actually 6 different values of ca0 and each ca0 has 12 data points.thts y the loop 6 and 12

Answer 1

[Star Strider]

Your ODE is not well-behaved. I had problems fitting it with synthesised data with even a small noise component, particularly with respect to k1 (or p(1) in my code). It is very sensitive to the initial parameter estimates. I have no idea what your parameters actually are, so you will have to experiment with the starting estimates to get a good fit to your data. I also used ‘0’ as an initial condition for x in the data synthesis and the DevyaniODEFIT function. Change this if it is not correct. I used a time vector of [0:12] for the independent variable. Change this to reflect the values of your actual independent variable.

I attached the objective function that integrates your ODE and is used by the main script as an objective function that the parameter estimation routine uses to fit your data. I used fminsearch here because I do not know if you have access to the Statistics Toolbox function nlinfit or the Optimization Toolbox function lsqcurvefit. They are more robust and will give a better fit than fminsearch, so you can easily change my code here to use them instead. Note that your paraesti function does not pass Ca0 as a parameter, so either include it in yours, or use my function instead.

The main code:

Ca0=[6 9 12 18];

DE = @(p,t,x,Ca0) p(1).*(3-x)-p(2).*Ca0.*x.^2;      % ODE function to create data

% CREATE DATA:
for k1 = 1:length(Ca0)
    p(:,k1) = rand(2, 1)*0.1;                       % Random parameters k1, k2
    [t x(:,k1)] = ode45(@(t,x) DE(p(:,k1),t,x,Ca0(k1)), [0:12], 0);
    x(:,k1) = x(:,k1) + 0.005*randn(size(x(:,k1))); % Add noise for realism
end

% ESTIMATE PARAMETERS: 
for k1 = 1:length(Ca0)
    xest = @(B) DevyaniODEFIT(B,t,Ca0(k1));         % Calls ODE function
    OLS = @(B) sum( (x(:,1) - xest(B) ).^2);        % Least squares objective function
    B0 = rand(2,1)*0.1;                             % Initial parameter estimates
    B(:,k1) = fminsearch(OLS, B0);                  % Optimise parameters to fit data
    xgrf(:,k1) = DevyaniODEFIT(B,t,Ca0(k1));        % Generated fitted data to plot
end

figure(1)
plot(t, x, '-x', 'MarkerSize',2, 'LineWidth',1)     % Plot original data
hold on                                             % Plot estimated data
plot(t, xgrf, '-p', 'LineWidth',1.5, 'MarkerSize',3);
hold off
grid

The function DevyaniODEFIT.m is attached.

The code definitely works. I documented it with comments as well as I can, so you can easily understand it. You will probably have to experiment with different starting values it to fit your data.

[Devyani]

thanku so much fr your answer, but I am trying to understand the code, but i have got the hang of it.Thank you for your time!!! :) :)

[Star Strider]

My pleasure!

[Devyani]

I am afraid but when i used the code, i am getting following error

Assignment has more non-singleton rhs dimensions than non-singleton
subscripts

Error in ==> Untitled3 at 6
    [t x(:,k1)] = ode45(@(t,x) DE(p(:,k1),t,x,Ca0(k1)), [0:12], 0);

[Star Strider]

I don’t know why you’re getting that error. It ran for me without problems or I’d not have posted it. (I have R2014a, but I don’t believe version differences could cause that error.)

I use that entire loop to create test data to use for the regression (since I didn’t have your data), and to demonstrate the way you need to orient your data for the parameter estimation. You can delete that loop, and substitute your own values of x and t in the regression.

[Devyani]

I dont knw ,although did my minimization using different program .I first solved the differential equation symbolically and then minimized the function by substituting the values of Ca0 and t :p :)

[Devyani]

I copied the same code and didnt change anything not even data

??? Error using ==> minus
Matrix dimensions must agree.

Error in ==> @(B)sum((x(:,1)-xest(B)).^2)

Error in ==> fminsearch at 320
    x(:) = xr; fxr = funfcn(x,varargin{:});

Error in ==> matlabanswer at 14
    B(:,k1) = fminsearch(OLS, B0);

[Star Strider]

Without seeing the function you used (from your solved differential equation), I can’t tell you how to change the code to do the curve-fit. Post it and I’ll see what I can do.

[Devyani]

Also while estimating parameters you have given initial values B0 after estimating it from ode45 , how is that possible?

 for k1 = 1:length(Ca0)
    xest = @(B) DevyaniODEFIT(B,t,Ca0(k1));         % Calls ODE function
    OLS = @(B) sum( (x(:,1) - xest(B) ).^2);        % Least squares objective function
    B0 = rand(2,1)*0.1;                             % Initial parameter estimates
    B(:,k1) = fminsearch(OLS, B0);                  % Optimise parameters to fit data
    xgrf(:,k1) = DevyaniODEFIT(B,t,Ca0(k1));        % Generated fitted data to plot
end

How are you calculating xest if no initial parameter p(1) is provided?

[Star Strider]

Your function is very sensitive to the initial estimates for the parameters, and so may not always be able to fit the test data (or your actual data). That generates the error if the ODE solver cannot integrate the ODE with the parameter estimates provided.

The fminsearch function calls the OLS function that calls the xest function. The initial parameter estimates are the elements of the B0 vector supplied to fminsearch.

[Devyani]

Matlab goes step by step i suppose.So if the line xest is given before the initial estimates are given ,it will show error and stop there before it can go ahead to see the initial estimates. I dont know if i m right. I will appreciate if you can elaborate on this more.

And i again copied the same code as it is, it is showing error, i m using ur data only. I dont know how is it working for you ? I am using Matlab 2008

Also in objective function you are subtracting elements of first column of x everytime with xest.They should also change as Ca0 is changing.If they are changing then how?

[Star Strider]

There could be version differences. I’m running R2014a and it runs for me without problems. I have no idea what your data are, so I have no idea what initial estimates might be most appropriate. Instead of using the random number generator, use:

B0 = [1; 1]*0.5;

It has very poor convergence with even synthetic data, something not common in my experience, so parameter updates that fminsearch generates that make your ODE unstable could be part of the problem.

[Devyani]

ok thanks ... i also think after two runs it is not able to solve the differential equation., may be that is the problem. thankyou for ur time!

[Star Strider]

My pleasure!

Take a look at the plots to get an idea of the problem. The differential equation is very sensitive to its initial conditions (so you may have to experiment to get the correct ones), and the parameter estimation is very slow to converge on a solution that is not always close to the actual synthetic data. The regression surface must be very shallow (low gradient). (This same approach has worked well in other situations and with other differential equations, even with fminsearch, although I prefer the Statistics Toolbox nlinfit.) Your actual data may be easier to fit than my test data, but I don’t have your data to fit.

It occasionally crashes for me as well, but restarting it usually solves the problem. It is successful with good initial parameter estimates, but the parameter updates issued by fminsearch can make it unstable. This is simply the nature of the differential equation and your data.

Curve fitting is heuristic, so it may take several initial estimates to get a good fit and reliable converged parameter estimates. I also used an initial condition of zero. Change it if that is not appropriate. You can even make the initial condition one of the parameters to be estimated if you know it is something other than zero. Pass it as one of the parameters in the parameter vector, and use it as an initial condition in your differenttial equation within DevyaniODEFIT.m.

[Devyani]

thanks star!! i got to know the prblem and have estimated the parameters with my data :) again thanku !

[Star Strider]

My pleasure!

I’m glad it worked.