How do I dynamically generate a model function for fitting or minimization?

For fitting response data with lsqcurvefit() or fitnlm(), I need a model function such as this Gaussian

G = @(A,X) A(1).*exp(-(X-A(2)).^2/(2*A(3).^2));

I work with response data with several features that I would like to fit. If I know the number of features beforehand, I can build a more complex model function "by hand", for example this one:

G6 = @(A,X) ...
    G(A(1:3),X)+...
    G(A(4:6),X)+...
    G(A(7:9),X)+...
    G(A(10:12),X)+...
    G(A(13:15),X)+...
    G(A(16:18),X);

However, I would like to be able to generate model functions with an arbitrary number of terms programatically. Typically my terms will be Gaussians or Lorentzians. I imagine calling a function to generate a function handle

modelfun = sum_of_n_terms(T,n)

, where T would be the function handle of a term (e.g., G in the above example), and n the number of times that this term should appear in modelfun. I read through the documentation on function handles and nested functions, but can't figure out how to apply this information to my problem, where the number of terms changes (rather than just some parameters of the same term). Any suggestions would be greatly appreciated.

Answer 1

[harsha001]

You can create a separate function file as :

function out = sum_of_n_terms( T, n, A, X )
% send fn handle T, and number of terms n, as well as inputs A and X
    out = zeros(size(X));
    for jj=1:n
    % Code for adding terms
         out = out + T(A((jj-1)*3 + (1:3)),X);
    end

end

And then create your fn handle and call:

G = @(A,X) A(1).*exp(-(X-A(2)).^2/(2*A(3).^2));

n = 1;

sum_of_n_terms(G, n, A, X)

Can be done for any arbitrary fn with different number of terms as long as you can describe it using the term number. Example: an expansion series approximation for sin(x):

function out = sine_exp( theta, n )
% get sin(theta) using upto n terms of maclurin expansion where theta is a vector of angles in radians
    out = zeros(size(theta));
    max_n = 100; n = min(n, max_n);    %maximum 100 terms
    theta = mod(theta, 2*pi );
    for jj=1:n
        jjth_term= (-1)^(jj-1) / factorial(2*jj-1) .* theta.^(2*jj-1)  ;
        out = out+jjth_term;
    end
end

And from cmdline:

fn1 = @(x,n) sine_exp(x,n);

t = 0:0.1:100;

plot(t, fn1(t, 20)) %get sin(t) using 20-term expansion

[Derk Joester]

Thank you and my apologies for the late response. It looks trivial in retrospect, but the crux seems to be the statement

out = zeros(size(theta));