How can I avoid nested for loops when varying multiple array elements?

I am trying to do a fit to 9 sets of data at once. They all have the same fit parameters with the exception of tau and n(k) which are different for each set. I am trying to find the best values for A, which is an array of four numbers. However, the only way I can think of to do this is to loop through the four values using nested loops (as I don't know which combinations of numbers will give the best fit). This is taking a very long time, so I would like to avoid the for loops (or at least the nested loops). I don't have much experience programming so this is proving to be a challenge. Any suggestions would be great!

The parts of the code I think are relevant are below:

% initialize parameters used in fitting
dbstop if error
tau = [4E-5 5E-5 7.5E-5 1E-4 1.5E-4 2E-4 3E-4 4E-4 5E-4];    
tc = [2E-5 3E-4 6E-3 1E-1];               
delM2 = 1E8*[.9 1 .2 .3];   
n = [41 41 41 31 21 16 11 8 7]; 

chisqr_B = [999 999 999 999 999 999 999 999 999];
A_fit = zeros(1,4);
comp_fit = zeros(41,9);
T2eAve = (T2elow + T2ehigh)/2;
impr = 0;

for T2eInv = T2eAve:1:T2ehigh

% test fit for different population amplitudes A
for w = 0.06:0.01:0.40
    for x = 0.06:0.01:0.40
        for y = 0.06:0.01:0.40
            for z = 0.06:0.01:0.40
                A = [w x y z];
                amp0 = A(1)+A(2)+A(3)+A(4); % fit amplitude for n=0
                chisqr = zeros(1,9);
                amp = zeros(9,41);
                comp = zeros(41,9);
                chi = zeros(9,41);    

                  for k = 1:9

                      amp(k,1) = amp0;   

                      % sum ampl contributions from 4 motions (i)
                      for m = 2:n(k)  
                          amp1 = A(1)*something;
                          amp2 = A(2)*something;
                          amp3 = A(3)*something;
                          amp4 = A(4)*something;
                          amp(k,m) = amp1 + amp2 + amp3 + amp4;
                          comp(m,k) = log(amp(k,m)/amp(1,2)); 
                          if m>2
                              chi(k,m) = abs(comp(m,k)-dat(m,2*k))^2/abs(comp(m,k));
                              chisqr(1,k) = chisqr(1,k) + chi(k,m);
                          end
                      end

                  end

                  % store fit and population ampl for smallest chi-squared
                  if (sum(chisqr) <= sum(chisqr_B))
                      chisqr_B = chisqr;
                      A_fit = A;
                      T2e_fit = T2eInv;
                      comp_fit = comp;
                      impr = 1;
                  end
              end
          end
      end
  end

Note: I don't want to get rid of the k = 1:9 or the m = 1:n(k) loops, so I say A(i)*something as the 'something' is quite long and not really relevant. Also T2e values are input to the function.

Cheers, Lauren

Answer 1

[Bjorn Gustavsson]

You should have a look at:

fminsearch

that is the first optimization function to use for problems like this. That way you'd have to write the contenst in the 2 inner-most loops in a function returning the chi-2 value for a given combination of parameters and data set:

   your_chi2fcn(par,data_set)

Then you call fminsearch with that function:

bestPars = fminsearch(@(par) your_chi2fcn(par,data_setNr1),par0);

Where par0 is some start-guess you have.

HTH,

[Lauren]

Thanks for your quick reply! This seems like the best option, but I'm getting an error. If my 'par' is an array, will this not work? I have:

function[chisqr] = chi2func(A, dat)

where chisqr is size (1,9) and A is [ A(1) A(2) A(3) A(4) ]. Then using fminsearch I have:

function[A] = myfitfun(dat)
A = fminsearch(@(A) chi2func(A,dat), [.4 .3 .3 .1]);

The error is in fminsearch at 205:

??? Subscripted assignment dimension mismatch.
fv(:,1) = funfcn(x,varargin{:});

I'm not sure if this is an issue with chisqr or with A (or maybe something else)

[Bjorn Gustavsson]

The way I thought you'd want to go about it was to sum up all the chisqr like you did in the if-condition in the loopy version above. I think that is what's causing you problems, fminsearch needs a scalar output from the objective function.

[Lauren]

I completely missed that part. It seems to be working fine now! Thanks a bunch!