How can I non linear fit this equation?

Dear all,

I know the equation which is used for fitting already,something like:

 y=E-109736.623./[x1.*(x-p1).^2+x2.*(x-p2).^2+x3.*(x-p3).^2]

I have x value, y value, x1, x2, x3 value, respectively, they are corresponding to fit.txt file column1, column2, column3, column4, column5, respectively.

How can I get E, p1, p2, p3 by using this equation?

Thanks

Answer 1

[Walter Roberson]

The easiest way is the Curvefitting Toolbox

[Pengju]

I tried but failed, as x1, x2, x3 can not be input to the data...

[Walter Roberson]

Yes, they can be. See use of anonymous functions. http://www.mathworks.com/help/curvefit/fit.html#bto2vuv-11

[Pengju]

I will try it again

[Pengju]

The code does not work...

data=xlsread('error_fit1.xlsx');
x=data(:,1);
y=data(:,2);
x1=data(:,3);
x2=data(:,4);
x3=data(:,5);
c=109736.623;
result = @(E, p1, p2, p3, x, y)(E-c./[x1.*(x-p1).^2+x2.*(x-p2).^2+x3.*(x-p3).^2]);
Algometryresult = fit(x, y, x1,x2,x3, result, ...
    'StartPoint', [33650, 10, 1, 0.8], ...
    'Lower', [-Inf, -Inf, -INf, -Inf], ...
    'Robust', 'LAR' )
plot( Algometryresult, x, y, x1,x2,x3);

[Walter Roberson]

Your anonymous function should exclude y, as y is your dependent. fit() knows that the calculated value of the objective needs to match the dependent variable, here passed in as the second parameter to fit.

Although that might affect the speed of calculation I am not sure it would affect the fit.

You should be using () instead of [] -- but it would not make any difference to this calculation.

I do not know at the moment why it might not have worked.

[Walter Roberson]

There is an excellent fit near

 [33690.9417099533603 .. 33690.9417107163827, 3.14040272505204721 .. 3.14040351360411218, 2.63427756239212973 .. 2.63427798866314777, 1.34093810668077085 .. 1.34093907346367525]

In the time I allowed my searching to run, the best location was

[33690.941710391533, 3.14040320712546128, 2.63427777512871719, 1.34093862330272717]

Code:

data = xlsread('error_fit1.xlsx');
x=data(:,1);
y=data(:,2);
x1=data(:,3);
x2=data(:,4);
x3=data(:,5);
c=109736.623;
NP = size(x,1);

Normk = @(Z,k) sqrt(sum(Z.^2,k));

repparm = @(var) repmat(permute(var, [ndims(var)+1, 1:ndims(var)]), [NP, 1]);
repinp = @(var1,var2) repmat(var1, [1, size(var2)]);
obj = @(E,p1,p2,p3) permute(Normk(repparm(E) - c./(repinp(x1,E).*(repinp(x,E) - repparm(p1)).^2 + repinp(x2,E) .* (repinp(x,E)-repparm(p2)).^2 + repinp(x3,E) .* (repinp(x,E) - repparm(p3)).^2) - repinp(y,E),1), [2:ndims(E)+1,1]);

obj is now a fully-vectorized objective function that can be applied to 4 dimensional grids of E, p1, p2, p3 values, giving one error residue value per grid location. You could wrap this in another function

obj1 = @(Ep1p2p3) obj(Ep1p2p3(1), Ep1p2p3(2), Ep1p2p3(3), Ep1p2p3(4))

and now obj1 can be used in a minimizer such as fminunc

[Pengju]

Very good....

how can run this code? I mean shall I set some initial value?

[Walter Roberson]

For example,

search_E = linspace(30000,40000, 50);
search_p1 = linspace(-10, 10, 50);
search_p2 = linspace(-10, 10, 50);
search_p3 = linspace(-10, 10, 50);
[gE, gp1, gp2, gP3] = ndgrid(search_E, search_p1, search_p2, search_p3);

residues = obj(gE, gp1, gp2, gp3);

Now residues would be a 50 x 50 x 50 x 50 array of results, sampled at each possible combination of search_E, search_p1, search_p2, and search_p3.

This will not in itself tell you exactly where the minima is, but by examining where the lower values are, you can find the lower lying areas that you can go back to with a finer mesh. Once you are reasonably sure that you are somewhere in the bowl of the global minima, you can use a point inside the bowl as the starting point for something like fmincon using obj1 as the objective function.

[Pengju]

??? Maximum variable size allowed by the program is exceeded.

Error in ==> repmat at 104
    B = A(subs{:});

Error in ==> @(var)repmat(permute(var,[ndims(var)+1,1:ndims(var)]),[NP,1])

Error in ==>
@(E,p1,p2,p3)permute(Normk(repparm(E)-c./(repinp(x1,E).*(repinp(x,E)-repparm(p1)).^2+repinp(x2,E).*(repinp(x,E)-repparm(p2)).^2+repinp(x3,E).*(repinp(x,E)-repparm(p3)).^2)-repinp(y,E),1),[2:ndims(E)+1,1])

Error in ==> error_fit at 24
residues = obj(gE, gp1, gp2, gp3);

[Pengju]

I guess I need decrease the array size, and I did it to20X20X20X20, but my memory size is not enough...

[Walter Roberson]

On my system, 18 x 18 x 18 x 18 came in less than 500 megabytes of memory used.

[Pengju]

I tried 19X19X19X19, it works....

[Pengju]

Dear Walter, I am sorry, I still can not figure out the minima.

Best

[Walter Roberson]

minidx = min(residues(:));
min_is_near = [gE(minidx), gp1(minidx), gp2(minidx), gp3(minidx)];

The larger the matrix you use originally, the more likely it is that you will find a good minima.

The code that I posted does make copies of these arrays, so you end up with about 6 arrays that are 84 by (19 x 19 x 19 x 19). The copies themselves add up, so when you do not have a lot of memory, you might want to break the grid up into sections to search and find the best of the results. For example, if you used 20 x 20 x 20 x 20 then you could search first on (1:5, :, :, :, :) then on (6:10, :, :, :, :) and so on.