Least squares linear regression with constraints

22 次查看(过去 30 天)
Hi everyone!
I am new to MATLAB and I have the following issue. I have the following data set:
x11 = [0.091, 0.068, 0.086, 0.091, 0.097, 0.064, 0.052, 0.066, 0.074];
x12 = [0.707, 0.920, 0.491, 0.527, 0.656, 0.474, 0.49, 0.267, 0.552];
y = [0.07, 0.029, 0.063, 0.077, 0.155, 0.043, 0.023, 0.057, 0.085];
where x11 and x12 are the measured inputs, and y is the known output.
I need to solve for the unknown coefficients a1, a2 and a3 in the following original equation:
y = a1*(x11^a2)*(x12^a3)
This is a non-linear problem of course, but I do linearize it by taking log on both sides, so that I have:
y2 = log10(y);
x1 = log10(x1);
x2 = log10(x2);
I then perform a linear regression analysis using the pinv function which makes use of a pseudoinverse matrix to find the coefficients.
My current code compiled based on the information I found on the Internet is following:
clear all close all
x11 = [0.091, 0.068, 0.086, 0.091, 0.097, 0.064, 0.052, 0.066, 0.074];
x12 = [0.707, 0.920, 0.491, 0.527, 0.656, 0.474, 0.49, 0.267, 0.552];
y = [0.07, 0.029, 0.063, 0.077, 0.155, 0.043, 0.023, 0.057, 0.085];
x1=log10(x11);
x2=log10(x12);
y2=log10(y);
n=length(x1);
a=[ones(n,1), x1', x2'];
c=pinv(a)*y2';
a1=c(1,:);
a2=c(2,:);
a3=c(3,:);
%here I take antilog to get back the original power relations:
y_pred = 10^(log10(a1)).*(x11.^a2).*(x12.^a3);
Although I do get a solution for all of the coefficients I am having two issues:
  1. I am not really satisfied with the output (poor R^2 coefficient)
  2. I want to make some constraints so that all of the coefficient stay positive, for example: using lsqlin didn't help much as the solution the is highly biased by my initial suggestions and the coefficient values either take the upper or lower bound value, which seems not to be correct.
I obtained a better solution using fminsearch approach without pre-linearizing the equation, but was curious if I can get comparable results taking the linear way.
I would be grateful for your help!

回答(2 个)

Torsten
Torsten 2023-6-30
编辑:Torsten 2023-6-30
You made a mistake in computing the correct coefficients (see below).
clear all close all
x11 = [0.091, 0.068, 0.086, 0.091, 0.097, 0.064, 0.052, 0.066, 0.074];
x12 = [0.707, 0.920, 0.491, 0.527, 0.656, 0.474, 0.49, 0.267, 0.552];
y = [0.07, 0.029, 0.063, 0.077, 0.155, 0.043, 0.023, 0.057, 0.085];
x1=log10(x11);
x2=log10(x12);
y2=log10(y);
n=length(x1);
a=[ones(n,1), x1', x2'];
c=a\y2';
a1_lin=10^c(1,:)
a1_lin = 34.9283
a2_lin=c(2,:)
a2_lin = 2.5909
a3_lin=c(3,:)
a3_lin = -0.4904
%here I take antilog to get back the original power relations:
y_pred = a1_lin.*(x11.^a2_lin).*(x12.^a3_lin);
norm(y_pred-y)
ans = 0.0695
fun = @(x)x(1)*x11.^x(2).*x12.^x(3)-y;
sol = lsqnonlin(fun,[a1_lin a2_lin a3_lin])
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
sol = 1×3
59.7121 2.7205 -0.1921
a1_nonlin = sol(1)
a1_nonlin = 59.7121
a2_nonlin = sol(2)
a2_nonlin = 2.7205
a3_nonlin = sol(3)
a3_nonlin = -0.1921
norm(fun(sol))
ans = 0.0666
  2 个评论
Alex Sha
Alex Sha 2023-8-4
移动:Bruno Luong 2023-8-4
If using direct nonlinear fitting
a1 59.737732722511
a2 2.72067588034148
a3 -0.192039313150924
While doing by linear transformation:
a1 34.9283448971434
a2 2.59089956584206
a3 -0.490446537714015
There are significant difference between two methods.
If want to all positive values:
a1 52.3606801051862
a2 2.62431938697266
a3 0

请先登录,再进行评论。


Bruno Luong
Bruno Luong 2023-6-30
编辑:Bruno Luong 2023-6-30
"This is a non-linear problem of course, but I do linearize it by taking log on both side"
That's the problem. When you transform problem by taking the log, if you have Gaussian noise on data, it will no longer be Gaussian, let alone normalized them with covariance matrix. So least-squares with transformed model is not correct way to sove the problem.
You need to use non-linear regression, fmincon or lsqnonlin (optimization toolbox required), that is a good way to go. fminsearch migh be not a good tool numerically, with poor convergence.
  2 个评论
Ace
Ace 2023-6-30
Thanks for your suggestion. But how do I put my contraint (positive values) not to force the model to show specific outcomes? Using upper and lower bounds seems too aggressive: the result might be actually equal to one of those. Thanks.
Bruno Luong
Bruno Luong 2023-6-30
If the constraints respect the physics then it is not agressive.
You have a range of ways to enter the contarints in fmincon and lsqnonlin, just read the doc.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by