How do I define a function to be used in an annonymous function?

Hello, I am trying to use the following approach to find the values p(1) and p(2) that gives the best aproximation of curves yreal and yfit:

%Function to calculate the sum of residuals for a given p1 and p2
fun = @(p) sum((yreal - yfit(p(1),p(2),signal1, signal2)).^2);

%starting guess
pguess = [1,0.2];

%optimise
[p,fminres] = fminsearch(fun,pguess)

where y_fit is a function that depends on those parameters and signals 1 and 2. I am building the function yfit because inside it I have to use ode45 to obtain the y_real values, I mean, it is not a function as easy as, for example:

fun = @(p) sum((yreal - p(1)*cos(p(2))*signal1+signal2)).^2);

The problem I am having with my function yfit is:

Undefined function 'minus' for input arguments of type 'struct'.

Error in fun = @(p) sum((yreal - yfit(p(1),p(2),signal1, signal2)).^2);

Any suggestion on how should I define fun or yfit?

Thanks

[Star Strider]

Your OLS cost function fun looks correct.

Post your code for yfit, since that seems to be where the problems are.

EDIT — Also please post your original equations, including your ODE and your call to ode45. We’re good here, but we are hopeless mind-readers, and can’t figure out what you want to do unless you tell us.

[dpb]

What's the definition of yfit look like and what are signal1/2 ? I'm puzzled about what generated the reference to a structure.

W/O any more info, I'd start at

fun = @(p,s1,s2) sum((yreal - yfit(p(1),p(2),s1, s2)).^2);

but not positive about it.

One thing I see as a problem in the initial is that for any variable in the anonymous function not in the argument list, the current value of that variable at the time of the creation of the handle is "builtin" to the function and changing its value (or even deleting it) won't change the values used in the function. Hence, in your initial case signal1 and -2 will never change during the search if they're supposed to.

[Fran]

Lets see. Vars "time", "yreal" and "e" are read from a file, but here I've written somethig to start with. This is the main program:

%---- MAIN
time = 0:0.1:10;                    % Example here
yreal = rand(length(t),1)'+sin(t);  % Example here
e = 0.1.*sin(t)                     % Example here

%Function to calculate the sum of residuals for a given p1 and p2
fun = @(p) sum((yreal - yfit(time, e, p(1),p(2))).^2);

%starting guess
pguess = [1,0.2];

%optimise
[p,fminres] = fminsearch(fun,pguess)
%----

and here you have the function yfit, which is saved in yfit.m and contains the two following functions:

%---- yfit
function act = yfit(time, e, PCSA, Cfast)

      dm = 1054;                
      lF = 0.0995828022467124;  
      m = PCSA * lF * dm;

      global Te Trise Tfall 
      Te = (25+0.1*m*(1-Cfast))/1000;
      Trise = (5+0.2*m*(1-Cfast)^2)/1000*15;
      Tfall = (30+0.5*m*(1-Cfast)^2)/1000*19;

      x0 = [0.00001, 0];
      tspan = time;

      a = ode45(@FES, [0;1],[0 0], [], e, time);
      act = a(:,1);

end

function dx = FES(t, x, e, t_ex)

      global Te Trise Tfall

      e_int= interp1(t_ex,e,t);
      e_int2=interp1(t_ex,e,t+0.01);

      if t<(t_ex(end))
          if e_int < e_int2
              k1 = Te*Trise;
              k2 = Trise+Te;
          else
              k1 = Te*Tfall;
              k2 = Tfall+Te;
          end
      else
          k1 = Te*Tfall;
          k2 = Tfall+Te;
      end

      dx = zeros(2,1);
      dx(1) = x(2);
      dx(2) = 1/k1*(-k2*x(2)-x(1)+e_int);
  end
  % ----

Thanks for your help

Answer 1

[Star Strider]

First, since yfit is an external function (we didn’t know that before), change fun to:

fun = @(p) sum((yreal - @yfit(time, e, p(1),p(2))).^2);

and see if that works. (Note the added ‘@’ sign.)

If you still have problems with it, the next step is to call yfit itself (outside of fun) and see what the results are. If yfit has problems in that situation, start troubleshooting by being sure your ODE in FES integrates.

[Fran]

Thanks for your proposals, but I haven't been able to achieve the results. I've been all day going forward and back withouth success. I have modified the script and now I have the following:

Main program:

%---- MAIN
time = 0:0.1:10;                    % Example here
yreal = rand(length(t),1)'+sin(t);  % Example here
e = 0.1.*sin(t)                     % Example here

save fitData.mat  time e
save realData.mat time yreal

%starting guess
pguess = [0.0062 0.52];

%optimise
options = optimset('Display', 'iter');
[p,fminres] = fminsearch(@sumatorio, pguess, options);
%----

Function sumatorio (this one is new), saved in sumatorio.m:

%----
function J = sumatorio(p)
    load realData.mat
    PCSA = p(1);
    Cfast = p(2);
    J = sum((act_TA - ex2act(time, PCSA, Cfast)).^2); 
end
%----

Function ex2act (in the previous comment was function yfit). It is practically the same but here I load some values. It is saved on ex2act.m

%----
function act = ex2act(timeSignal, PCSA, Cfast)
      x0 = [0.00001, 0];
      tspan=timeSignal;
      % Gfolher et al. J. of Mech Med Biol, 2004, 4, 77-92
      dm = 1054;             
      lF = 0.0995828022467124;
      m = PCSA * lF * dm;
      global Te Trise Tfall
      Te = (25+0.1*m*(1-Cfast))/1000;
      Trise = (5+0.2*m*(1-Cfast)^2)/1000*15;
      Tfall = (30+0.5*m*(1-Cfast)^2)/1000*19;
      [time, a] = ode45(@FES, tspan, x0);
      act = a(:,1);
  end

% FES
function dx = FES(t,x)
    load fitData.m
    global Te Trise Tfall       
    e_int= interp1(time,e,t);
    e_int2=interp1(time,e,t+0.01);
      if t<(time(end))
          if e_int < e_int2
              k1 = Te*Trise;
              k2 = Trise+Te;
          else
              k1 = Te*Tfall;
              k2 = Tfall+Te;
          end
      else
          k1 = Te*Tfall;
          k2 = Tfall+Te;
      end
      dx = zeros(2,1);
      dx(1) = x(2);
      dx(2) = 1/k1*(-k2*x(2)-x(1)+e_int);
  end
  %----

With this mess it seems to work, at least it enters on fminsearch, but it doesn't goes out. There is a a big problem of inneficiency. Right now I am completely lost. I wait your comments. Thanks

[Star Strider]

At least we’re making some progress.

What do you mean by ‘it doesn't goes out’? What is not occurring that should be? Are you getting any errors?

If it will not produce parameters that fit to your data, there could be any number of reasons, the most likely being a choice of initial parameter estimates that are quite far from the correct ones.

If you want to upload the PDF the paper you referred to (Gfolher et al. J. of Mech Med Biol, 2004, 4, 77-92) as an attachment to your original question (use the ‘paper-clip’ icon), it might help me understand what you are doing. (No promises!)

[Fran]

I think it's easier and faster if I explain it to you.

I have two signals, one from experiments and another one calculated. This last one depends on two parameters and I want to know which are those parameters to obtain the best fit. That is, I am trying to minimize the sum of squared differences varying those parameters:

min sum(act_TA - act(p1,p2))^2

act_TA comes from experiments and act is the signal that I get after integration of a second order differential equation:

k1 * act_dd + k2 * act_d + act = e

where _dd indicates second derivativa, _d first derivative, k1 and k2 are calculated according the previous comments (Gohfler et. al paper) and e is a prescribed signal. Theory is easy, efficient practice is complicated.

When I said it doesn't goes out I meant it takes too too long. Right now I've left fminsearch and I am doing the following:

p1 = 0.0045:0.0005:0.0085;
p2 = 0.1:0.1:0.9;

for ind1 = 1:length(p1)
    P(1) = p1(ind1);
    for ind2 = 1:length(p2)
        P(2) = p2(ind2);
        J(ind1,ind2) = sumatorio(P);
        toc
    end
end

It has been running during an hour on a pentium i3 core and right now it is on:

ind1 = 4; ind2 = 4;

[Star Strider]

First question: If you give your DE the appropriate parameter estimates and initial conditions, how long does it take it to integrate?

Second question: I don’t understand the ind1 ... ind2 loop. It seems to me that you just have to do the fminsearch call once comparing your two signals, be sure you have a good fit to your data and parameter estimates that make sense, and you’ve done what you set out to do.

I admit to not understanding what you are doing, most significantly with respect to your loop. I thought you are doing a relatively straightforward curve-fitting parameter estimation, but apparently I missed something.

I suggest you output and save the p and fminres variables at each iteration. If it’s taking as long as it is, you don’t want to risk losing any information.