Need help with least squares fit (lsqonlin)

This should be an easy one for someone I am a matlab newbie (have been for about a year now!) working with pre-existing code. The code was a genetic algorithm that found a bunch of parameters for me, and now I need to see how closely it fits data that I already have without re-running the whole thing, with the ability to manually change some of the parameters(for a sensitivity analysis)

The equations are in a file called WTTNF.m

function dy = WTTNF(~,y,par)
dy = zeros(4,1);
dy(1) = par.alpha*y(1)*(1-y(1)/y(4))-par.beta*y(1)-par.gamma*y(3)*y(1);  
dy(2) = par.beta*y(1)+par.zeta*par.gamma*y(3)*y(1)-par.delta*y(2); 
dy(3) = par.theta*y(2)+par.epsilon;                               
dy(4) = par.eta_c*y(3)*y(1)+par.eta_h*y(3)-par.psi*y(4)*(y(4)/par.Ko)^2;

The parameters are stored in par , as par.alpha, par.beta, etc.

The data is stored in a data , with time values and four sets of data points as follows:

>> data.t

ans =

     0     1     2     3     6     8
>>data.val

ans =

         132         183         358         448         563         526
         246         312         741         665         996        1126
         193         258         267         295         478         676
         185         207         489         574         758         292

Now, this function (WTTNF) works fine within the context of a much larger and more complicated (GA) genetic algorithm. At the end of the GA, it spits out a fitness value, but now I want to start manually messing with some of the parameters (for a sensitivity analysis) and see what fitness value I get, without running the whole GA again. Should be simple

My best guess was to do this,

>> lsqonlin(WTTNF,data)

??? Input argument "par" is undefined.

Error in ==> WTTNF at 7 dy(1) = par.alpha*y(1)*(1-y(1)/y(4))-par.beta*y(1)-par.gamma*y(3)*y(1);

At the time I run this, I do have par in memory with all its associated values for each parameter.

To see what I mean, I can type:

>> par

par =

      alpha: 0.1000
       beta: 0.0437
      gamma: 0.0273
       zeta: 0.4935
      delta: 6.2036
      eta_h: 6.6211
      eta_c: 0.0150
      theta: 0.0235
    epsilon: 5.0568
        psi: 0.0243
         Ko: 400

>> par.alpha

ans =

    0.1000

>> par.beta

ans =

    0.0437

This is probably something simple. Thanks for any help.

Regards

-Dave K

Answer 1

[Matt J]

but now I want to start manually messing with some of the parameters (for a sensitivity analysis) and see what fitness value I get, without running the whole GA again.

It can't be deduced from your post what the fitness function optimized by GA actually was, though a good guess is that it was norm(dy) or norm(dy)^2 where dy is the output of WTTNF. In any case, there is no clear reason why you would be using lsqnonlin if all you want to do is re-evaluate the fitness function at a given set of tweaked parameters. Just feed the new parameters to WTTNF again and calculate the fitness based on its output.

[David]

Here's how it's run in the GA:

[x, fvalN] = lsqnonlin(@F2,lb+rand(1,length(lb)).*(ub-lb),lb, ub,options, par,paramsToFit, data);

So the original fitness function is lsqonlin. If I'm not saying this correctly, I apologize.

When I try to run WTTNF on it's own I get:

>> WTTNF ??? Input argument "par" is undefined.

Error in ==> WTTNF at 7 dy(1) = par.alpha*y(1)*(1-y(1)/y(4))-par.beta*y(1)-par.gamma*y(3)*y(1);

And yet this function runs fine as part of the larger program.

Does that help? Pardon my ignorance. My math ability exceeds my programming skills by a longshot!

-Dave K

[Matt J]

It's still not at all clear what is being minimized. Since lsqnonlin, like GA, is also a function minimization routine, you are using GA to "minimize the minimum" of some function F2 that we still cannot see. If the whole idea is to solve the equations in WTTNF, I'm not sure why you wouldn't just apply fsolve() or lsqnonlin() to WTTNF directly.

When I try to run WTTNF on it's own I get: WTTNF ??? Input argument "par" is undefined.

Well, you can't just type WTTNF and be done with it. From the code you posted, WTTNF is expecting 3 input arguments. The error message is complaining that you did not give it all 3.

[David]

Ok, well obviously I'm using the completely wrong thing, or asking incorrectly, and I apologize for not clarifying how I was running WTTNF. But that's beside the point now. Let me just ask this way...

What I'm trying to do is see how well the equations in WTTNF match the given data when I start messing with the parameters. I thought lsqonlin would do that, since that's what is used in the GA. So the new question (and maybe I should just post a new question) is what should I be using to accomplish this?

Thanks for your patience.

Regards,

Dave K

[Matt J]

Again, it sounds like you should simply be calculating norm(dy) where dy is the output of WTTNF, called with whatever set of input data y and parameter data par that you are interested in. From the code you posted, dy appears to be the difference between the left and right hand side of your equations and therefore a very natural measure of how well they are satisfied.

[David]

Thanks Matt. I'll give it a try and respond again shortly. Sorry I'm not more fluent or I probably would have listened to you the first time!

[David]

I'm going to go ahead and mark this as answered and start a new question, as I am not sure how to use norm. Thanks for your help.