You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to fit 6 curves simultaneously to solve for 2 unknowns?
8 views (last 30 days)
Show older comments
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.
2 Comments
Delores Davis
on 7 Jun 2015
My bad. I copied the data to an excel file and then forgot to attach it. It's here now.
Answers (1)
Walter Roberson
on 7 Jun 2015
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?
46 Comments
Delores Davis
on 7 Jun 2015
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
on 7 Jun 2015
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
on 8 Jun 2015
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
on 9 Jun 2015
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
on 9 Jun 2015
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
on 9 Jun 2015
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
on 9 Jun 2015
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
on 9 Jun 2015
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'.
Delores Davis
on 10 Jun 2015
My main point was that if Igor can do it, Matlab can too. I just need to know how.
Walter Roberson
on 11 Jun 2015
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
on 11 Jun 2015
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
on 11 Jun 2015
Edited: Walter Roberson
on 12 Jun 2015
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
on 12 Jun 2015
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
on 12 Jun 2015
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
on 15 Jun 2015
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
on 17 Jun 2015
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
on 17 Jun 2015
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
on 17 Jun 2015
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
on 17 Jun 2015
You might have to code with .* instead of *.
Also, you should use "syms" instead of "sym"
syms Gain TimeDelay
Walter Roberson
on 17 Jun 2015
You might have to code with .* instead of *.
Also, you should use "syms" instead of "sym"
syms Gain TimeDelay
Delores Davis
on 18 Jun 2015
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
on 18 Jun 2015
Edited: Walter Roberson
on 18 Jun 2015
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
on 18 Jun 2015
Edited: Walter Roberson
on 19 Jun 2015
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
on 19 Jun 2015
You do have the symbolic toolbox; you are using "syms" which is provided by the symbolic toolbox.
Walter Roberson
on 19 Jun 2015
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
on 19 Jun 2015
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
on 19 Jun 2015
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
on 20 Jun 2015
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
on 20 Jun 2015
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
on 20 Jun 2015
Edited: Walter Roberson
on 20 Jun 2015
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
on 20 Jun 2015
>> 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
on 21 Jun 2015
Edited: Walter Roberson
on 21 Jun 2015
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
on 21 Jun 2015
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
on 21 Jun 2015
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
on 21 Jun 2015
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
on 22 Jun 2015
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
on 23 Jun 2015
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
on 24 Jun 2015
>> 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
on 24 Jun 2015
Edited: Walter Roberson
on 26 Jun 2015
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
on 26 Jun 2015
'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
on 26 Jun 2015
Edited: Walter Roberson
on 26 Jun 2015
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
on 27 Jun 2015
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
on 27 Jun 2015
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
on 29 Jun 2015
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
on 30 Jun 2015
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).
See Also
Categories
Find more on Audio Processing Algorithm Design in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)