MATLAB单调三次多项式拟合

Question

使用MATLAB编写单调三次多项式拟合函数。
x为实数，范围是[0,5]。
y为实数，范围是[0,5]。
编写x到y的单调三次多项式拟合函数。

MATLAB中有多项式拟合函数polyfit，可以拟合三次多项式，但是不能拟合『单调』三次多项式。
基于均方误差的拟合。

Answer 1

『单调』这个要求比较奇怪。
如果保证了单调性，原则上是不能保证比较好的拟合效果的。
最好把数据贴出来再看看。追问

要求是单调三次多项式。
要求代码。
例如x=rand(20,1); y=rand(20,1);
求它们之间的，单调三次多项式的fit代码。

追答

1、明确一下：你要的“单调”是指在给定数据x区间内单调，还是在整个实数域内单调？这二者是不一样的，思路也不同。

2、在保证单调性的前提下，拟合的效果可能不会太理想，只能是保证总误差最小（最小二乘）。

追问

1）在给定区间单调。区间范围是[0,5]。
2）保证总误差最小即可。

追答

基本思路：

要求多项式单调，只会有单调增和单调减两种情况，分别按照单调减和单调增两种限制进行拟合，然后再从二者中取方差较小的那个。

 

程序有点长，不太好贴。先上几张效果图（数据随机生成）：

 

花了不少时间研究这个问题，楼主能告诉我之所以提出单调这个要求的背景吗？



追问

提出这个『单调』限制的背景，主要跟应用有关。因为我现在做的一个事情，是提出一个模型预测一个prediction_score。然后和true_score的真实值进行比较。在比较前，要做一个『单调』的拟合，这样可以使得拟合后的fit_score和true_score更加接近。『单调』可以保证拟合的时候，不会改变数据的单调性，就是预测的顺序。

追答

原来如此。

由于字数限制，这里只能贴出部分代码，完整的代码以附件的形式上传，供参考。

function zd
x=rand(20,1)*5;
y=rand(20,1)*5;
f = @(c,x) c(1)*x.^3 + c(2)*x.^2 + c(3)*x +c(4);
% 基本思想：要求多项式单调，只会有单调增和单调减两种情况，分别按照单调减和单调
% 增两种限制进行拟合，然后再从二者中取方差较小的那个
 
% 目标函数为方差最小
se = @(c) sum( (y - f(c,x)).^2 );
c0 = [1 1 1 1];
c1 = fmincon(se,c0,[],[],[],[],[],[],@nlcon1);
c2 = fmincon(se,c0,[],[],[],[],[],[],@nlcon2);
% 不考虑单调性，直接拟合（也可以用polyfit）
c3 = lsqcurvefit(f,c0,x,y);
% 绘图对比拟合效果
clf
xs = linspace(0,5,100);
h= plot(x,y,'.',xs,f(c3,xs),'r--',xs,f(c1,xs),'c:',xs,f(c2,xs),'m:');
% 比较单调减和单调增两种拟合函数，取方差较小的（加粗显示）
if se(c1) < se(c2)
    c = c1;
    set(h(3),'LineStyle','-','linewidth',3)
else
    c = c2;
    set(h(4),'LineStyle','-','linewidth',3)
end
legend('原始数据','无单调性限制','单调减','单调增',0)
% 显示拟合多项式
disp(vpa(poly2sym(c,'x'),4))

function [c, ce] = nlcon1(c0)
% 三次多项式单调减的约束条件：导数g(x)=f'(x)<=0，即
% (1) 在区间的端点g(0)<=0, g(5)<=0
% (2) 如果对称轴x=x0在0-5区间内，则g(x0)<=0
g = @(x) 3*c0(1)*x^2 + 2*c0(2)*x + c0(3);
ce = [];
c = [g(0); g(5); -1];
x0 = -c0(2)/(3*c0(1));
if x0 > 0 && x0 < 5
    c(3) = g(x0);
end