拟合参数的最小二乘法答案

【问题标题】：Least Squares Method to fit parameters拟合参数的最小二乘法
【发布时间】：2018-12-05 19:00:22
【问题描述】：

要求我使用最小二乘法来拟合y = α*exp(-β*x) 中的参数α 和β，

给出分数：

x = [1 2 3 4 5 6 7]
y = [9 6 4 2 4 6 9]

我无法确定我的矩阵应该是什么样子。我知道我应该取函数两边的自然对数以摆脱指数，并获得y值的自然对数，它们是：

ln_y = [2.19 1.79 1.39 0.69 1.39 1.79 2.19]

但是我的矩阵应该是什么样子，因为我剩下的是 ln(y) = ln(α) - β*x?

所以-β 列由1 组成，x 列将是我的x 值，但α 列应该包含什么？

这是我认为我应该得到的：

A = [1 1 1 1 1 1 1; 1 2 3 4 5 6 7]

我的想法正确吗？

【问题讨论】：

你为什么要尝试用这个函数来拟合你的实验数据？它看起来一点也不像指数函数。顺便说一句，你可以使用lsqnonlin。

标签： matlab least-squares non-linear-regression

【解决方案1】：

我们可以做的第一件事是在等式两边取自然对数ln（log in Matlab））：

y = α * e^(-β * x)

变成：

ln(y) = ln(α * e^(-β * x))
// Law of logarithms
ln(x * y) = ln(x) + ln(y) 

// thus:
ln(y) = ln(α) + ln(e^(-β * x))
Simplifying:
ln(y) = -β * x + ln(α)

现在我们将ln(y) 作为x 的线性函数，问题简化为找到最小二乘意义上的线性回归。让我们定义lny = log(y)和A = ln(α)，我们可以将问题重写为

lny = -β * x + A

在哪里

x = [1 2 3 4 5 6 7]
lny = [2.19 1.79 1.39 0.69 1.39 1.79 2.19]

对于 x 中的每个 x_i，我们可以如下计算 lny（重写为 x 的升幂）：

lny(x1) = A - β * x1
lny(x2) = A - β * x2
...
lny(xn) = A - β * xn

矩阵形式

LNY = X * [A β]'
Or,
X * [A β]' = LNY
// let Coefs = [A β]'
Coefs = X^-1 * LNY

在 Matlab 中

x = [1 2 3 4 5 6 7];
y = [9 6 4 2 4 6 9];
lny = log(y);
X = [ones(length(y), 1), -x']; % design matrix
coefs = X\lny'
% A = coefs(1) and β = coefs(2)
% ln(α) = A thus α = exp(A)
alpha = exp(coefs(1));
beta = coefs(2)

【讨论】：

【解决方案2】：

你几乎拥有它。第二行应该是-x。

x = [1 2 3 4 5 6 7]
y = [9 6 4 2 4 6 9]

logy = log(y)

n = length(x);
A = [ones(1,n); -x]

c = logy/A; %Solve for coefficients

alpha = exp(c(1))
beta = c(2);

【讨论】：

【解决方案3】：

在此示例中，导出最小二乘估计量是一个好主意。其他答案采用这种方法。

有一种快速而肮脏的方法，既灵活又方便。

只是数字而已。您可以使用fminsearch 来完成这项工作。

% MATLAB R2019a
x = [1 2 3 4 5 6 7];
y = [9 6 4 2 4 6 9];

% Create anonymous function (your supposed model to fit)
fh =@(params) params(1).*exp(-params(2).*x);

% Create anonymous function for Least Squares Error with single input
SSEh =@(params) sum((fh(params)-y).^2);     % Sum of Squared Error (SSE)

p0 = [1 0.5];                   % Initial guess for [alpha, beta]
[p, SSE] = fminsearch(SSEh,p0);
alpha = p(1);           % 5.7143
beta = p(2);            % 1.2366e-08  (AKA zero)

将结果绘制为健全性检查总是一个好主意（我经常搞砸，这一次又一次地节省了我）。

yhath=@(params,xval) params(1).*exp(-params(2).*xval);

Xrng = min(x)-1:.2:max(x)+1;
figure, hold on, box on
plot(Xrng,p(1).*exp(-p(2).*Xrng),'r--','DisplayName','Fit')
plot(x,y,'ks','DisplayName','Data')
legend('show')

关于非线性的说明：
由于convexity，这适用于线性模型。如果您的误差函数是非线性但凸的，如平方误差和 (SSE)，那么这也会返回全局最优值。

请注意，非凸函数将需要多个起点来尝试捕获许多局部最优值，然后采用最好的一个仍然不能保证最优性。为解决方案添加约束将需要惩罚函数或切换到约束求解器，因为fminsearch 解决了无约束问题（除非您正确惩罚它）。

易于修改：
很容易修改模型和误差函数。例如，如果您想最小化绝对误差的总和，则可以直接使用abs。

% Create anonymous function for Least Absolute Error with single input
SAEh =@(params) sum(abs(fh(params)-y));     % Sum of Absolute Error

【讨论】：