参数拟合的最小二乘法

Question

我被要求使用最小二乘法来拟合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]

我的想法正确吗？

Answer 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)

Answer 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);

Answer 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')

关于非线性的注解：

由于使用了，这可以很好地处理线性模型。如果误差函数是非线性但凸的，如平方误差总和(SSE)，则这也会返回全局最优值。

请注意，一个非凸函数需要多个起始点来尝试捕获许多局部最优，然后取最好的一个仍然不能保证最优性。将约束添加到解决方案将需要惩罚函数或切换到受约束的求解器，因为fminsearch解决了无约束问题(除非您对其进行了适当的惩罚)。

易于修改的：

模型和误差函数易于修改。例如，如果您想要最小化绝对误差的总和，那么使用abs就很简单。

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