How to call a function with vector input in the fit type function?

I am trying to fit an equation to a model and I need to call a function "kkrebook2" in my fit type funtion (For "kkrebook2", the inputs are two vectors and its output is a vector). This is the snippet of my fit type function code where I call "kkrebook2" function. But it doesn't work when I call my fit type function "SCR" in my fit code. Could you please let me know what is the issue of my code?

I think it is related to the way I am calling "kkrebook2" in my fit type function "SCR" but I don't figure out how to fix it. I have attached the "kkrebook2" function as well.

Thank you in advance!

function p = SCR(nu,numGaussians,a,center,sigma)
for j = 1:length(nu)
    for i = 1:length(nu)
        for k = 1 : numGaussians
            if k==1   
                b = Start;
            else 
                b = Start + (center * (k-1));
            end
            thisGaussian(i) = a.*exp(-((nu(i)-b).^2)/(2.*(sigma.^2)));
            % Add into accumulator array:
            gaussEqn1(i) = gaussEqn1(i) + thisGaussian(i);
            z(i)= gaussEqn1(i);   
        end
    end
    Rn(:) = kkrebook2(nu,z,0);   
    hi2(j) = (4*pi*n.*nu(j)*L)/c;
    hi1(j)=(Rn(j));
    f(j)=(-r1+r2*exp(-(z(j))).*exp(-1i*hi1(j)).*exp(-1i*hi2(j)));
    g(j)=(1-r1*r2*exp(-(z(j))).*exp(-1i*hi1(j)).*exp(-1i*hi2(j)));
    h(j)= f(j)./g(j);
    m(j)= abs(h(j));
    p(j)= (m(j).^2);       
end
end

[dpb]

What is the model you're trying to fit -- having other than a vector output of the functional is not valid in fit--it doesn't handle the multivariate case at all; it will and can only accept a vector as the reponse data to fit.

Unless you can arrange your function to output the results of the two vectors as one, you simply will not be able to use fit

It's possible depending upon just what your model is that you might be able to write the function in an acceptable format, but without understanding what it is you're trying to fit as the underlying model and what are the independent and dependent variables and coefficients to estimate, it's not possible to say how, or even if it is possible.

[Catalytic]

The parameter "Start" is never passed to your function. Also, it looks like the inner loop over i is independent of j. You therefore can and, for efficiency's sake, should hoist it out.

function p = SCR(nu,numGaussians,a,center,sigma)

    for i = 1:length(nu)
        for k = 1 : numGaussians
            if k==1   
                b = Start;
            else 
                b = Start + (center * (k-1));
            end
            thisGaussian(i) = a.*exp(-((nu(i)-b).^2)/(2.*(sigma.^2)));
            % Add into accumulator array:
            gaussEqn1(i) = gaussEqn1(i) + thisGaussian(i);
            z(i)= gaussEqn1(i);   
        end
    end

  Rn(:) = kkrebook2(nu,z,0);   
    
 for j = 1:length(nu)


    hi2(j) = (4*pi*n.*nu(j)*L)/c;
    hi1(j)=(Rn(j));
    f(j)=(-r1+r2*exp(-(z(j))).*exp(-1i*hi1(j)).*exp(-1i*hi2(j)));
    g(j)=(1-r1*r2*exp(-(z(j))).*exp(-1i*hi1(j)).*exp(-1i*hi2(j)));
    h(j)= f(j)./g(j);
    m(j)= abs(h(j));
    p(j)= (m(j).^2);       
 end
end

[Shaily_T]

Thanks @Catalytic. I have already tried writing my fit type function similar to the above code which raised an issue I have posted here https://www.mathworks.com/matlabcentral/answers/1729655-inconsistency-in-the-fitted-curve-and-the-plot-of-obtained-fitting-parameters?s_tid=mlc_ans_email_view#comment_2188115

[William Rose]

@Shaily_T,

When ever I write a non-trivial function in Matlab, I also write a relatively simple script to test it. I encourage you to do the same. I am attaching examples of what I mean. You can modify the test scripts to pass parameters that are more realistic for your situation. I do not know what typical values are for the arguments to these functions.

The test script for kkrebook2.m works fine, when omega nd imchi are vector of length 2 or vectors of length 10. The script kkrebook2test.m fails if omega and imchi are scalars. This is to be expected, because if omega is a scalar then deltaomega cannot be computed.

The test script for SCR fails, because SCR fails before it even tries to call kkrebook2(). SCR crashes before that, because it does not know what "Start" is. This was already correctly pointed out by @Catalytic. Fix that first.

Run the simple test scripts, and make adjustments to the functions SCR() and kkrebook2(), until you have the functions returning reasonable results. Then try calling SCR from your more complicated parent program.

Good luck.

[Shaily_T]

@William Rose

Thanks for your response and consideration!

Yes, that is correct I have "Start" in my code but I have forgotten to include it in the code that I have posted here.

Thanks for your nice suggestion. This is what I have tried when writing "SCR" function. I tried to write a simple code that call this function but still I cannot figure out what is the issue with this code.

Trying to run the simple code that calls "SCR" I came up with a solution for how to write "SCR" and it works well when I call it in another code. But the issue with this code is the weird result for the fitted curve when I call it as a fit type function. And this raised the question I have posted here:

https://www.mathworks.com/matlabcentral/answers/1729655-inconsistency-in-the-fitted-curve-and-the-plot-of-obtained-fitting-parameters?s_tid=srchtitle

I think the fit result is correct for that code but I don't undrestand why the fitted curve is not the same as the plot of the same function with the obtained fitting parameters.

Again, thank you!

[Matt J]

It is hard to comment without being sure what code you are actually working with. The code you have posted for us does not run:

nu0=1;
SCR(nu0,3,1,1,1)

Not enough input arguments.

Error in solution>SCR (line 10)
b = Start;

[Matt J]

Even when I give your code an input value for the missing "Start" Parameter and fix the missing declaration of gaussEqn1, it still does not produce an output at the given nu0, which I believe is what @Catalytic is talking about in his answer below.

nu0=1;
SCR(nu0,3,1,1,1,0)

Index exceeds the number of array elements. Index must not exceed 1.

Error in kkrebook2 (line 19)
deltaomega=omega(2)-omega(1);

Error in solution>SCR (line 22)
Rn(:) = kkrebook2(nu,z,0);

function p = SCR(nu,numGaussians,a,center,sigma,Start)

[gaussEqn1,thisGaussian]=deal(zeros(1,length(nu)));

for j = 1:length(nu)
    for i = 1:length(nu)
        for k = 1 : numGaussians
            if k==1   
                b = Start;
            else 
                b = Start + (center * (k-1));
            end
            thisGaussian(i) = a.*exp(-((nu(i)-b).^2)/(2.*(sigma.^2)));
            % Add into accumulator array:
            gaussEqn1(i) = gaussEqn1(i) + thisGaussian(i);
            z(i)= gaussEqn1(i);   
        end
    end
    Rn(:) = kkrebook2(nu,z,0);   
    hi2(j) = (4*pi*n.*nu(j)*L)/c;
    hi1(j)=(Rn(j));
    f(j)=(-r1+r2*exp(-(z(j))).*exp(-1i*hi1(j)).*exp(-1i*hi2(j)));
    g(j)=(1-r1*r2*exp(-(z(j))).*exp(-1i*hi1(j)).*exp(-1i*hi2(j)));
    h(j)= f(j)./g(j);
    m(j)= abs(h(j));
    p(j)= (m(j).^2);       
end
end

Answer 1

[Catalytic]

SCR has to be a 1D function of the independent variable nu, meaning it has to give valid output when nu is a scalar. That is not the case for your function.

Essentially, you appear to be fitting some N-dimensional surface function , given only a single sample, .

[Matt J]

Further to what @Catalytic points out, I've done a small test below using a version of your code from previous threads. Notice that at nu0=3.778651600000000e+14, a different value of SCR(nu0) is obtained depending on what other values of nu accompany nu0 as input. That should not happen for a 1D function. A curve sample f(x0) should be completely defined by the input value x0, not any other x.

format long

a =0.228738094823143;
center = 2.428378513756761e+07;
sigma = 3.833690154350237e+06;
d0 =0.047400255501257;

fun=@(nu) SCR(nu,9,a,center,sigma,d0);

nuA=linspace(3.778651600000000e14 ,  3.778653400000000e14,100);
nuB=linspace(3.778651600000000e14 ,  3.778653400000000e14,200);

scrA=fun(nuA); 
scrB=fun(nuB);

nuA(1)

ans = 
     3.778651600000000e+14

scrA(1)

ans = 
   0.296831567474829

nuB(1)

ans = 
     3.778651600000000e+14

scrB(1)

ans = 
   0.330814936483596

function p = SCR(nu,numGaussians,a,center,sigma,d0)

r1= 0.667746135604857;
r2= 0.998500100870347;
L= 0.00400002911721321;
c= 2.99792458E8;
n = 1.79997756058945;

StartComb = 377865161049434;

p = zeros(size(nu));
d = zeros(size(nu));
gaussEqn1 = zeros(1,length(nu));
thisGaussian = zeros(1,length(nu));
% alphad = zeros(1,length(nu));
% T = zeros(1,length(nu)); 
% Imag_n = zeros(1,length(nu)); 
% Real_n = zeros(1,length(nu)); 
phi = zeros(1,length(nu));
f = zeros(1,length(nu));
g = zeros(1,length(nu));
h = zeros(1,length(nu));
m = zeros(1,length(nu));

for i = 1:length(nu)
    
    for k = 1 : numGaussians
        
        if k==1 
            
        b = StartComb;
        
        else 
    
        b = StartComb + (center * (k-1));
        
        end
   
        thisGaussian(i) = a.*exp(-((nu(i)-b).^2)/(2.*(sigma.^2)));
        
        % Add into accumulator array:
        gaussEqn1(i) = gaussEqn1(i) + thisGaussian(i);

        d(i)= gaussEqn1(i);
        
    end

 end
       
        alphad = d/L;
        T = exp(-alphad);
        Imag_n =(-log(T)*c)./(nu*(4*pi));  %imaginary index of reffraction
        Real_n=kkrebook2(nu,Imag_n,0) + n; 

for i = 1:length(nu)

        phi(i)=(4*pi*Real_n(i).*nu(i)*L)/c;
        
        Finesse=center/(sigma*sqrt(8*log(2)));

        f(i)=(-r1+r2*exp(-d0)*exp(-(d(i)/Finesse)).*exp(-1i*phi(i)));
        g(i)=(1-r1*r2*exp(-d0)*exp(-(d(i)/Finesse)).*exp(-1i*phi(i)));
        h(i)= f(i)./g(i);
        m(i)= abs(h(i));
        p(i)= (m(i).^2);
   
        
end

 p=abs(p);

end

[Shaily_T]

Thanks for your response! If you omit the part related to "kkrebook2" in my code, you can get a result even with one nu.

nu0=1;

SCR(nu0,3,1,1,1,0)

The reason the "SCR" output is dependent on what other values of nu accompany nu0 as input is that "kkrebook2" function get a vector of "nu" and so its output in each "nu" is dependent on what was the input "nu" vector for that. This is the whole spirit of my question that how can I call "kkrebook2" such that I can use "SCR" as a fit type function. I have came up with a solution for how to write "SCR" and it works well when I call it in a code (The same as what you have tried for in this comment). But the issue with this code is the weird result for the fitted curve when I call it as a fit type function. And this raised the question I have posted here:

https://www.mathworks.com/matlabcentral/answers/1729655-inconsistency-in-the-fitted-curve-and-the-plot-of-obtained-fitting-parameters?s_tid=srchtitle

I think with your explanations here now I undrestand why the fitted curve in the above link doesn't make sense (because for fit it needs one scalar for each nu and it should not depend on the input vector of "nu"). But I am wondering is there a way that I can use "kkrebook2" in a fit type function ("SCR") and get reasonable results.

Thank you!

[Catalytic]

In order to answer that, you have to first explain to us if and why SCR should be considered a curve. A curve is a function of a single independent variable, whereas your SCR function is a function of many independent variables.

[Shaily_T]

My "SCR" function inputs are a vector of "nu" which is the independent variable and some scalar parameters like numGaussians, a center, and sigma, and the output of "SCR" is a vector. I think for fitting a function to data we have a vector of the dependent variable and a vector of the independent variable for data. And so the function to be fitted should also get a vector of independent variables and give back a vector of dependent variables. The way I have written it seems to do the same but when I call "SCR" it doesn't work.

[Catalytic]

"I think for fitting a function to data we have a vector of the dependent variable and a vector of the independent variable for data."

That's true in general, but the Curve Fitting Toolbox doesn't support general N-dimensional fitting. The Curve Fitting Toolbox only supports the fitting of functions with a 1-dimensional or 2-dimensional domain, whereas your function has an N-dimensional domain.

To put it more formally, if your code implements a mapping y=f(x): , then for f() to be considered a 1D curve, each y(i) can depend only on the corresponding x(i). However, in your SCR function each y(i) depends on all of the x(j) simultaneously.

For more general function fitting problems, you could look at lsqnonlin.