曲线拟合和分布拟合

此示例说明如何执行曲线拟合和分布拟合，并讨论每种方法适用的情况。

在曲线拟合与分布拟合之间进行选择

曲线拟合和分布拟合是不同类型的数据分析。

• 当您要将某个响应变量建模为预测变量的函数时，请使用曲线拟合。

• 当您要为单一变量的概率分布建模时，请使用分布拟合。

曲线拟合

在以下试验数据中，预测变量为 period，即服用药物之后的时间。响应变量为 conc，即血液中的药物浓度。假设只有响应数据 conc 受试验误差的影响。

using TyPlot
using TyCurveFitting
using TyStatistics
period = [0.1;0.1;0.3;0.3;1.3;1.7;2.1;2.6;3.9;3.9;5.1;5.6;6.2;6.4;7.7;8.1;8.2;8.9;9.0;9.5;9.6;10.2;10.3;10.8;11.2;11.2;11.2;11.7;12.1;12.3;12.3;13.1;13.2;13.4;13.7;14.0;14.3;15.4;16.1;16.1;16.4;16.4;16.7;16.7;17.5;17.6;18.1;18.5;19.3;19.7];
conc = [0.01;0.08;0.13;0.16;0.55;0.90;1.11;1.62;1.79;1.59;1.83;1.68;2.09;2.17;2.66;2.08;2.26;1.65;1.70;2.39;2.08;2.02;1.65;1.96;1.91;1.30;1.62;1.57;1.32;1.56;1.36;1.05;1.29;1.32;1.20;1.10;0.88;0.63;0.69;0.69;0.49;0.53;0.42;0.48;0.41;0.27;0.36;0.33;0.17;0.20];

假设您要将血液浓度建模为时间的函数。绘制 conc 对 period 的图。

figure()
plot(period, conc, "o")
xlabel("Period");
ylabel("Blood Concentration");

假设 conc 作为 period 的函数，遵循双参数 Weibull 曲线。Weibull 曲线的形式和参数如下

其中，a 是水平缩放，b 是形状参数，c 是垂直缩放。

使用非线性最小二乘来拟合 Weibull 模型。

modelFun =  (x,p)->p[3] * (x/p[1])^(p[2]-1) * exp(-(x/p[1])^p[2]);
startingVals = [10, 2, 5];
nlModel = fit(modelFun, period, conc, "initial_param", startingVals);

在数据上绘制 Weibull 曲线。

hold("on")
xgrid = range(0,20,100);
line(xgrid, nlModel(xgrid), color="r");

拟合的 Weibull 模型存在问题。fitnlm 假设试验误差为加性误差，并且来自具有常量方差的对称分布。但散点图显示，误差方差与曲线高度成正比。而且，加性对称误差意味着可能出现负的血液浓度测量值。

更加现实的假设是，乘性误差在对数尺度上对称。依据该假设，通过取两侧的对数，对数据进行 Weibull 曲线拟合。使用非线性最小二乘来拟合曲线：

nlModel2 = fit((p, x) -> log(modelFun(p, x)), period, log.(conc), "initial_param", startingVals);
line(xgrid,exp.(nlModel2(xgrid)),color=[0 .5 0],linestyle="--");
legend(["Raw Data","Additive Errors Model","Multiplicative Errors Model"]);
hold("off")

模型对象 nlModel2 包含精度估计。最佳做法是检查模型的拟合优度。例如，在对数尺度上绘制残差图，以检验假设乘性误差具有常量方差是否正确。

分布拟合

假设您要对电子元件使用寿命的分布进行建模。变量 life 用于测量 50 个完全相同的电子元件的失效时间。

 life =[6.2;16.1;16.3;19.0;12.2;8.1;8.8;5.9;7.3;8.2;16.1;12.8;9.8;11.3;5.1;10.8;6.7;1.2;8.3;2.3;4.3;2.9;14.8;4.6;3.1;13.6;14.5;5.2;5.7;6.5;5.3;6.4;3.5;11.4;9.3;12.4;18.3;15.9;4.0;10.4;8.7;3.0;12.1;3.9;6.5;3.4;8.5;0.9;9.9;7.9];

用直方图可视化数据。

figure()
binWidth = 2;
lastVal = ceil(Int,Base.maximum(life));
binEdges = 0:binWidth:(lastVal+1);
_,_,h = histogram(life,binEdges);
xlabel("Time to Failure");
ylabel("Frequency");
ylim([0 10]);

由于使用寿命数据通常遵循 Weibull 分布，因此可以使用前一个曲线拟合示例中的 Weibull 曲线来拟合直方图。要尝试此方法，请将直方图转换为一个 (x,y) 点集，其中 x 是 bin 中心，y 是 bin 高度，然后对这些点进行曲线拟合。

counts = histcounts(life,binEdges);
binCtrs = binEdges[1:(end-1)] .+ binWidth/2;
h.set_facecolor([.9 .9 .9])
hold("on")
plot(binCtrs,counts,"o");
hold("off")

但是，对直方图进行曲线拟合容易出现问题，通常不建议这样做。

1.该过程违背了最小二乘拟合的基本假设。bin 计数是非负值，这意味着测量误差不对称。而且，bin 计数在分布的尾部和中心具有不同的变异性。最后，bin 计数的总和是固定的，这意味着它们不是独立的测量值。

2.如果您对条形高度进行 Weibull 曲线拟合，则必须对曲线进行约束，因为直方图是经过缩放的经验概率密度函数 (pdf)。

3.对于连续数据，对直方图（而不是数据）进行曲线拟合会丢弃信息。

4.直方图的条形高度依赖于 bin 边界和 bin 宽度的选择。

对于许多参数化分布，最大似然是更好的参数估计方法，因为它能避免这些问题。Weibull pdf 与 Weibull 曲线具有几乎相同的形式：

只不过 b/a 取代了尺度参数 c，因为该函数的积分必须为 1。要使用最大似然对数据进行 Weibull 分布拟合，请使用 wblfit 拟合，用 Weibull 生成分布。与最小二乘不同，最大似然会找到与缩放后的直方图最匹配的 Weibull pdf，而无需最小化 pdf 与条形高度之差的平方和。

res, = wblfit(life);
pd = Weibull(res[2],res[1]);

绘制经过缩放的数据直方图，并叠加绘制拟合的 pdf。

figure()
_,_,h = histogram(life,binEdges,normalization="pdf",facecolor=[.9 .9 .9]);
xlabel("Time to Failure");
ylabel("Probability Density");
ylim([0 0.1]);
xgrid = range(0,20,100);
pdfEst = pdf.(pd,xgrid);
hold("on")
line(xgrid, pdfEst, color=[0.00, 0.45, 0.74])
hold("off")