How to fit 6 curves simultaneously to solve for 2 unknowns?

y=(1-exp(-0.5*Phi*Gain*(x-TimeDelay)^2))*Amax

Gain and TimeDelay are unknown. x is time (T1-T6) and y is response (F1-F6)(source of noise). Phi is different for each curve, and Amax is constant.

Attached is a sample dataset. The values for Phi and the Amax for this dataset are found to the right of the data.

I tried the Curve Fitting Toolbox, but it seems to be designed for one trace at a time. I can get the curves to fit for the first 4 traces but not the last 2. However, I need to get one gain value and one TimeDelay that is the best (or least bad) fit across all traces. I don't know what function to use.

Thank you! Any help is appreciated.

***note: I am editing the question with the comments from Walter taken into consideration. I initially typed the equation incorrectly (wrote Amax was negative and it isn't) and I also thought TimeDelay was constant and it isn't.

[Walter Roberson]

The dataset didn't make it to attachment :(

[Delores Davis]

My bad. I copied the data to an excel file and then forgot to attach it. It's here now.

Answer 1

[Walter Roberson]

Solve for Gain to get

Gain = -2*ln((Amax+y)/Amax)/(Phi*(x-TimeDelay)^2)

then make a single table of all of the values over all of the datasets, and run a least-squared fit. Looks like the solution to that would just be the mean of those values.

Per-variable confidence bounds is not as easy to compute. If I understand your setup correctly, "y" should be treated as having noise, but not the others?

[Delores Davis]

You are correct, Y is the noisy variable. When you say to make a single table of all the values, do you mean like the sample dataset that I attached, or a table of gain values? How would I plot that and how do I run a least-squared fit? Do you have any idea how to calculate the confidence interval?

[Walter Roberson]

Each of your columns has an associated Phi. For each column use the associated Phi and the two constants and the x and y data and the above formula to calculate estimated Gain, resulting in a column of Gain for each of your input columns. Then create a single column by concatenating all of those together into one big list. With that on hand, mean() of it all gives "best" gain, at least under least-squares model.

To get some confidence levels you can use polyfit() of all the y for one coordinate (concatenated together, same as would be used in calculating the master list of Gain), with the master list of Gain as the dependent variable, and ask for 1 as the degree, and ask for the optional output arguments mu and S. Those give you information about the confidence bounds.

Your Gain depends on x in the form 1/(x-TimeDelay)^2 which implies it goes to infinity as x approaches the time delay. If you have x values near TimeDelay then the contribution in the list of Gain is going to be high, leading to a high fitted gain, and the effect is going to be stronger the closer x gets to TimeDelay. Your TimeDelay is 0.005 and you do have some values in the range of 0.0056, so you are going to get big spikes in value that are going to mess up the calculation of best Gain by mean. This should lead you to distrust the calculated values.

As I look at your table of values it looks to me as if it does not support your model. Your model predicts that for any fixed Gain, as x approaches TimeDelay that the y value should go to 0, but there is no sign of that in your data.

[Delores Davis]

I just found out a new fact about this equation: TimeDelay is a constant FOR EACH ANIMAL but it varies slightly with a range between 4-10 ms. The curve fit equation is supposed to give this value to me.

As far as the model, it is widely accepted in my field for interpretation of this data type and has been published multiple times.

I attached a truncated dataset from the flash point to an inflection point around the middle of the maximum response. This is the area where a good fit is mandatory. I was told by a coworker that if I include data past the point where I want to fit, Matlab will try to fit the other data too and skew the results (and I don't know if that is true).

So keeping in mind that I can't plug in a number for TimeDelay, how do I approach this equation?

Should I send an unaltered dataset?

[Walter Roberson]

To check:

Is the question now to find a single Gain for all of the experiments together, and also find the per-group TimeDelay?

[Walter Roberson]

If we examine the formula we see that the minimum for any curve with fixed Gain and TimeDelay occurs when x = TimeDelay, and that the value is larger on either side of TimeDelay.

As we examine the data you attached above, the F1 goes to 0 at 0.0041 (4.1 ms), which fits with your 4 ms - 10 ms figure. F2 goes to 0 at 0.0011, 1.1 ms, significantly below that time range. F3 an F4 go to 0 at 0.0001, 0.1 ms, not even close to the range you indicated. F3 to F6 are strictly increasing and so if they fit the model then the implication would be that the TimeDelay is no greater than 0.1 ms for them.

There are three ways in which the model could be rescued:

• T1 might be biased by at least 4 ms (0.004), so x = T1 + 0.004 (or some other constant... or does it differ per dataset?) instead of x = T1 . This would make the smallest x 0.0041 which together with a TimeDelay of 0.0040 would allow a non-zero response "before" the data shown
• We could assume that the "noise" in y is potentially at least as large as 85.5, to allow the "true" y for F6 time 0.0041 to be "corrected" to 0 to get the 0 entry we need at x = TimeDelay and minimum 0.004 TimeDelay. This would require building a noise model for y that would have rather extraordinary characteristics, perhaps time-based.
• We could postulate that you attached the wrong data

[Walter Roberson]

Another possibility: you gave the wrong formula.

Consider the first portion

(1-exp(-0.5*Phi*Gain*(x-TimeDelay)^2))

that is non-negative until exp(-0.5*Phi*Gain*(x-TimeDelay)^2) > 1. Log both sides and that is -0.5*Phi*Gain*(x-TimeDelay)^2 > 0 . But Phi and Gain are positive values and 0.5 is a positive multiplier so that reduces to -(x-TimeDelay)^2 > 0 which reduces to (x-TimeDelay)^2 < 0 . With x being positive the only way that can happen is if TimeDelay is complex valued, which is known not to be the case. Therefore -0.5*Phi*Gain*(x-TimeDelay)^2 <= 0 and so exp(-0.5*Phi*Gain*(x-TimeDelay)^2) <= 1 and so (1-exp(-0.5*Phi*Gain*(x-TimeDelay)^2)) >= 0

We thus have a non-negative quantity being multiplied by -Amax , where Amax is a positive quantity and so -Amax is negative. non-negative multiplied by negative non-positive, so the formula gives non-positive values.

We now examine the F1, F2, F3, F4 values and see that they are non-negative. If we do not assume that the data is wrong, then the other possibility is that the formula is wrong. For example it might be

(1-exp(-0.5*Phi*Gain*(x-TimeDelay)^2)) * Amax

with no negative sign before Amax. That would at least allow us to get positive values out...

[Delores Davis]

Amax in the equation is positive. (My bad, I am so sorry)

The Gain and the TimeDelay should be treated equally in the equation, as far as we need the one best value of each (with confidence bounds) to fit the data across all runs F1 thru F6. The entire dataset given is from the same animal, (ERG a-waves from low flash intensity to high F1->F6).

I am not getting e-mail notifications when you leave a comment, I checked the box, but I'm getting nothing. Is there a way to fix that?

And kudos for noticing my incorrect sign.

[Delores Davis]

Using the curve fitting toolbox, I can get F1, F2, F3, and F4 to fit to the equation but not F5 and F6. I've tried manipulating the StartPoint of the fit on these traces and I get the error:

'Warning: A negative R-square is possible if the model does not contain a constant term and the fit is poor (worse than just fitting the mean). Try changing the model or using a different StartPoint.'

The values for TimeDelay and Gain are nowhere close to each other for each trace (but honestly I don't know if that is expected or not) and I don't know how to relate the traces and/or outputs to each other to get one value.

Do you by chance understand IgorPro's GlobalFit tool? It plots all of the traces, shows the fit for each trace, gives TimeDelay and Gain for each trace, and also shows a TimeDelay and Gain for the entire dataset with Confidence Intervals all simultaneously. There is a feature to link variables and from what I recall, all the variables are selected to be linked to each other, but I am scared of the results because it asks for your guess for the TimeDelay and Gain...and the output changes considerably depending on how you guess.

We have one license on one computer in lab for IgorPro and since I am the first user of a different instrument for ERG data collection, I want to change the method for data analysis to something less dependent on 'guesswork'.

[Walter Roberson]

Sorry, I am not familiar with IgorPro at all.

[Delores Davis]

My main point was that if Igor can do it, Matlab can too. I just need to know how.

[Walter Roberson]

Meaningful documentation about IgorPro appears to only be available to people who have paid for IgorPro, so I have to be satisfied with trying to invent mathematics to handle the situation.

[Walter Roberson]

Sometimes in data fitting, what you need to do is not look for parameters that give you a "good" fit for everything, but instead look for parameters that give you a least equally bad fit for everything

[Delores Davis]

I wrote this code based on answers to previous questions from other users in the forum.

syms g1 t1;

solve((T1==(1-exp(-0.5*38.01*g1*(F1-t1).^2))*453),t1);

syms g2 t2;
solve((T2==(1-exp(-0.5*114.03*g2*(F2-t2).^2))*453), t2);

syms g3 t3;

solve((T3==(1-exp(-0.5*380.1*g3*(F3-t3).^2))*453), t3);
syms g4 t4;
solve((T4==(1-exp(-0.5*1140.3*g4*(F4-t4).^2))*453), t4);

syms g5 t5;

solve((T5==(1-exp(-0.5*3801*g5*(F5-t5).^2))*453), t5);

syms g6 t6;
solve((T6==(1-exp(-0.5*9502.5*g6*(F6-t6).^2))*453), t6);

g1=(log(1-(T1/453))/(-0.5*38.01*(F1-t1).^2));

g2=(log(1-(T2/453))/(-0.5*114.03*(F2-t2).^2));

g3=(log(1-(T3/453))/(-0.5*380.1*(F3-t3).^2));

g4=(log(1-(T4/453))/(-0.5*1140.3*(F4-t4).^2));

g5=(log(1-(T5/453))/(-0.5*3801*(F5-t5).^2));

g6=(log(1-(T6/453))/(-0.5*9502.5*(F6-t6).^2));

hold on

scatter (T1,F1)

scatter (T2,F2)

scatter (T3,F3)

scatter (T4,F4)

scatter (T5,F5)

scatter (T6,F6)

plot (T1,((1-exp(-0.5*38.01*g1*(F1-t1).^2))*453))

plot (T2,((1-exp(-0.5*114.03*g2*(F2-t2).^2))*453))

plot (T3,((1-exp(-0.5*380.1*g3*(F3-t3).^2))*453))

plot (T4,((1-exp(-0.5*1140.3*g4*(F4-t4).^2))*453))

plot (T5,((1-exp(-0.5*3801*g5*(F5-t5).^2))*453))

plot (T6,((1-exp(-0.5*9502.5*g6*(F6-t6).^2))*453))

I get the following errors:

    Warning: Cannot find explicit solution. 
    > In solve at 319
      In PLOT_IT_ALL at 3 
    Warning: Cannot find explicit solution. 
    > In solve at 319
      In PLOT_IT_ALL at 5 
    Warning: Cannot find explicit solution. 
    > In solve at 319
      In PLOT_IT_ALL at 7 
    Warning: Cannot find explicit solution. 
    > In solve at 319
      In PLOT_IT_ALL at 9 
    Warning: Cannot find explicit solution. 
    > In solve at 319
      In PLOT_IT_ALL at 11 
    Warning: Cannot find explicit solution. 
    > In solve at 319
      In PLOT_IT_ALL at 13 
    Warning: The system is rank-deficient. Solution is not unique. 
    > In C:\Program Files\MATLAB\R2014b\toolbox\symbolic\symbolic\symengine.p>symengine at 56
      In sym.sym>sym.privBinaryOp at 835
      In sym.sym>sym.mrdivide at 272
      In PLOT_IT_ALL at 19 
    Warning: The system is rank-deficient. Solution is not unique. 
    > In C:\Program Files\MATLAB\R2014b\toolbox\symbolic\symbolic\symengine.p>symengine at 56
      In sym.sym>sym.privBinaryOp at 835
      In sym.sym>sym.mrdivide at 272
      In PLOT_IT_ALL at 20 
    Warning: The system is rank-deficient. Solution is not unique. 
    > In C:\Program Files\MATLAB\R2014b\toolbox\symbolic\symbolic\symengine.p>symengine at 56
      In sym.sym>sym.privBinaryOp at 835
      In sym.sym>sym.mrdivide at 272
      In PLOT_IT_ALL at 21 
    Warning: The system is rank-deficient. Solution is not unique. 
    > In C:\Program Files\MATLAB\R2014b\toolbox\symbolic\symbolic\symengine.p>symengine at 56
      In sym.sym>sym.privBinaryOp at 835
      In sym.sym>sym.mrdivide at 272
      In PLOT_IT_ALL at 22 
    Warning: The system is rank-deficient. Solution is not unique. 
    > In C:\Program Files\MATLAB\R2014b\toolbox\symbolic\symbolic\symengine.p>symengine at 56
      In sym.sym>sym.privBinaryOp at 835
      In sym.sym>sym.mrdivide at 272
      In PLOT_IT_ALL at 23 
    Warning: The system is rank-deficient. Solution is not unique. 
    > In C:\Program Files\MATLAB\R2014b\toolbox\symbolic\symbolic\symengine.p>symengine at 56
      In sym.sym>sym.privBinaryOp at 835
      In sym.sym>sym.mrdivide at 272
      In PLOT_IT_ALL at 24 
    Error using plot
    A numeric or double convertible argument is expected

    Error in PLOT_IT_ALL (line 33)
    plot (T1,((1-exp(-0.5*38.01*g1*(F1-t1).^2))*453))

:(

[Walter Roberson]

If your T* and F* are vectors then solve() will not be able to find a solution as you are asking it to solve a vector (probably longer than 2 elements) of equations in two variables. solve() is going to try to find the g1 and t1 that solve all of the elements simultaneously. solve() is not going to look for a least-squares solution or anything similar: it is going to try to solve all of the equations simultaneously, and is almost certain to fail in doing so.

[Delores Davis]

Do you know what function or command to use? I am now thinking nlinfit, but I don't understand how to 'word' the command for my purposes. Looking up tutorials and examples, would be good to know if this is the right command before I get too deep into the wrong rabbit hole.

[Walter Roberson]

For all of the data considered together, the optimal gain is between 12.4148263515933 and 12.4148299319728, and the optimal time delay is between 0.00214757638888889 and 0.00214757725694444

This is not to say that the fit for any of the individual curves is good, only that these values minimize the sum of the squares of y minus predicted y.

Method for this particular solution: put all the data together, constructed a symbolic expression of the sum of squares of y - predicted y, did a 3D plot, kept refining the area of interest.

There is an alternate method that can get pretty close; I will write about it later, but I need to sleep now.

[Delores Davis]

Are the two values given above the 95% confidence bounds?

When you 'put all the data together' is it the response data in a list then the times in a list? Essentially, what I want to know is did you filter anything out?

What is the predicted y? Is it the actual response or a mathematically derived y? If derived? How??

When you do a 3D plot, are you using the curve fitting toolbox? I don't understand what the 3 dimensions are. x=time, y=response, z=???

Just to clarify the sum of squares of y is this for all of the data together or each group separately?

Double checking the code for the sum of squares calculation:

F1_squared=F1.^2

sum_F1_squared=sum(F1_squared)

then, if I understand you correctly next would be something like... sum of squares - response

delta_F1=sum_F1_squared - F1

Would delta_F1 be Z in the 3D plot? How would I integrate delta_F1 into the original equation?

How do you refine an area of interest?

By the way, thank you for your help.

[Walter Roberson]

For each of the experiments, construct a table of x (time), y, and repeat the experiment's phi for each entry, so three columns. Those three items together now contain all of the "identity" of the samples. Do the same thing for each of the other experiments. Now do a vertical concatenation of all of the data, making a single large table

x1_1 y1_1 phi1
x1_2 y1_2 phi1
x1_3 y1_3 phi1
...
x1_25 y1_25 phi1
x2_1 y2_1 phi2
x2_2 y2_2 phi2
...
x2_25 y2_25 phi2
x3_1 y3_1 phi3
x3_2 y3_2 phi3
...
x3_24 y3_24 phi3

and so on, so that all of the entries are together in one table. There will be 126 rows based on the information you posted

Now using symbolic Gain and TimeDelay, evaluate

 (y) - ((1-exp(-0.5*Phi*Gain*(x-TimeDelay)^2))*Amax)

for each of the table entries, getting one result for each row that is the difference between actual y and the y that would be predicted by that x and phi value with symbolic gain and timedelay.

Square each of those differences. Then add them all together to produce one long formula that has a bunch of constants (the experimental values) and symbolic gain and timedelay all over the place. Call this formula F. Our goal in fitting is to minimize the sum of squares of the differences between actual and prediction, so that translates to minimizing that long single formula, F, in two variables.

As F is now a formula in two variables, you can do a 3D plot, along the lines of

Fn = matlabFunction(F,'vars', [Gain, TimeDelay]);  %convert symbolic to anonymous function

gain = linspace(0,40,100);
timedelay = linspace(0,.01,100);
[Gain_grid, td_grid] = ndgrid(gain, timedelay);
predict_diff = Fn(Gain_grid, td_grid);
surf(Gain_grid, td_grid, predict_diff)

and now you can visually look to see the shape and figure out the useful ranges.

At the first few scales you might be skipping over some important features, but once you have moved in a bit, you can use

[mindiff, minidx] = min(predict_diff);
[r, c] = ind2sub(size(Gain_grid), minidx);

and then the "real" minimum should occur somewhere between gain(r-1) and gain(r+1), and timedelay(c-1) and timedelay(c+1) . You can use this and a few iterations to narrow down the minimum positions.

I have an appointment now, but I can post the extracted data later.

[Delores Davis]

Made the table with column one corresponding to response, column 2 corresponding to time and column 3 corresponding to Phi.

Then I entered the following code:

sym Gain 
sym TimeDelay
y=All_data(:,1);
x=All_data(:,2);
Phi=All_data(:,3);
Amax=453;
y_p=y-((1-exp(-0.5*Phi*Gain*(x-TimeDelay)*(x-TimeDelay)))*Amax);

And got the following errors:

>> y_predicted

ans =

Gain

ans =

TimeDelay

Error using mupadmex
Error in MuPAD command: The dimensions do not match. [(Dom::Matrix(Dom::ExpressionField()))::_mult2]

Error in sym/privBinaryOp (line 835)
            Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in  *  (line 219)
        X = privBinaryOp(A, B, 'symobj::mtimes');

Error in y_predicted (line 7)
y_p=y-((1-exp(-0.5*Phi*Gain*(x-TimeDelay)*(x-TimeDelay)))*Amax);

This is as far as I was able to go. What did I do wrong?

[Walter Roberson]

You might have to code with .* instead of *.

Also, you should use "syms" instead of "sym"

syms Gain TimeDelay

[Walter Roberson]

You might have to code with .* instead of *.

Also, you should use "syms" instead of "sym"

syms Gain TimeDelay

[Delores Davis]

Changing * to .* and using syms instead of sym I was able to make a virtual taco. I assume this is what the graph was supposed to look like.

I don't know how to move in on an area of interest. I am assuming it would be in the area where pred_diff is the smallest. I changed the axis of the figure and ran the last bit of code and the output is a row of 100 zeros or 100 ones for r,c, mindiff, and minidx. I am assuming this is because refining the ROI is more involved than changing the axis range.

Also, I need a gain value and a time delay value with 95% confidence intervals, which I don't know how to create from the 'taco'.

[Walter Roberson]

To be sure we are working with the same code:

syms Gain TimeDelay
y=All_data(:,1);
x=All_data(:,2);
Phi=All_data(:,3);
Amax=453;
y_p=y-((1-exp(-0.5*Phi*Gain*(x-TimeDelay)*(x-TimeDelay)))*Amax);

F = sum((y_p).^2);  %objective function in symbolic form

Fn = matlabFunction(F,'vars', [Gain, TimeDelay]);  %convert symbolic to anonymous function

%starting point for search
lowgain = 0;
highgain = 30;
lowtd = 0;
hightd = 0.01;

refine_factor = 100;
refine_steps = 5;

for K = 1 : refine_steps
  gain = linspace(lowgain, highgain, refine_factor+1);
  timedelay = linspace(lowtd, hightd, refine_factor+1);
  [Gain_grid, td_grid] = ndgrid(gain, timedelay);
  predict_diff = Fn(Gain_grid, td_grid);
  surf(Gain_grid, td_grid, predict_diff);
  drawnow();
 [mindiff, minidx] = min(predict_diff(:));
 [r, c] = ind2sub(size(Gain_grid), minidx);
 bestgain = Gain_Grid(r,c);
 besttd = td_grid(r,c);
 fprintf('Step #%d, minimum value %g near (Gain = %g, TimeDelay = %g)\n', K, mindiff, bestgain, besttd);

   lowgain = Gain_grid(r-1,c);
   highgain = Gain_grid(r+1,c);
   lowtd = td_grid(r, c-1);
   hightd = td_grid(r, c+1);
 end

Note: the above is untested as I do not have the symbolic toolbox.

[Delores Davis]

Step #1, minimum value 120131 near (Gain = 14.7, TimeDelay = 0.0025)
Step #2, minimum value 120097 near (Gain = 14.586, TimeDelay = 0.002474)
Step #3, minimum value 120097 near (Gain = 14.592, TimeDelay = 0.00247468)
Attempted to access Gain_grid(102,68); index out of bounds because size(Gain_grid)=[101,101].

Error in Walter (line 29)
   highgain = Gain_grid(r+1,c);

I got this output from your code. I don't have the symbolic toolbox. Can I download it on the mathworks site?

[Walter Roberson]

You do have the symbolic toolbox; you are using "syms" which is provided by the symbolic toolbox.

[Walter Roberson]

I did encounter something like that myself at one point. As the function is not entirely smooth, it is possible that when it is examined at one scale the minimum is in a different region than if it is examined at a different scale, especially as you make your TimeDelay more refined. Let me see...

Change the

   lowgain = Gain_grid(r-1,c);
   highgain = Gain_grid(r+1,c);
   lowtd = td_grid(r, c-1);
   hightd = td_grid(r, c+1);

to

   if r > 1
     lowgain = Gain_grid(r-1, c);
   else
     lowgain = Gain_grid(1, c) - 2 * (Gain_grid(2,c) - Gain_grid(1,c));
     fprintf('Step #%d: widening Gain search lower\n', K);
   end
   if r <= refine_factor
     highgain = Gain_grid(r+1,c);
   else
     highgain = Gain_grid(r,c) + 2 * (Gain_grid(r,c) - Gain(Grid(r-1,c));
     fprintf('Step #%d: widening Gain search higher\n', K);
   end

   if c > 1
     lowtd = td_grid(r, c-1);
   else
     lowtd = td_grid(r,1) - 2 * (td_grid(r,2) - td_grid(r,1));
     fprintf('Step #%d: widening TimeDelay search lower\n', K);
   end
   if c <= refine_factor
     hightd = td_grid(r,c+1);
   else
     hightd = td_grid(r,c) + 2 * (td_grid(r,c) - td_grid(r,c-1));
     fprintf('Step #%d: widening TimeDelay search higher\n', K);
   end

Again untested.

In this code, when it is detected that the minima was along one of the sides, then I project by twice the current spacing rather than by just the current spacing, with the idea being that the detected narrow grove might be long.

[Delores Davis]

Here is the output. Kudos on fixing that problem!

 >> Walter
  Step #1, minimum value 120131 near (Gain = 14.7, TimeDelay = 0.0025)
  Step #2, minimum value 120097 near (Gain = 14.586, TimeDelay = 0.002474)
  Step #3, minimum value 120097 near (Gain = 14.592, TimeDelay = 0.00247468)
  Step #3: widening Gain search higher
  Step #4, minimum value 120097 near (Gain = 14.5922, TimeDelay = 0.00247472)
  Step #4: widening Gain search higher
  Step #5, minimum value 120097 near (Gain = 14.5922, TimeDelay = 0.00247472)
  Step #5: widening Gain search higher

How would I find 95% confidence bounds?

[Delores Davis]

Also, the symbolic toolbox doesn't appear under Apps, so I didn't realize I have it. I thought 'syms' was general Matlab terminology.

[Walter Roberson]

Posting the data here to make it easier to access:

All_data = [
  0.111e1 0.1e-3 0.3801e2;
  0.116e1 0.6e-3 0.3801e2;
  0.111e1 0.11e-2 0.3801e2; 
  0.114e1 0.16e-2 0.3801e2; 
  0.85e0 0.21e-2 0.3801e2; 
  0.62e0 0.26e-2 0.3801e2; 
  0.37e0 0.31e-2 0.3801e2; 
  0.1e-1 0.36e-2 0.3801e2; 
  0 0.41e-2 0.3801e2;
  0.9e-1 0.46e-2 0.3801e2; 
  0.67e0 0.51e-2 0.3801e2; 
  0.127e1 0.56e-2 0.3801e2; 
  0.231e1 0.61e-2 0.3801e2; 
  0.325e1 0.66e-2 0.3801e2; 
  0.484e1 0.71e-2 0.3801e2; 
  0.634e1 0.76e-2 0.3801e2; 
  0.861e1 0.81e-2 0.3801e2; 
  0.1093e2 0.86e-2 0.3801e2; 
  0.1354e2 0.91e-2 0.3801e2; 
  0.1621e2 0.96e-2 0.3801e2; 
  0.1925e2 0.101e-1 0.3801e2;
  0.2242e2 0.106e-1 0.3801e2; 
  0.2618e2 0.111e-1 0.3801e2; 
  0.2953e2 0.116e-1 0.3801e2; 
  0.3362e2 0.121e-1 0.3801e2; 
  0.15e0 0.1e-3 0.11403e3; 
  0.5e-1 0.6e-3 0.11403e3; 
  0 0.11e-2 0.11403e3; 
  0.6e-1 0.16e-2 0.11403e3; 
  0.32e0 0.21e-2 0.11403e3; 
  0.68e0 0.26e-2 0.11403e3; 
  0.132e1 0.31e-2 0.11403e3; 
  2 0.36e-2 0.11403e3; 
  0.301e1 0.41e-2 0.11403e3; 
  0.409e1 0.46e-2 0.11403e3; 
  0.567e1 0.51e-2 0.11403e3; 
  0.747e1 0.56e-2 0.11403e3; 
  0.1002e2 0.61e-2 0.11403e3; 
  0.1293e2 0.66e-2 0.11403e3; 
  0.1677e2 0.71e-2 0.11403e3; 
  21 0.76e-2 0.11403e3; 
  0.2617e2 0.81e-2 0.11403e3; 
  0.3157e2 0.86e-2 0.11403e3; 
  0.3765e2 0.91e-2 0.11403e3; 
  0.4407e2 0.96e-2 0.11403e3; 
  0.5088e2 0.101e-1 0.11403e3; 
  0.5797e2 0.106e-1 0.11403e3; 
  0.6535e2 0.111e-1 0.11403e3; 
  0.7298e2 0.116e-1 0.11403e3; 
  0.8088e2 0.121e-1 0.11403e3; 
  0 0.1e-3 0.3801e3; 
  0.108e1 0.6e-3 0.3801e3; 
  0.262e1 0.11e-2 0.3801e3; 
  0.407e1 0.16e-2 0.3801e3; 
  0.628e1 0.21e-2 0.3801e3; 
  0.842e1 0.26e-2 0.3801e3; 
  0.1162e2 0.31e-2 0.3801e3; 
  0.1468e2 0.36e-2 0.3801e3; 
  0.1924e2 0.41e-2 0.3801e3; 
  0.236e2 0.46e-2 0.3801e3; 
  0.3012e2 0.51e-2 0.3801e3; 
  0.3633e2 0.56e-2 0.3801e3; 
  0.4546e2 0.61e-2 0.3801e3; 
  0.5391e2 0.66e-2 0.3801e3; 
  0.6586e2 0.71e-2 0.3801e3; 
  0.7643e2 0.76e-2 0.3801e3; 
  0.9069e2 0.81e-2 0.3801e3; 
  0.10275e3 0.86e-2 0.3801e3; 
  0.11835e3 0.91e-2 0.3801e3; 
  0.13103e3 0.96e-2 0.3801e3; 
  0.1469e3 0.101e-1 0.3801e3; 
  0.15942e3 0.106e-1 0.3801e3; 
  0 0.1e-3 0.11403e4; 0.338e1 
  0.6e-3 0.11403e4;
  0.662e1 0.11e-2 0.11403e4; 
  0.979e1 0.16e-2 0.11403e4; 
  0.1345e2 0.21e-2 0.11403e4; 
  0.1717e2 0.26e-2 0.11403e4; 
  0.2181e2 0.31e-2 0.11403e4; 
  0.2708e2 0.36e-2 0.11403e4; 
  0.3412e2 0.41e-2 0.11403e4; 
  0.4249e2 0.46e-2 0.11403e4; 
  0.538e2 0.51e-2 0.11403e4; 
  0.6689e2 0.56e-2 0.11403e4; 
  0.837e2 0.61e-2 0.11403e4; 
  0.10194e3 0.66e-2 0.11403e4; 
  0.12367e3 0.71e-2 0.11403e4; 
  0.14553e3 0.76e-2 0.11403e4;
  0.16962e3 0.81e-2 0.11403e4;
  0.19218e3 0.86e-2 0.11403e4;
  0.21412e3 0.91e-2 0.11403e4;
  0.23546e3 0.96e-2 0.11403e4;
  0.52e1 0.1e-3 3801;
  0.969e1 0.6e-3 3801;
  0.1415e2 0.11e-2 3801;
  0.1774e2 0.16e-2 3801;
  0.2128e2 0.21e-2 3801; 
  0.2446e2 0.26e-2 3801; 
  0.286e2 0.31e-2 3801; 
  0.3343e2 0.36e-2 3801; 
  0.4155e2 0.41e-2 3801; 
  0.5141e2 0.46e-2 3801; 
  0.6764e2 0.51e-2 3801; 
  0.8559e2 0.56e-2 3801; 
  0.11198e3 0.61e-2 3801; 
  0.1382e3 0.66e-2 3801; 
  0.17129e3 0.71e-2 3801; 
  0.20302e3 0.76e-2 3801; 
  0.23716e3 0.81e-2 3801; 
  0.26764e3 0.86e-2 3801; 
  0.29534e3 0.91e-2 3801; 
  0.3625e2 0.1e-3 0.95025e4;
  0.4035e2 0.6e-3 0.95025e4; 
  0.4444e2 0.11e-2 0.95025e4; 
  0.474e2 0.16e-2 0.95025e4; 
  0.509e2 0.21e-2 0.95025e4; 
  0.5446e2 0.26e-2 0.95025e4; 
  61 0.31e-2 0.95025e4; 
  0.6969e2 0.36e-2 0.95025e4; 
  0.855e2 0.41e-2 0.95025e4; 
  0.10498e3 0.46e-2 0.95025e4; 
  0.13504e3 0.51e-2 0.95025e4; 
  0.16578e3 0.56e-2 0.95025e4; 
  0.21008e3 0.61e-2 0.95025e4; 
  0.24829e3 0.66e-2 0.95025e4; 
  0.29631e3 0.71e-2 0.95025e4;];

And I have attached the objective function Fn, with the numbers expressed as rationals (rationals are easier to do symbolic work on).

These values should not be anything new, but I am more comfortable with having the expression posted so we are using the same thing for sure. Also because I do not have the symbolic toolbox I needed to convert from the symbolic package I am using to put into MATLAB form and might as well post that function; others might want to apply other optimizations to it.

[Delores Davis]

    All_data =

         1.1600e+000   600.0000e-006    38.0100e+000
         1.1100e+000     1.1000e-003    38.0100e+000
         1.1400e+000     1.6000e-003    38.0100e+000
       850.0000e-003     2.1000e-003    38.0100e+000
       620.0000e-003     2.6000e-003    38.0100e+000
       370.0000e-003     3.1000e-003    38.0100e+000
        10.0000e-003     3.6000e-003    38.0100e+000
         0.0000e+000     4.1000e-003    38.0100e+000
        90.0000e-003     4.6000e-003    38.0100e+000
       670.0000e-003     5.1000e-003    38.0100e+000
         1.2700e+000     5.6000e-003    38.0100e+000
         2.3100e+000     6.1000e-003    38.0100e+000
         3.2500e+000     6.6000e-003    38.0100e+000
         4.8400e+000     7.1000e-003    38.0100e+000
         6.3400e+000     7.6000e-003    38.0100e+000
         8.6100e+000     8.1000e-003    38.0100e+000
        10.9300e+000     8.6000e-003    38.0100e+000
        13.5400e+000     9.1000e-003    38.0100e+000
        16.2100e+000     9.6000e-003    38.0100e+000
        19.2500e+000    10.1000e-003    38.0100e+000
        22.4200e+000    10.6000e-003    38.0100e+000
        26.1800e+000    11.1000e-003    38.0100e+000
        29.5300e+000    11.6000e-003    38.0100e+000
        33.6200e+000    12.1000e-003    38.0100e+000
        50.0000e-003   600.0000e-006   114.0300e+000
         0.0000e+000     1.1000e-003   114.0300e+000
        60.0000e-003     1.6000e-003   114.0300e+000
       320.0000e-003     2.1000e-003   114.0300e+000
       680.0000e-003     2.6000e-003   114.0300e+000
         1.3200e+000     3.1000e-003   114.0300e+000
         2.0000e+000     3.6000e-003   114.0300e+000
         3.0100e+000     4.1000e-003   114.0300e+000
         4.0900e+000     4.6000e-003   114.0300e+000
         5.6700e+000     5.1000e-003   114.0300e+000
         7.4700e+000     5.6000e-003   114.0300e+000
        10.0200e+000     6.1000e-003   114.0300e+000
        12.9300e+000     6.6000e-003   114.0300e+000
        16.7700e+000     7.1000e-003   114.0300e+000
        21.0000e+000     7.6000e-003   114.0300e+000
        26.1700e+000     8.1000e-003   114.0300e+000
        31.5700e+000     8.6000e-003   114.0300e+000
        37.6500e+000     9.1000e-003   114.0300e+000
        44.0700e+000     9.6000e-003   114.0300e+000
        50.8800e+000    10.1000e-003   114.0300e+000
        57.9700e+000    10.6000e-003   114.0300e+000
        65.3500e+000    11.1000e-003   114.0300e+000
        72.9800e+000    11.6000e-003   114.0300e+000
        80.8800e+000    12.1000e-003   114.0300e+000
         1.0800e+000   600.0000e-006   380.1000e+000
         2.6200e+000     1.1000e-003   380.1000e+000
         4.0700e+000     1.6000e-003   380.1000e+000
         6.2800e+000     2.1000e-003   380.1000e+000
         8.4200e+000     2.6000e-003   380.1000e+000
        11.6200e+000     3.1000e-003   380.1000e+000
        14.6800e+000     3.6000e-003   380.1000e+000
        19.2400e+000     4.1000e-003   380.1000e+000
        23.6000e+000     4.6000e-003   380.1000e+000
        30.1200e+000     5.1000e-003   380.1000e+000
        36.3300e+000     5.6000e-003   380.1000e+000
        45.4600e+000     6.1000e-003   380.1000e+000
        53.9100e+000     6.6000e-003   380.1000e+000
        65.8600e+000     7.1000e-003   380.1000e+000
        76.4300e+000     7.6000e-003   380.1000e+000
        90.6900e+000     8.1000e-003   380.1000e+000
       102.7500e+000     8.6000e-003   380.1000e+000
       118.3500e+000     9.1000e-003   380.1000e+000
       131.0300e+000     9.6000e-003   380.1000e+000
       146.9000e+000    10.1000e-003   380.1000e+000
       159.4200e+000    10.6000e-003   380.1000e+000
         3.3800e+000   600.0000e-006     1.1403e+003
         6.6200e+000     1.1000e-003     1.1403e+003
         9.7900e+000     1.6000e-003     1.1403e+003
        13.4500e+000     2.1000e-003     1.1403e+003
        17.1700e+000     2.6000e-003     1.1403e+003
        21.8100e+000     3.1000e-003     1.1403e+003
        27.0800e+000     3.6000e-003     1.1403e+003
        34.1200e+000     4.1000e-003     1.1403e+003
        42.4900e+000     4.6000e-003     1.1403e+003
        53.8000e+000     5.1000e-003     1.1403e+003
        66.8900e+000     5.6000e-003     1.1403e+003
        83.7000e+000     6.1000e-003     1.1403e+003
       101.9400e+000     6.6000e-003     1.1403e+003
       123.6700e+000     7.1000e-003     1.1403e+003
       145.5300e+000     7.6000e-003     1.1403e+003
       169.6200e+000     8.1000e-003     1.1403e+003
       192.1800e+000     8.6000e-003     1.1403e+003
       214.1200e+000     9.1000e-003     1.1403e+003
       235.4600e+000     9.6000e-003     1.1403e+003
         9.6900e+000   600.0000e-006     3.8010e+003
        14.1500e+000     1.1000e-003     3.8010e+003
        17.7400e+000     1.6000e-003     3.8010e+003
        21.2800e+000     2.1000e-003     3.8010e+003
        24.4600e+000     2.6000e-003     3.8010e+003
        28.6000e+000     3.1000e-003     3.8010e+003
        33.4300e+000     3.6000e-003     3.8010e+003
        41.5500e+000     4.1000e-003     3.8010e+003
        51.4100e+000     4.6000e-003     3.8010e+003
        67.6400e+000     5.1000e-003     3.8010e+003
        85.5900e+000     5.6000e-003     3.8010e+003
       111.9800e+000     6.1000e-003     3.8010e+003
       138.2000e+000     6.6000e-003     3.8010e+003
       171.2900e+000     7.1000e-003     3.8010e+003
       203.0200e+000     7.6000e-003     3.8010e+003
       237.1600e+000     8.1000e-003     3.8010e+003
       267.6400e+000     8.6000e-003     3.8010e+003
       295.3400e+000     9.1000e-003     3.8010e+003
        40.3500e+000   600.0000e-006     9.5025e+003
        44.4400e+000     1.1000e-003     9.5025e+003
        47.4000e+000     1.6000e-003     9.5025e+003
        50.9000e+000     2.1000e-003     9.5025e+003
        54.4600e+000     2.6000e-003     9.5025e+003
        61.0000e+000     3.1000e-003     9.5025e+003
        69.6900e+000     3.6000e-003     9.5025e+003
        85.5000e+000     4.1000e-003     9.5025e+003
       104.9800e+000     4.6000e-003     9.5025e+003
       135.0400e+000     5.1000e-003     9.5025e+003
       165.7800e+000     5.6000e-003     9.5025e+003
       210.0800e+000     6.1000e-003     9.5025e+003
       248.2900e+000     6.6000e-003     9.5025e+003
       296.3100e+000     7.1000e-003     9.5025e+003

This is my dataset. I truncated some of the baseline on a couple of the curves while I was playing with the curve fitting toolbox. When I tried to enter your data I got an error.

Error using vertcat
Dimensions of matrices being concatenated are not consistent.

Then I got another error trying to run Fn.m

   function sqrerr  = Fn(Gain, TimeDelay)
   |
  Error: Function definitions are not permitted in this context.

I may be the worst operator of Matlab, ever!

This is Fn copied and pasted from my workspace using the data above.

@(Gain,TimeDelay)(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-1.9005e2)).*4.53e2-3.5025e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-1.9005e3)).*4.53e2-3.148e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-1.9005e1)).*4.53e2-4.4975e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-1.9005e2)).*4.53e2-4.294e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-1.9005e3)).*4.53e2-1.8536e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-5.7015e1)).*4.53e2-4.32e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-5.7015e1)).*4.53e2-4.51e2).^2+(exp(Gain.*(TimeDelay-1.01e-2).^2.*(-1.9005e1)).*4.53e2-4.3375e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-1.9005e3)).*4.53e2-4.244e2).^2+(exp(Gain.*(TimeDelay-1.01e-2).^2.*(-1.9005e2)).*4.53e2-3.061e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-1.9005e2)).*4.53e2-4.3832e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-1.9005e3)).*4.53e2-2.4998e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-1.9005e3)).*4.53e2-2.1584e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-1.9005e2)).*4.53e2-4.5192e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-1.9005e1)).*4.53e2-4.5184e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-1.9005e2)).*4.53e2-3.3465e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-4.75125e3)).*4.53e2-4.056e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-1.9005e3)).*4.53e2-1.5766e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-1.9005e3)).*4.53e2-4.1145e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-1.9005e3)).*4.53e2-4.3885e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-5.7015e2)).*4.53e2-2.1754e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-1.9005e1)).*4.53e2-4.5215e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-5.7015e2)).*4.53e2-3.992e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-1.9005e3)).*4.53e2-3.8536e2).^2+(exp(Gain.*(TimeDelay-1.06e-2).^2.*(-1.9005e2)).*4.53e2-2.9358e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-1.9005e2)).*4.53e2-4.2288e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-5.7015e2)).*4.53e2-4.2592e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-1.9005e3)).*4.53e2-4.3172e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-1.9005e2)).*4.53e2-4.3376e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-1.9005e2)).*4.53e2-4.4672e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-5.7015e2)).*4.53e2-3.693e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-1.9005e1)).*4.53e2-4.4816e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-5.7015e1)).*4.53e2-4.5295e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-1.9005e3)).*4.53e2-4.3526e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-1.9005e1)).*4.53e2-4.5186e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-5.7015e2)).*4.53e2-2.3888e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-5.7015e1)).*4.53e2-4.5232e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-1.9005e1)).*4.53e2-4.4666e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-4.75125e3)).*4.53e2-3.92e2).^2+(exp(Gain.*(TimeDelay-1.11e-2).^2.*(-5.7015e1)).*4.53e2-3.8765e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-5.7015e2)).*4.53e2-2.6082e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-4.75125e3)).*4.53e2-7.35e2./2.0).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-5.7015e1)).*4.53e2-4.1535e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-5.7015e2)).*4.53e2-4.3955e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-1.9005e3)).*4.53e2-4.2854e2).^2+(exp(Gain.*(TimeDelay-1.06e-2).^2.*(-1.9005e1)).*4.53e2-4.3058e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-1.9005e3)).*4.53e2-3.4102e2).^2+(exp(Gain.*(TimeDelay-1.21e-2).^2.*(-5.7015e1)).*4.53e2-3.7212e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-1.9005e2)).*4.53e2-4.4458e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-1.9005e1)).*4.53e2-4.5238e2).^2+(exp(Gain.*(TimeDelay-1.01e-2).^2.*(-5.7015e1)).*4.53e2-4.0212e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-5.7015e2)).*4.53e2-4.1888e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-4.75125e3)).*4.53e2-4.1265e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-5.7015e1)).*4.53e2-4.5168e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-5.7015e1)).*4.53e2-4.5268e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-4.75125e3)).*4.53e2-4.021e2).^2+(exp(Gain.*(TimeDelay-1.16e-2).^2.*(-5.7015e1)).*4.53e2-3.8002e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-1.9005e2)).*4.53e2-3.8714e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-5.7015e2)).*4.53e2-3.5106e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-1.9005e2)).*4.53e2-4.0754e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-4.75125e3)).*4.53e2-2.8722e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-5.7015e1)).*4.53e2-4.5294e2).^2+(exp(Gain.*(TimeDelay-1.21e-2).^2.*(-1.9005e1)).*4.53e2-4.1938e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-4.75125e3)).*4.53e2-2.4292e2).^2+(exp(Gain.*(TimeDelay-1.11e-2).^2.*(-1.9005e1)).*4.53e2-4.2682e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-5.7015e2)).*4.53e2-2.8338e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-1.9005e2)).*4.53e2-4.4138e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-1.9005e1)).*4.53e2-4.3946e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-1.9005e2)).*4.53e2-4.5038e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-1.9005e2)).*4.53e2-3.2197e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-4.75125e3)).*4.53e2-3.1796e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-5.7015e2)).*4.53e2-4.4962e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-4.75125e3)).*4.53e2-4.0856e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-4.75125e3)).*4.53e2-3.4802e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-1.9005e3)).*4.53e2-3.6741e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-1.9005e3)).*4.53e2-2.8171e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-5.7015e2)).*4.53e2-4.4638e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-5.7015e1)).*4.53e2-4.4298e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-4.75125e3)).*4.53e2-3.9854e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-1.9005e2)).*4.53e2-3.7657e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-4.75125e3)).*4.53e2-2.0471e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-5.7015e2)).*4.53e2-3.0747e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-4.75125e3)).*4.53e2-1.5669e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-1.9005e2)).*4.53e2-4.1667e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-1.9005e3)).*4.53e2-4.1957e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-1.9005e1)).*4.53e2-4.3679e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-1.9005e2)).*4.53e2-3.9909e2).^2+(exp(Gain.*(TimeDelay-1.16e-2).^2.*(-1.9005e1)).*4.53e2-4.2347e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-1.9005e3)).*4.53e2-4.0159e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-1.9005e2)).*4.53e2-4.4893e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-1.9005e2)).*4.53e2-3.6231e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-1.9005e1)).*4.53e2-4.5173e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-5.7015e2)).*4.53e2-3.8611e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-1.9005e1)).*4.53e2-4.5299e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-5.7015e1)).*4.53e2-4.0893e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-1.9005e3)).*4.53e2-4.4331e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-1.9005e1)).*4.53e2-4.4207e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-1.9005e1)).*4.53e2-4.5291e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-5.7015e2)).*4.53e2-3.2933e2).^2+(exp(Gain.*(TimeDelay-1.06e-2).^2.*(-5.7015e1)).*4.53e2-3.9503e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-5.7015e2)).*4.53e2-4.4321e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-5.7015e1)).*4.53e2-4.4553e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-5.7015e2)).*4.53e2-4.1051e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-1.9005e1)).*4.53e2-4.4439e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-5.7015e1)).*4.53e2-4.2143e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-1.9005e1)).*4.53e2-4.5069e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-1.9005e1)).*4.53e2-4.5189e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-1.9005e1)).*4.53e2-4.5233e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-1.9005e1)).*4.53e2-4.5263e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-4.75125e3)).*4.53e2-3.8331e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-5.7015e2)).*4.53e2-4.3583e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-5.7015e1)).*4.53e2-4.4007e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-5.7015e1)).*4.53e2-4.4891e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-5.7015e1)).*4.53e2-4.2683e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-5.7015e2)).*4.53e2-4.3119e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-5.7015e1)).*4.53e2-4.3623e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-5.7015e1)).*4.53e2-4.4733e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-5.7015e1)).*4.53e2-4.4999e2).^2

I don't understand why I have so much trouble copying and pasting in a Matlab friendly manner. You wanna work with my slightly truncated numbers or show me how to bypass the error issues??

[Walter Roberson]

Darn I had an editing error when I formatted that data for posting.

You will need to save that Fn definition into a file named Fn.m that is on your MATLAB path.

All_data = [1.11 0.0001 38.01;1.16 0.0006 38.01;1.11 0.0011 38.01;1.14 0.0016 38.01;0.85 0.0021 38.01;0.62 0.0026 38.01;0.37 0.0031 38.01;0.01 0.0036 38.01;0 0.0041 38.01;0.09 0.0046 38.01;0.67 0.0051 38.01;1.27 0.0056 38.01;2.31 0.0061 38.01;3.25 0.0066 38.01;4.84 0.0071 38.01;6.34 0.0076 38.01;8.61 0.0081 38.01;10.93 0.0086 38.01;13.54 0.0091 38.01;16.21 0.0096 38.01;19.25 0.0101 38.01;22.42 0.0106 38.01;26.18 0.0111 38.01;29.53 0.0116 38.01;33.62 0.0121 38.01;0.15 0.0001 114.03;0.05 0.0006 114.03;0 0.0011 114.03;0.06 0.0016 114.03;0.32 0.0021 114.03;0.68 0.0026 114.03;1.32 0.0031 114.03;2 0.0036 114.03;3.01 0.0041 114.03;4.09 0.0046 114.03;5.67 0.0051 114.03;7.47 0.0056 114.03;10.02 0.0061 114.03;12.93 0.0066 114.03;16.77 0.0071 114.03;21 0.0076 114.03;26.17 0.0081 114.03;31.57 0.0086 114.03;37.65 0.0091 114.03;44.07 0.0096 114.03;50.88 0.0101 114.03;57.97 0.0106 114.03;65.35 0.0111 114.03;72.98 0.0116 114.03;80.88 0.0121 114.03;0 0.0001 380.1;1.08 0.0006 380.1;2.62 0.0011 380.1;4.07 0.0016 380.1;6.28 0.0021 380.1;8.42 0.0026 380.1;11.62 0.0031 380.1;14.68 0.0036 380.1;19.24 0.0041 380.1;23.6 0.0046 380.1;30.12 0.0051 380.1;36.33 0.0056 380.1;45.46 0.0061 380.1;53.91 0.0066 380.1;65.86 0.0071 380.1;76.43 0.0076 380.1;90.69 0.0081 380.1;102.75 0.0086 380.1;118.35 0.0091 380.1;131.03 0.0096 380.1;146.9 0.0101 380.1;159.42 0.0106 380.1;0 0.0001 1140.3;3.38 0.0006 1140.3;6.62 0.0011 1140.3;9.79 0.0016 1140.3;13.45 0.0021 1140.3;17.17 0.0026 1140.3;21.81 0.0031 1140.3;27.08 0.0036 1140.3;34.12 0.0041 1140.3;42.49 0.0046 1140.3;53.8 0.0051 1140.3;66.89 0.0056 1140.3;83.7 0.0061 1140.3;101.94 0.0066 1140.3;123.67 0.0071 1140.3;145.53 0.0076 1140.3;169.62 0.0081 1140.3;192.18 0.0086 1140.3;214.12 0.0091 1140.3;235.46 0.0096 1140.3;5.2 0.0001 3801;9.69 0.0006 3801;14.15 0.0011 3801;17.74 0.0016 3801;21.28 0.0021 3801;24.46 0.0026 3801;28.6 0.0031 3801;33.43 0.0036 3801;41.55 0.0041 3801;51.41 0.0046 3801;67.64 0.0051 3801;85.59 0.0056 3801;111.98 0.0061 3801;138.2 0.0066 3801;171.29 0.0071 3801;203.02 0.0076 3801;237.16 0.0081 3801;267.64 0.0086 3801;295.34 0.0091 3801;36.25 0.0001 9502.5;40.35 0.0006 9502.5;44.44 0.0011 9502.5;47.4 0.0016 9502.5;50.9 0.0021 9502.5;54.46 0.0026 9502.5;61 0.0031 9502.5;69.69 0.0036 9502.5;85.5 0.0041 9502.5;104.98 0.0046 9502.5;135.04 0.0051 9502.5;165.78 0.0056 9502.5;210.08 0.0061 9502.5;248.29 0.0066 9502.5;296.31 0.0071 9502.5];

[Delores Davis]

>> Walter
Step #1, minimum value 129808 near (Gain = 12.3, TimeDelay = 0.0021)
Step #2, minimum value 129678 near (Gain = 12.414, TimeDelay = 0.002148)
Step #3, minimum value 129678 near (Gain = 12.4147, TimeDelay = 0.00214756)
Step #4, minimum value 129678 near (Gain = 12.4148, TimeDelay = 0.00214758)
Step #5, minimum value 129678 near (Gain = 12.4148, TimeDelay = 0.00214758)

So are we on the same page? And don't worry about editing mishaps, I've made several in the few minutes I've known you, it's reassuring. Did you notice I named the function after you?

[Walter Roberson]

It turns out that fminsearch does a good job, provided that you give it a starting point with a gain that is not too close to 0:

[min_g_td, minval] = fminsearch(@(x) Fn(x(1), x(2)), [1, 4/1000]);

This will give [12.4148129957773, 0.00214757370191429] as the location and 129677.636841265 as the best value. I have shown that this is not the best point possible, that there is a better one at [12.4148292375643, 0.00214757688127498] with value 129677.636840818, and basically you can still do some narrowing in around the 6th decimal place of the Gain and the 8th decimal place of the TimeDelay in order to get lower in around the 6th decimal place of minimum. There is still notable slope as you vary Gain around the 6th decimal place, but the TimeDelay is pretty flat near that 8th decimal place.

I had done an fminsearch() earlier but I had not trusted its result, thinking that it was caught in a local minima. However, I did a lot of stocastic testing (millions of random points) and found that Yes, the global minima is pretty very close to that location. fminsearch does, though, get caught in local minima near that point so you need to use other techniques if you need the Gain and TimeDelay more precisely.

I have included an improved implementation of Fn.m . This one can act on vectors or 2D arrays (grids).

[Delores Davis]

When reading previous papers regarding the gain and delay time, it wasn't reported past 3-4 decimal places, so we are all good there.

Now, the pesky business of cranking out confidence intervals. Previously, CI's were reported for the gain and time delay, so my boss expects the same.

[Walter Roberson]

For the time delay, would 3-4 decimal places be 0.00x and 0.000x seconds, or would it be xx.xxx and xx.xxxx ms ? For the gain, which is around 12.4, would the 12.4 be the 3 decimal places or would 12.415 be the 3 decimal places? Is it "decimal places" or "significant figures" being described?

I have been doing some testing on how sensitive the results are to measurement error.

The major sensitivity is to the measurement times, x: a change of 1/20 ms leads to a time delay change of 1/20 ms, so anything beyond the 1/10th millisecond position of the time delay is suspect. Which makes sense when you think about it: if the times are off systematically then the timedelay should be off by much the same amount. The times are given to the 1/10th millisecond so 1/20 millisecond measuring error would be standard for numeric analysis. A normally distributed time jitter with a standard deviation of 1/20th millisecond changes the time delay by about 1/40th millisecond so 3 significant figures should not really be reported unless you can justify that the timing precision is better than 1/10th millisecond.

random jitter near 0.005 in Y values or near 0.05 in the Phi values preserve the fundamental shapes of the curves, but might change the positions of the minima slightly. In particular, the third decimal place (fifth significant digit) of the gain can change by about 2.

[Delores Davis]

Significant figures. I haven't thought about such simplistic mathematical terms since elementary school. Essentially i need the results to look like this:

x.xx plus/minus x.x (p<0.5)

The paper I am looking at did not report the time delay, although it was used to calculate the gain, it isn't reported. (I attribute my initial belief that the time delay was constant to this lack of reporting.)

I am not terribly interested at this point in small errors in either the time delay or gain.

My gameplan is to manually clean up my traces by removing obvious noise manually, and if the recording is very noisy exclude the dataset completely. I also plan to normalize the data as far as aligning baselines to zero, and ensuring I have no negative numbers (which I don't think I did for the attached dataset).

[Walter Roberson]

I worked for a while investigating whether fminsearch() was definitely going to hit the right area for the global minima. I saw that it did not always do so for gains close to 0. I found the second best local minima and confirmed that fminsearch() escapes it even if started right in the best position, but that could be explained by the way it generates the initial simplex, and I think it is still plausible that for some data sets and some random starting positions it might find a second-best minima.

The solution to the known difficulty for the random starting point near gain 0, and for the potential difficulty with a lesser global minima, would be to take a variety of random starting points for the search; the search is fast enough to be able to afford that easily. Generate a couple of dozen starting points over the plausible ranges, fminsearch from each one, and you could be pretty sure that you had avoided any possible local minima.

But with it being essentially certain that the global minima had been found to within several decimal places, the question became how to define the confidence intervals.

I thought about that more and I realized that really one should do testing leaving out one of the 6 tests to see how much each individual test influenced the results. The larger the influence, the less confidence we could have in the results as the more important rejecting "bad" tests becomes.

I quickly generalized the process to test with Leave One Out: a robust result should change very little if you leave out any one of the samples.

I was astonished to see that leaving out the last sample of test 4, or leaving out the first sample of test 6, had quite an influence on the results.

If you run the tests repeatedly with Leave One Out, then you get a set of values, and you can calculate mean and standard deviation from those, and so reach a confidence interval on that basis.

For that particular set of data, the global minimum is near Gain = 12.4148536104427, TimeDelay = 0.00214758035802553, with mean Gain = 12.4196061450642 and 1 standard deviation 0.168669220360332; with 1.96 standard deviations giving 95% confidence, that is +/- 0.33059167190625. The minimum is 11.9063433093434 when leaving out the last sample of test 4, and the maximum is 13.8666325209502 when leaving out the first sample of test 6. If you leave out test 4 entirely, the Gain falls to 10.4556646070149 and if you leave out test 6 entirely, the Gain rises to 15.9517706896603.

It took me some further testing to realize that a lower Gain when a test is left out means that the test is pulling the Gain up, and that a higher Gain when a test is left out means that the test is pulling the Gain down.

I have included the latest version of the objective function, Fn2. The parameter order is:

1) Gain - scalar, or vector, or ndgrid or meshgrid
2) TimeDelay - scalar, or vector, or ndgrid or meshgrid
3) deltax - omit, or [], or scalar, or ndgrid or meshgrid
4) deltay - omit, or [], or scalar, or ndgrid or meshgrid
5) deltaphi - omit, or [], or scalar, or ndgrid or meshgrid

any of the above that are vectors or grids must be compatible with the other sizes.

6) selector - omit, or [], or numeric vector, or logical vector of length 6, or logical vector of length 126

deltax, deltay, deltaphi can be used as "jitter" on the measurements, to explore the effects of measurement error. For example,

Fn2(12.4147760266371, 0.00214757248582734, randn(10000,1)*0.00005)

could be used to explore jitter in the accuracy of the timing at the 50 microsecond standard deviation level.

A selector which is a logical vector of length 6 determines whether the particular test is included or not. A selector which is a logical vector of length 126 determines whether individual samples are included or not. A selector which is numeric (scalar or vector) gives the indices of samples to exclude, so you can easily do leave-one-out by passing the index of the sample to exclude.

   LOO = zeros(126, 2);
   for K = 1 : 126; LOO(K,:) = fminsearch(@(x) Fn(x(1), x(2), [], [], [], K, [12, 0.01]); end

[Delores Davis]

Interesting. Some of what you've written isn't quite making sense to me, right now, but I haven't played with Fn2 yet.

Giving a talk tomorrow, so I am completely swamped until that's done. I will talk to you after that is over, that is...if I don't go swimming in a vat of margaritas (and you are more than welcome to join me!)

[Delores Davis]

>> Fn2
Error using Fn2 (line 29)
Not enough input arguments.

I ran this as written:

Fn2(12.4147760266371, 0.00214757248582734, randn(10000,1)*0.00005)

And the output is a 10000 point answer that I assume is the selector. I don't know how to use a selector to choose which samples to exclude. But, since Fn2 worked in that syntax, maybe I made a mistake by trying to run Fn2 alone.

Also, this was an error.

   for K = 1 : 126; LOO(K,:) = fminsearch(@(x) Fn(x(1), x(2), [], [], [], K, [12, 0.01])); end
  Error using fminsearch (line 84)
  FMINSEARCH requires at least two input arguments or a structure with valid fields.

Is LOO the Leave out one command? I believe I am supposed to put something in the brackets, but I don't know what.

[Walter Roberson]

I missed a ')' when I posted it.

LOO = zeros(126, 2);
for K = 1 : 126; LOO(K,:) = fminsearch(@(x) Fn2(x(1), x(2), [], [], [], K), [12, 0.01]); end

LOO is the resulting array of values. LOO(K,1) is the best gain found when leaving out the K'th sample, and LOO(K,2) is the corresponding best time delay.

I programmed in several ways of indicating which samples are to be selected. The 6th argument can be any of the following:

• left out completely. In this case, all samples are used.
• a scalar or vector of positive integers such as 17 or 23:48 or [17, 23:29, 58]. When you use this form, the sample numbers you list are to be left out. The numbers are in the range 1 to 126
• a logical vector exactly 6 elements long, such as [true, false, true, true, false, true] . When you use this form, the tests corresponding to true are used and the tests corresponding to false are left out. With the example vector here, the 2nd and 5th test would be left out
• a logical vector exactly 126 elements long, such as (rand(126,1)<=0.98) . When you use this form, the samples corresponding to true are used and the samples corresponding to false are left out. With the example vector here, about 98% of the samples would be kept and about 2% would be ignored.

The line

Fn2(12.4147760266371, 0.00214757248582734, randn(10000,1)*0.00005)

has nothing to do with selecting samples. That code has to do with adding noise to the measurements. You would do that in order to determine how much the results depend upon measurement accuracy. In the above example, 10000 normally distributed random changes are generated with a standard deviation of 0.00005, and those changes will be added to the time measurements that are currently numbers such as 0.0016. The standard deviation of 0.00005 (50 microseconds) was chosen to correspond to the rounding range of the 0.0001 (100 microseconds) -- e.g., since anything between 0.00155 and 0.00165 would round to 0.0016, we need to examine the stability as we alter the timings.

The current implementation of Fn2 uses the same offset for all of the samples. At some point I should enhance the code to allow a different jitter for each sample. Not today.

Another jitter example:

jitterx = linspace(-2E-8,2e-8, 10001);
ptsjx = Fn2(12.4147760266371, 0.00214757248582734, jitterx);
[minval, minidx] = min(ptsjx);
jitterx(minidx)

The result, 1.796e-09, indicates that the combination Gain = 12.4147760266371, TimeDelay = 0.00214757248582734 would be an even better match if all of the measured times are 1.796E-09 early -- e.g., so the time 0.0016 actually represented 0.001600001796. std(ptsjx) is only about 1E-05 on the sum-of-squares, so over the 2E-08 time scale the results are very stable, which one would hope for. But you don't know until you check ;-)

[Delores Davis]

'Jitter' is giving me a headache. My instrument gives a datapoint every .5 ms soooo....I'm kinda not concerned about anything at the E-09 level for numbers that I will report to at one or two decimal places (but I may be thinking about things too simplistically). I might need some convincing that jitter is an issue to concern myself with.

The LOO equation is still not working correctly. I don't know why.

I believe that you may have misunderstood me when I talked about excluding data if it is too noisy. I meant that I would skip all six trials (the whole animal not just a singular run from an animal), so I may not need the leave one out analysis, unless, I am mis-comprehending its purpose.

>> LOO = zeros(126, 2);
for K = 1 : 126; LOO(K,:) = fminsearch(@(x) Fn(x(1), x(2), [], [], [], K), [12, 0.01]); end
Error using
symengine>@(Gain,TimeDelay)(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.16e-2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.06e-2).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.06e-2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0./6.25e2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0e-4).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0e-4).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-1.9005e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.01e-2).^2.*(-1.9005e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0e-4).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.01e-2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.11e-2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.21e-2).^2.*(-1.9005e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.16e-2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.6e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.06e-2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0./6.25e2).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.6e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0e-4).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-5.7015e2)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.6e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0e-4).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-8.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-9.1e-3).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.01e-2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.11e-2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.21e-2).^2.*(-5.7015e1)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.0e-4).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.6e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.6e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.6e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.0e-4).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-1.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-2.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-3.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-4.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-5.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-6.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2+(exp(Gain.*(TimeDelay-7.1e-3).^2.*(-4.75125e3)).*4.53e2-4.53e2).^2
Too many input arguments.

Error in @(x)Fn(x(1),x(2),[],[],[],K)

Error in fminsearch (line 190)
fv(:,1) = funfcn(x,varargin{:});

What are the terms: symengine, varargin, and funfcn? I also saw nargin in Fn2, and I don't know what that does either.

Playing around with Confidence intervals, and I did this, but the output makes no sense, based on our estimate for gain, as 14 isn't in the range.

SEM = std(gain)/sqrt(length(gain));           % Standard Error
ts = tinv([0.025  0.95],length(gain)-1);      % T-Score
CI = mean(gain) + ts*SEM;                     % Confidence Intervals

[17.674137609318844,21.946279785842890]

Wrote this last night and wow didn't click submit (surprised it was still here!)

Thank you for all the work you are putting into my question. I appreciate you so much.

[Walter Roberson]

You accidentally invoked Fn instead of Fn2. Looks like I had posted the wrong name at the time.

LOO = zeros(126,2);
for K = 1 : 126; LOO(K,:) = fminsearch(@(x) Fn2(x(1), x(2), [], [], [], K), [12, 0.01]); end

[Delores Davis]

Changing it to Fn2 worked. The output is a two column matrix of what I assume is gain and 0.0021 or 0.0022. Is that what I was supposed to get? What is the 2nd column?

[Walter Roberson]

My instrument gives a datapoint every .5 ms

But it doesn't. It gives you a datapoint approximately every .5 ms. There is going to be some variation in the time it gathers the data. Sometimes it will be earlier than the 0.5 ms, sometimes it will be later than 0.5 ms. For example one time it might be 0.49732 ms after the previous sample. But the calculations depend on it being 0.5000000000000000 ms, with even 0.4999999999999732 ms for one reading changing the result. 0.5 ms is one significant figure, so anywhere between 0.4950000000000000 ms and 0.5499999999999999 ms would report the same 0.5 ms. We need to know how important the precision of the timing is, so we need to study the effect of varying the timing. Which is what the jitter in the third parameter is for: it allows the timings to be advanced or delayed in unison (in the current version of Fn2).

This is a matter of establishing confidence boundaries. If you get results that depend upon the times being exactly such and such, then how much do the results vary if you sweep the alteration of the times through the range of possibilities? Your confidence in your results cannot exceed the variation you observe over that range. If the variation over the range is smaller than your reporting requirements then all is well -- but you won't know until you do the test. The code is there to allow the test to be run.

For example:

First establish by repeated fittings with fminsearch() that there is a computed minima at

bestp = [12.4148560859473 0.00214757857698551];

then

jitterx = linspace(-5E-5,5E-5, 10001);
r = zeros(length(jitterx),2); s = zeros(length(jitterx),1);
for K = 1 : length(jitterx); [r(K,:), s(K,:)] = fminsearch(@(x) Fn2(x(1), x(2), jitterx(K)), bestp); end
plot(jitterx,r(:,1))

and observe how the best gain that fminsearch finds varies quite irregularly. And then observe

plot(jitterx,r(:,2))

and see the very straight line: the fminsearch calculated best time delay is nearly completely linear changing very close to 1:1 with the TimeDelay. This makes some sense: if all of the times are delayed by deltaT then a TimeDelay change of deltaT should give the same result offset and lead to the same result. This clean result gives us some confidence in the search mechanism. You can use

c = polyfit(jitterx(:), r(:,2), 1)

to see that the linear change is 0.999997891696971. And std(r(:,2)) is almost exactly the same as std(jitterx) as would be expected for such a linear correlation.

But now return our attention to the notable variation of r(:,1). As we know that the differences between the adjusted times and the timedelay are staying nearly constant, we know that the adjusted surface should be the same except shifted in the timeward direction. And the timeward part is being accurately found, so the gain should stay the same, but obviously does not. What we must have, then, is found limitations on fminsearch: given a point very near to the minima, how well does fminsearch come in finding the minima? And the answer can be seen in the variation in r(:,1) values. std(r(:,1)) is on the order of 3E-5 so the 95% confidence level is on the order of 6E-5. Which is well within your reporting requirements, but now you have information on how well fminsearch does on finding minima from close by.

Likewise you have output readings, your y, that are measured to a certain number of decimal places. The true reading is +/- 0.005 of the reported reading. How much does that variation matter to the outcome? You can explore that through a jitter value in the 4th input parameter. The 5th input parameter is to explore the possibilities that the true Phi are different than the recorded Phi due to limitations on the number of decimal places.

[Delores Davis]

I am sucking at writing the code to add jitter to the other parameters. Need your help.

jitterx = linspace(-5E-5,5E-5, 10001);
r = zeros(length(jitterx),2); s = zeros(length(jitterx),1);
for K = 1 : length(jitterx); [r(K,:), s(K,:)] = fminsearch(@(x) Fn2(x(1), x(2), jitterx(K)), bestp); end
plot(jitterx,r(:,1))
plot(jitterx,r(:,2))
c = polyfit(jitterx(:), r(:,2), 1)

jittery = linspace(-5E-5,5E-5, 10001);
u = zeros(length(jittery),2); v = zeros(length(jittery),1);
for L = 1 : length(jittery); [u(L,:), v(L,:)] = fminsearch(@(y) Fn2(y(1), y(2), jittery(L)), bestp); end
plot(jittery,u(:,1))
plot(jittery,u(:,2))
d = polyfit(jittery(:), u(:,2), 1)

jitterz = linspace(-5E-5,5E-5, 10001);
p = zeros(length(jitterz),2); q = zeros(length(jitterz),1);
for M = 1 : length(jitterz); [p(M,:), q(M,:)] = fminsearch(@(z) Fn2(z(1), z(2), jitterz(M)), bestp); end
plot(jitterz,p(:,1))
plot(jitterz,p(:,2))
e = polyfit(jitterz(:), p(:,2), 1)

Essentially, I am learning how to make Matlab stay 'busy' and not learning much else.

[Delores Davis]

Step #1, minimum value 918155 near (Gain = 0, TimeDelay = 0)
Step #1: widening Gain search lower
Step #1: widening TimeDelay search lower
Step #2, minimum value 2.05049e+07 near (Gain = 0.003, TimeDelay = 0.0001)
Step #2: widening TimeDelay search higher
Step #3, minimum value 6.31486e+06 near (Gain = 0.00012, TimeDelay = 0.000106)
Step #3: widening TimeDelay search higher
Step #4, minimum value 498239 near (Gain = 4.8e-06, TimeDelay = 0.00010618)
Step #4: widening TimeDelay search higher
Step #5, minimum value 489409 near (Gain = 5.808e-06, TimeDelay = 0.000106185)
Step #5: widening TimeDelay search higher

Tried to test the reproducibility of this test on another animal, and it is obviously horrible, for un-obvious reasons. Do you have any idea why? I have attached the data for this other animal (wildtype so should fall into range with the data we've been working on).