1

我想为绘制的数据拟合一个衰减指数。我没有曲线拟合或优化工具箱。

x = [0    0.0036    0.0071    0.0107    0.0143    0.0178    0.0214    0.0250    0.0285    0.0321    0.0357    0.0392    0.0428    0.0464    0.0464];
y = [1.3985    1.3310    1.2741    1.2175    1.1694    1.1213    1.0804    1.0395    1.0043    0.9691    0.9385    0.9080    0.8809    0.7856    0.7856];

figure()
plot(x,y,'*')

我如何在 MATLAB 中实现这一点？

4

2 回答 2

8

假设您在输入和输出点之间存在高斯分布误差，并且假设这些误差是相加的，您可以通过经典的最小二乘法来解决这个问题。它归结为拥有一个超定线性方程组，其中每个约束定义一个输入-输出观察。以最少的残差解决这个超定线性系统是您正在寻找的解决方案。

警告

Jubobs 在下面对我的评论中提出了一个非常有趣的观点。最小化此残差的参数通常不会最小化原始问题的残差。这个线性化步骤允许我们以更简单的方式求解参数，但这不是等效的问题。但是，它通常在实践中被接受，因为解决方案足够好。

为了把它变成一个线性系统，我们需要做一些巧妙的重新排列。因为要使用指数模型拟合一系列点，所以输入x和输出之间的关系y是：

为了使其成为“线性”，我们可以取两边的自然对数：

通过使用 的事实ln(ab) = ln(a) + ln(b)，我们有：

也知道，这简化为：

如您所见，上述等式现在相对于对数空间是“线性的”。给定一堆x和y值，所以(x_1, x_2, ..., x_n)和(y_1, y_2, ..., y_n)，我们可以在一个线性系统中将一堆方程连接在一起：

如果我们让ln(A) = A'符号方便，并重新排列它以使其为矩阵形式，我们得到：

因此，我们只需要求解A'和b，您可以通过伪逆来解决。具体来说，上述问题的形式为：

因此，我们需要求解X，所以：

M^{+}是矩阵的伪逆。完成后，只需带exp操作员A'获取原始A. MATLAB 具有非常高效的线性系统求解器和最小二乘求解器。具体来说，您可以使用\orldivide运算符。您所要做的就是根据值创建M矩阵x，创建值向量y并求解您的系统。真的很简单：

x = ...; %// Define numbers here - either row or column vectors
y = ...;
M = [ones(numel(x),1), x(:)]; %// Ensure x is a column vector
lny = log(y(:)); %// Ensure y is a column vector and take ln

X = M \ lny; %// Solve for parameters
A = exp(X(1)); %// Solve for A
b = X(2); %// Get b

因此，使用您的x和y价值观，这就是我得到的：

A =

    1.3882

b = 

   -11.508

如果我们绘制上述点以及适合该线的指数曲线，我们可以这样做：

xval = linspace(min(x), max(x));
yval = A*exp(b*xval);
plot(x,y,'r.',xval,yval,'b');

第一行代码定义了一组介于我们数据集x的最小值和最大值之间的x值。对于下一行，我们获取这些x值并通过我们的指数模型运行它们。最后，我们将原始数据点以及通过上述过程找到的参数与指数曲线一起绘制。点为红色，而线为蓝色。

我们得到：

在此处输入图像描述

我觉得这看起来很不错！对于那些注意到的人来说，上面的图看起来与 MATLAB 生成的普通图和图形窗口有点不同。该图是在 Octave 中生成的，因为我目前正在使用的计算机上没有 MATLAB。但是，上面的代码应该仍然可以在 MATLAB 中运行。

于 2015-03-27T16:37:04.290 回答
1

对于像我这样更愚蠢的人，整个代码（感谢 rayryeng）是：

x = [0    0.0036    0.0071    0.0107    0.0143    0.0178    0.0214    0.0250    0.0285    0.0321    0.0357    0.0392    0.0428    0.0464    0.0464];
y = [1.3985    1.3310    1.2741    1.2175    1.1694    1.1213    1.0804    1.0395    1.0043    0.9691    0.9385    0.9080    0.8809    0.7856    0.7856];

M = [ones(numel(x),1), x(:)]; %// Ensure x is a column vector
lny = log(y(:)); %// Ensure y is a column vector and take ln

X = M\lny; %// Solve for parameters
A = exp(X(1)); %// Solve for A
b = X(2); %// Get b

xval = linspace(min(x), max(x));
yval = A*exp(b*xval);
plot(x,y,'r.',xval,yval,'b');
于 2015-03-27T16:50:22.840 回答