Nonlinear Regression using a gaussian in a lattice

3 次查看(过去 30 天)
I have a system of equations where the relationship between input and output is derived from a pixel lattice:
where and stand for the distance between row and column of pixel i and j
In matrix form, this results in
where, for example, for a 3x3 lattice, A is:
I need to find α and σ for known x(0) and x(K) vectors that solve the linear system in K steps. This is, perform this linear regression:
K can be one (single time step - if possible or more than one timesteps)
I thought about applying some natural logarithms or use lsqn fucntion, but it is not straightforward. Any ideas?
Best
  2 个评论
Torsten
Torsten 2023-8-17
编辑:Torsten 2023-8-17
Shouldn't the regression be
x(K) = alpha^K * A^K * x(0)
?
Thus @Star Strider 's code should be changed to
syms a sigma
K = 3; % e.g.
dr = sym('dr',[3 3]);
dc = sym('dc',[3 3]);
A = exp((dr.^2+dc.^2)./(2*sigma^2));
A = a^K * A^K;
... etc

请先登录,再进行评论。

采纳的回答

Star Strider
Star Strider 2023-8-17
One approach (using smaller matrices for efficiency and to demonstrate and test the code) —
syms a sigma
dr = sym('dr',[3 3]);
dc = sym('dc',[3 3]);
A = exp((dr.^2+dc.^2)./(2*sigma^2));
A = A.*a
A = 
Afcn = matlabFunction(A, 'Vars',{a,sigma,[dr],[dc]}) % 'dr' = 'in3', 'dc' = 'in4'
Afcn = function_handle with value:
@(a,sigma,in3,in4)reshape([a.*exp((1.0./sigma.^2.*(in4(1).^2+in3(1).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(2).^2+in3(2).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(3).^2+in3(3).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(4).^2+in3(4).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(5).^2+in3(5).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(6).^2+in3(6).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(7).^2+in3(7).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(8).^2+in3(8).^2))./2.0),a.*exp((1.0./sigma.^2.*(in4(9).^2+in3(9).^2))./2.0)],[3,3])
dr = rand(3) % Provide The Actual Matrix
dr = 3×3
0.6108 0.6612 0.5732 0.4048 0.9912 0.1329 0.8820 0.9308 0.6280
dc = rand(3) % Provide The Actual Matrix
dc = 3×3
0.4187 0.9831 0.8027 0.0961 0.5272 0.7765 0.5988 0.7329 0.0857
x = randn(3,1) % Provide The Actual Vector
x = 3×1
0.2330 -0.4001 -0.6323
LHS = randn(3,1) % Provide The Actual Vector
LHS = 3×1
0.3999 0.0783 0.0984
B0 = rand(2,1);
B = lsqcurvefit(@(b,x)Afcn(b(1),b(2),dr,dc)*x, B0, x, LHS)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
B = 2×1
-0.0047 0.3826
fprintf('\ta = %8.4f\n\tsigma = %8.4f\n',B)
a = -0.0047 sigma = 0.3826
To create the (9x9) matrix, use:
dr = sym('dr',[9 9]);
dc = sym('dc',[9 9]);
and then the appropriate vectors.
I am not confident that I understand what you want to do, however this should get you started.
See if this works with your ‘x(k)’ and ‘x(k+1)’ (that I call ‘LHS’) vectors and produces the desired result.
.
  4 个评论
L
L 2023-8-22
Hi @Star Strider, is there anyway to also constraint alpha to be positive?
Thanks,
L
Star Strider
Star Strider 2023-8-22
Yes.
The lsqcurvefit function permits parameter range bounds. Since ‘alpha’ is the first parameter, the lower bounds would be set to [0 -Inf]. to constrain only ‘alpha’. You can set the upper bounds as well, however it is not necessary if you want to leave them unconstrained.
The revised lsqcurvefit call would then be:
B = lsqcurvefit(@(b,x)Afcn(b(1),b(2),dr,dc)*x, B0, x, LHS, [0 -Inf])
That should work.
.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 MATLAB 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by