I need to find the value of tree variables: a
, b
and c
, by finding a global minimum for least squares method. My function is as follows:
f = (1/a)*(asinh((Z(i)/b)^(1/c))^(-1)
where i
is the index of vector Z
. The vector Z
has 9 values given in the task.
I also have a vector with 9 values for s
.
The least squares method need to sum differences between values computed by function f
with the values from vector s
. It should look something like this:
((s(i)-f(i))/s(i))^2
I also have the boundaries for value a
, b
and c
:
10e10>a>10e19
, 10e-7>b>50
, 10e-15>c>10
.
I tried to use the lsqnonlin
but I have no idea how to do it. I will appreciate any of your help!
I tried to do something like this:
function f=Fsigma(x, Z, sigma)
f=0;
for i=1:length(sigma)
f=f+((sigma(i)-((1/x(1,:))*(asinh((Z(i)/x(2,:))^(1/x(3,:)))^(-1))))/sigma(i))^2
end
end
and call this function in lsqnonlin
like this:
Z= [1.49E+18 1.49E+19 1.49E+20 1.99E+15 1.99E+16 1.99E+17 1.49E+13 1.49E+14 1.49E+15];
sigma = [55.1705 79.1016 105.636 25.4809 40.8572 61.7238 12.8147 21.4054 34.8319];
a=linspace(10e10,10e19);
b=linspace(10e-7,50);
c=linspace(10e-15,10);
x=[a; b; c];
p=lsqnonlin(Fsigma(x,Z,sigma));
OK @Agata, I'll go through it in detail, so you can learn some basics.
First, how to hand functions to functions: function handles aka. the @
-operator
fnc = @(x) Fsigma(x,Z,sigma);
fnc
is an object -- in fact, it is a something that points at the function Fsigma
. However, it even overshadows the additional inputs of Fsigma
by declaring that its only input is x
(this is called an anonymous function handle, indicated by the ()
, in between which you define the inputs that someone calling fnc
may use. The other inputs to Fsigma
in this line are the values of the variables that they have in this line.
boundaries Boundaries should be provided as vectors:
% bounds
lb = [ 10e10;
10e-7;
10e-15];
ub = [ 10e19;
50;
10];
initial guess + optimization call if you read the docs of lsqnonlin
, it requires an initial guess
% initial guess
x0 = ones(3,1);
% optimization call: x = lsqnonlin(fun,x0,lb,ub)
[p,fval] = lsqnonlin(fnc,x0,lb,ub)
Note to the cost function
lsqnonlin
is better if you do not provide a summed cost but an array of errors (see the docs, again). So I adjust you Fsigma
-function
function f = Fsigma(x, Z, sigma)
f = ((sigma-((1/x(1))./ asinh((Z./x(2)).^(1/x(3)))) )./sigma).^2;
end
If you want to have use a single output -- as most optimization algorithms require from their cost function -- you can use fmincon
Have a look at the full code
Z = [1.49E+18 1.49E+19 1.49E+20 1.99E+15 1.99E+16 1.99E+17 1.49E+13 1.49E+14 1.49E+15];
sigma = [55.1705 79.1016 105.636 25.4809 40.8572 61.7238 12.8147 21.4054 34.8319];
% bounds
lb = [ 10e10;
10e-7;
10e-15];
ub = [ 10e19;
50;
10];
% initial guess
x0 = ones(3,1);
% create an anyonymous function handle (using @(x)
fnc = @(x) Fsigma(x,Z,sigma);
% optimization call: x = lsqnonlin(fun,x0,lb,ub)
[p,fval] = lsqnonlin(fnc,x0,lb,ub);
% optimization call: x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)
fnc2 = @(x) sqrt(sum(fnc(x).^2));
[p2,fval2] = fmincon(fnc2,x0,[],[],[],[],lb,ub);
% cost function
function f = Fsigma(x, Z, sigma)
f = ((sigma-((1/x(1))./ asinh((Z./x(2)).^(1/x(3)))) )./sigma).^2;
end
PS: note that this community is not a coding service, so do your readings the next time!