How to fit 6 curves simultaneously to solve for 2 unknowns?
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
Walter Roberson
on 7 Jun 2015
The dataset didn't make it to attachment :(
Delores Davis
on 7 Jun 2015
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
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
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
Delores Davis
on 9 Jun 2015
Walter Roberson
on 10 Jun 2015
Sorry, I am not familiar with IgorPro at all.
Delores Davis
on 10 Jun 2015
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
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
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
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
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
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
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
Delores Davis
on 19 Jun 2015
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
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 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
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
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
Delores Davis
on 24 Jun 2015
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
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
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
Delores Davis
on 30 Jun 2015
Categories
Find more on Audio Processing Algorithm Design in Help Center and File Exchange
on 7 Jun 2015
on 30 Jun 2015
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
