How to loop a function though multiple datasets

Hello everyone,

I have a bit of a basic question, and I've tried reading other answers on this, but they each seem to have their own specific problem (nor have I been able to successfully implement their solutions to my own issues).

In a nutshell, I have successfully been able to use lsqccurvefit for one of my datasets. I am now trying to create a loop, so I can run the same function on multiple datasets. My code is this:

Data1=[0.596421471	0.06
0.5859375	0.119284294
0.566037736	0.29296875
0.530035336	0.622641509
0.418994413	1.219081272
0.388619855	3.058659218
];
Data2=[5.00E+04	3.82E+04	3.45E+04	3.42E+04	3.74E+04	3.21E+04	2.81E+04	2.88E+04
5.00E+04	3.82E+04	3.45E+04	3.42E+04	3.74E+04	3.21E+04	2.81E+04	2.88E+04
3.08E+04	3.07E+04	2.19E+04	2.23E+04	2.53E+04	2.05E+04	1.98E+04	1.89E+04
2.05E+04	1.87E+04	1.30E+04	1.10E+04	1.62E+04	1.31E+04	1.05E+04	1.05E+04
8963	1.11E+04	6243	3504	6454	4346	4448	4357
0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00
];
Pt=Data1(:,1);
Lt=Data1(:,2);
Int=Data2(:,1);
plot(Lt,Int, 'ro');
xdata=[Pt,Lt]
F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(2))-((xdata(:,1)+xdata(:,2)+x(2)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1)/2.*xdata(:,1))
x0=[1 1]
options=optimoptions('lsqcurvefit','MaxFunEvals',5e2)
lb=-Inf
up=+Inf
[x,resnorm,~,exitflag,output]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
hold on 
plot (Lt, F(x,xdata))
hold off

Data 1 is constant and does not change, it is what my function is using. Data2 is my data. I am also plotting my data and the output of the function so I can see how well they match (which has so far been extremely poor matching). What I want to do now is have my [Int] cycle through each column independently, so for the first round it runs the lsqc through column ;1, 2nd round through column ;2, etc. At the same time, I would also like to individuallly (or combined) plot each individual column with the output of the function (for each run it goes through). I know I need to use the "for" function, but I don't quite know how to cycle through each column individually using the Int that I've already defined.

Finally, I had gotten this to work using one column set, but now that I've inputed multiple column sets, I seem to have broken my script. I now get the error

Index in position 1 is invalid. Array indices must be positive integers or logical
values.

Error in Testscript4loop (line 18)
plot(Lt,Int, 'ro');

Which is strange, since before I added all the extra data (Data2 used to be only 1 column), it worked fine.

[Adam]

for n = 1:size( Data2 )
   Int=Data2(:,n);
  ...
end

should work, depending how many plots you want at once and what results you may want to store for outside the loop (just put these into an array indexed by n).

The error you report doesn't seem to be in any way related though and looks like the kind of error you would get if you have declared a variable named plot which is hiding the function of the same name (never do this!).

[Sam Mahdi]

Thank you, it worked perfectly. Although, 3 things from this:

1) Could you explain exactly what's going on here? I'm trying to understand the process.

It appears when you state n=1:aize ( Data2 ), you're telling it to start with the first column in Data2. The next line Int=Data2(:,n); appears to be defining Int (as I had before), but now defining it to be [n] number of columns of Data2 (previously I had it as Int=(:,1) so I only defined it as one). So is n in the 2nd line (Int+...) being defined by the n int he 1st line (n=1:size...)?

2) For plotting this, I'm a bit confused what you mean put it into an array indexed by n. For example, I have it written as

plot (Lt,Int, 'ro'),

but since Int is defined as an array already, do I have to change it? Shouldn't the plot change as Int changes? Shouldn't this apply to the 2nd plot as well? Since the fitting is dependent on Int?

3) Finally, is there a way to improve the fitting by telling it to predict the starting values based off the results of the first run? I don't have a lot of data points, but I think it would be more useful if instead of having the same function start from my initial values (x0=[1 1]), have it actually predict the initial values for the 2nd data set, based off the results from the first one. And then have this continue for each iteration, having the starting values converge hopefully closer and closer to the true values (I'm thinking the easiest way this can be done is using the results of the previous run to be the initial values for the current run).

[Adam]

Actually

1:size( Data2 )

was a typo and should read

1:size( Data2, 2 )

so I'm slightly surprised that works as it was!

Putting this in a loop will tell it to run the loop as many times as there are columns in Data2.

doc for

should help you understand for loops.

Int=Data2(:,n);

will take the nth column of Data2 and because it is in a loop defined as above n will increment from 1 up to the number of columns of Data2 each time round the loop.

Putting results in an array indexed by n depends on whether you want to keep the results from a loop outside of the loop (if not they will be overwritten next time round the loop). If you don't then your plot instruction can stay the same (though you will need a hold on instruction to ensure next time round the loop your plot is added). If you do want to keep the results then e.g.

Int(:,n) = Data2(:,n)

will do this. Though you should predefine Int as:

Int = zeros( size( Data2 ) );

In the case of Int this is pointless though as you would just end up with a copy of Data2, but later results can be stored in this way.

I don't know anything about lsqcurvefit though so I can't help with improving its results.

[Sam Mahdi]

Thank you for the speedy detailed reply!

A few things:

1) I'm curious as to why the 2 needs to be included? My run was already running through each column in Data 2, so I don't see what that changed. Furthermore, why the number 2 specifically. I.E. Will it be different if I put 1 or 3?

2) I tried those commands, put it just popped up errors saying YDATA dimensions were off. I do have a hold in my script

hold on 
plot (Lt, F(x,xdata))
plot (Lt, F(x2,xdata))
hold off
end

Ignore (Ignore F(x2) that's another algorithm I'm testing). I don't have a hold on the initial plot using Lt and Int simply because the plot for that shouldn't change much. So I don't know if this is overwriting all my previous results and only displaying the plot of the last run (as you were saying), or if all the plots simply overlap perfectly.

[Adam]

doc size

tells you about the arguments to the size function and what it returns. The documentation, combined with the command line is by far your best resource for understanding how functions work. Get used to using both and it will help your learning immeasurably. I've been using Matlab for 13 years now, but I still always have multiple tabs open in the help as I am always referring to something or other and am constantly trying things out on the command line with test matrices or whatever other inputs I need. For example:

>> a = rand( 3, 4 );
>> size( a )

ans =

     3     4

>> size( a, 1 )

ans =

     3

>> size( a, 2 )

ans =

     4

Very simple to test in under a minute and helps you understand what results are giving.

Specifically, here if you don't give a second argument to the size function you get an array of results, giving the size for each dimension (rows first, then columns). Passing '2' in says 'give me the size of dimension 2 - i.e. columns'.

I've never really worked out what a for loop or if statement does if you give it something unexpected like a vector when it wants a scalar - I just make sure not to do it!

But in this case again you can test very quickly and see :

>> 1:size( a )

ans =

     1     2     3

which is to be expected really. It just takes the first output of the size function (the number of rows), so unless you happen to have the same number of rows as columns then your for loop would either crash or miss some data out doing this while indexing columns.

As for your plotting, I can't work out, just at a glance, (and especially not seeing the updated code with the loop), what the sizes are of x and xdata when you plot them, but again, a breakpoint on the relevant line and

size( x )
size( xdata )

on command line will soon tell you this. They need to both be vectors, of equal length, in order to plot.

[Sam Mahdi]

I understand thank you. As for the plotting, I've just been using the same as before, but using hold to include all the runs.

for n = 1:size( Data2, 2 )
   Int=Data2(:,n);
  ...
hold on
plot(Lt,Int, 'ro');
xdata=[Pt,Lt]
F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(2))-((xdata(:,1)+xdata(:,2)+x(2)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1)/2.*xdata(:,1))
x0=[1 1]
options=optimoptions('lsqcurvefit','MaxFunEvals',5e2)
lb=[];
up=[]
[x,resnorm,~,exitflag,output]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [];
ub = [];
[x2,resnorm2,~,exitflag2,output2]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
M=mean(x) 
M=mean(x2) 
plot (Lt, F(x,xdata))
plot (Lt, F(x2,xdata))
hold off
end

So this plots both the fit and the data for each run. This is why I was confused as to your earlier post because I don't need to redefine anything, I can just move the hold earlier and inside the loop. The only problem I'm having is I have a lot of data right now, so I'm using the mean function to just take an average and give me a singular value. I'm trying to plot that singular value, so I've done something like plot (Lt, M), and assumed this would plot the mean of x, but doesn't appear to be that straightforward.

[Adam]

Again I've lost track of exactly what dimensionality various results are, but the mean function also has an (optional) dimension argument to tell it which dimension to calculate the mean over. It may be that you are doing it long the wrong dimension and may need

M = mean( x, 2 )

instead.

[Sam Mahdi]

It calculates the mean properly, I just want to find a way to plot the mean. I.E. Ideally I'd like to remove plot (Lt, F(x,xdata)) and the 2nd plot as well, and simply have it plot the mean the mean (I.E. plot (Lt, M) or something along those lines). This way instead of getting 5 different lines on my plot, I get a singular line representing the mean values of my function.

To do this, I defined M2 to to be the average of my function, and then tried something like

M2=mean(F(x,xdata))
plot (Lt, M2)

While this calculates the mean of my function well, it doesn't actually plot anything. In a nuthsell, instead of graphing the values of my [for] loop, I'd like to graph the combined mean value of the values of that loop.

[Adam]

So you will need to store whatever data is relevant in the loop in an array instead of plotting it. As it stands you will end up with a mean value per column of your data, each of which gets thrown away apart from the one for the final column.

[Sam Mahdi]

I see, yes you're right. All M is doing is taking the average of the 2 x2s that I'm solving (x(1) and x(2). It is not actually averaging my loop. I initially thought the difference betwee M and my output was because it was averaging, but it's only different because I have 2 variables, and it's only averaging the values of the 2 variables from each individual run.