请教高手曲线拟合问题(椭圆拟合)

babecheer

想把图中这些数据点拟合为一条如红色所示的曲线，不知道该用什么函数来拟合，请大虾帮忙啊！

具体数据如下：
x=[0 0.099 0.15 0.1975 0.25 0.3265 0.411 0.5619 0.7265 0.951 1.126 1];y=[0.9 1.1 0.9 1.1 0.9 1.1... %绘图数值
0.975 0.9859 0.877 0.724 0.4152 0];...

bainhome

直接用cftool试试看，不过在你的图形中所表示的曲线估计是无法用单一函数拟合的，后半部分甚至根本不是一般意义上的“函数”

linji526

楼主所画的哪种曲线应该不可能是一个函数，因为函数不可能“一对多”，但图中在x=1附近一个x对应两个y了
下面是用1stopt的快速公式拟合搜索得到的一个结果：
Function: y = Exp(p1+p2*x^p3)
Algorithms: 准牛顿法(BFGS) + 通用全局优化法
Root of Mean Square Error (RMSE): 0.165323408211398
Sum of Square Error (SSE): 0.327981951631592
Correlation Coef. (R): 0.846385252269932
R-Square: 0.716367995260037
Determination Coef. (DC): 0.715567379672821
Parameters Name Parameter Value
=============== ===============
p1       0.00721837247197112
p2       -0.879713674792157
p3       4.67925360238175
======== 输出结果 =========
No. Observed Y Calculated Y
1 0.9 1.00724448772126
2 1.1 1.00722679499353
3 0.9 1.00712084412137
4 1.1 1.00679660997584
5 0.9 1.0058955438356
6 1.1 1.00254773319365
7 0.975 0.993517730119734
8 0.9859 0.949267218769211
9 0.877 0.826928998610224
10 0.724 0.502482705726298
11 0.4152 0.217481430851873
12 0 0.417907441698515

linji526

这种问题的困难之处我想就在于该选取什么样的函数进行拟合，通常最好是建立在理论基础之上，比如人口模型就是这样，如果单纯从这些点决定用什么进行拟合往往结果不是很好。

babecheer

恩，谢谢大家的热心回复，我主要是想用这些点来确定一条包络线，觉得有点像椭圆的形状，但是拟合出来觉得不好，最主要是想让两个端点都经过1，中间的曲线大致呈向外凸的形状，不过觉得难以找到合适的方法，还请大家多出良策^_^。不知能否用椭圆来直接拟合呢？

xianggui7895

我对这个问题也不是很理解，用最小二乘法进行椭圆的拟合不知道行不行？

babecheer

各位大虾，小弟再请教一个基础的问题，如果想用椭圆来拟合的话，能否直接用椭圆的标准方程形式x^2/a^2+y^2/b^2=1（隐函数）来拟合，我是把它转换成y=f(x)来拟合的，但是正如楼上所说的，这条曲线不是一个标准意义上的函数，所以用y=f(x)形式拟合出来效果很不好，请问各位用隐函数的方程该怎么拟合？

bainhome

椭圆拟合目前比较可靠流行的方法是6楼所说的最小二乘法，不过比较没技术含量：
椭圆标准函数为A x^2 + B x y + C y^2 + D x + E y + F = 0
并且 B^2 - 4 A C < 0 （否则是双曲线或者抛物线）。
建立f(A,B,C,D,E,F)=∑(A x^2 + B x y + C y^2 + D x + E y + F)
强制B^2 - 4 A C=-1，使之成为等式约束条件。
利用拉格朗日乘子法建立函数：
F(A,B,C,D,E,F)=∑(A x^2 + B x y + C y^2 + D x + E y + F)+λ(B^2 - 4 A C+1)
分别对A～F六个自变量求导，极值点处数值为0，加上刚才的等式约束，正好建立恰定非线性方程组，用fsolve求解。A~F和λ七个未知数的数值。
据说还有人为此方法申请了专利，那可真是TNND需要B4一下...

babecheer

非常感谢，再请教下bainhome：建立f(A,B,C,D,E,F)=∑(A x^2 + B x y + C y^2 + D x + E y + F)，是否是对所有离散点进行求和呢？此公式是否有一点的依据？谢谢︿＿︿

cumt123

看这种拟合意义是什么，再确定方法，否则就少量应用价值

bainhome

是否是对所有离散点进行求和呢

当然

此公式是否有一点的依据

最小二乘法和拉格朗日乘子法都是最基本的方法，"依据"的事情感觉就不用谈了吧，随便一google，很多的哦。

babecheer

再次感谢。

babecheer

再请教下bainhome：建立f(A,B,C,D,E,F)=∑(A x^2 + B x y + C y^2 + D x + E y + F)，个人对这里有点疑惑，比如说在上述公式的基础上最终求出Ａ～Ｆ，但是如果ｘ１和ｙ１这组数据下A x^2 + B x y + C y^2 + D x + E y + F＝１，ｘ２和ｙ２这组数据下A x^2 + B x y + C y^2 + D x + E y + F＝－１的时候，f(A,B,C,D,E,F)=∑(A x^2 + B x y + C y^2 + D x + E y + F)＝０依然满足，但显然这时建立的f(A,B,C,D,E,F)=∑(A x^2 + B x y + C y^2 + D x + E y + F)这个方程本身就存在一定的问题了，个人愚见，还请bainhome继续指教。

bainhome

指教不敢当，在这里千万不要这样讲，志同道合者一起讨论问题而已。
今天简单试了一下MATLAB的几个函数，隐函数的拟合似乎比较困难，我提的方法也尝试失败，正想尝试转成优化问题看看，但时间太晚了，作罢。最后用1stopt简单做了一下，在给定初值范围之后，效果不大理想，主要是没有经过给定点，难度与预料相同：在于最后两个点。不过像这种需要经过定点的拟合问题shamohu好像做过，我已经PM他，请他出手做一下。

babecheer

恩，希望大家能多出良策，最终解决这个问题，呵呵。我再找找相关方面的资料，希望能有所突破。

linji526

我提个比较直接的方法、
观察提供的数据，椭圆的中心大致在(0，0.4)，短轴长等于0.6，长轴长取1.35。这样可以得到一个椭圆方程为x^2/1.35^2+(y-0.4)^2/0.6^2=1
以下是Matlab命令：
>> x=[0 0.099 0.15 0.1975 0.25 0.3265 0.411 0.5619 0.7265 0.951 1.126 1];
>> y=[0.9 1.1 0.9 1.1 0.9 1.1 0.975 0.9859 0.877 0.724 0.4152 0];
>> plot(x,y,'r*')
>> hold on
>> ezplot('x^2/1.35^2+(y-0.4)^2/0.6^2-1',[0,1.5])
>>
结果为

bainhome

已解决。
主函数

1. function ff=MainElipse(x)
   
2. ff=12.*x(6)^2+18.8966*x(3)*x(6)+18.3815*x(3)*x(5)+19.9542*x(5)*x(6)+9.10511*x(3)^2+...
   
3.     9.44832*x(5)^2+8.85028*x(1)*x(6)+7.96096*x(1)*x(4)+11.5988*x(4)*x(6)+...
   
4.     3.84964*x(1)^2+4.42514*x(4)^2+3.71853*x(1)*x(2)+3.64590*x(1)*x(3)+...
   
5.     4.73504*x(1)*x(5)+5.94353*x(2)*x(3)+4.73504*x(2)*x(4)+6.53216*x(2)*x(5)+...
   
6.     6.53216*x(3)*x(4)+7.58636*x(4)*x(5)+7.58636*x(2)*x(6)+1.82295*x(2)^2;

复制代码

约束函数

1. function [c,ceq]=NonElipse(x)
   
2. c=-(12.*x(6)^2+18.8966*x(3)*x(6)+18.3815*x(3)*x(5)+19.9542*x(5)*x(6)+9.10511*x(3)^2+...
   
3.     9.44832*x(5)^2+8.85028*x(1)*x(6)+7.96096*x(1)*x(4)+11.5988*x(4)*x(6)+...
   
4.     3.84964*x(1)^2+4.42514*x(4)^2+3.71853*x(1)*x(2)+3.64590*x(1)*x(3)+...
   
5.     4.73504*x(1)*x(5)+5.94353*x(2)*x(3)+4.73504*x(2)*x(4)+6.53216*x(2)*x(5)+...
   
6.     6.53216*x(3)*x(4)+7.58636*x(4)*x(5)+7.58636*x(2)*x(6)+1.82295*x(2)^2);
   
7. ceq=x(2)^2-4*x(1)*x(3)+1;

复制代码

运行：

1. %% 调用函数
   
2. clc
   
3. x0=round(10*rand(1,6));                    % Starting guess
   
4. xdata=[0 0.099 0.15 0.1975 0.25 0.3265 0.411 0.5619 0.7265 0.951 1.126 1];
   
5. ydata=[0.9 1.1 0.9 1.1 0.9 1.1 0.975 0.9859 0.877 0.724 0.4152 0];
   
6. A=[];b=[];
   
7. Aeq=[0 0 1 0 1 1;1 0 0 1 0 1;1.126^2 1.126*.4152 .4152^2 1.126 .4152 1];
   
8. beq=zeros(3,1);
   
9. [x,fval,exitflag,output] = fmincon(@MainElipse,x0,A,b,Aeq,beq,[],[],@NonElipse)
   
10. f=@(x2,y2) 0.609856620799381*x2.^2+0.904625201010936*x2.*y2+...
    
11.     0.745399237479023*y2.^2-1.216767995653643*x2-...
    
12.     1.352310612333285*y2+0.606911374854262;
    
13. ezplot(f,[0,1.2])
    
14. hold on
    
15. plot(xdata,ydata,'r*')

复制代码

结果图形：

优化结果：

1. x =
   
2.   Columns 1 through 5
   
3.   -0.609891418321910  -0.904934994768739  -0.745586510558483   1.216995510549203   1.352690602785776
   
4.   Column 6
   
5.   -0.607104092227293
   
6. fval =
   
7.    0.010647923391415

复制代码

其中：设定椭圆方程为一般式：A x^2 + B x y + C y^2 + D x + E y + F ＝0
以数据∑(A x^2 + B x y + C y^2 + D x + E y + F)^2最小为目标函数，以定点为约束条件，转为非线性约束优化问题求解，优化结果中：[x1,x2,...,x6]=[A,B,C,D,E,F]

[ 本帖最后由 bainhome 于 2008-5-31 15:38 编辑 ]

linji526

原帖由 bainhome 于 2008/5/31 14:45 发表
已解决。
主函数function ff=MainElipse(x)
ff=12.*x(6)^2+18.8966*x(3)*x(6)+18.3815*x(3)*x(5)+19.9542*x(5)*x(6)+9.10511*x(3)^2+...
    9.44832*x(5)^2+8.85028*x(1)*x(6)+7.96096*x(1)*x(4)+11.5988*x(4)*x ...

赞!!!
不过主函数和约束函数的形式是如何得到的没有具体说明。

babecheer

非常感谢版主的热心帮忙，希望能更进一步的交流，我也想知道主函数和约束函数的形式是如何得到的，因为类似的曲线我后面还有三四条需要拟合，不知能否告知解决的思路或方法，再次感谢。

bainhome

8楼和17楼我的两个帖子其实讲得很清楚了：目标函数是将所有点值代入设定公式，先平方后求和，这个步骤用符号工具箱推一下即可。约束函数感觉其实c项可以不用写，直接c=[];即可，保险起见我还是写上了，就是让主目标函数值大于0，就这么简单。
运行部分中的Aeq和beq是等式约束，即：(0,1)、(1,0)、(1.126,.4152)三个指定点代入椭圆方程，其数值必须为0。