h****a 发帖数: 580 | 1 向各位牛人请教一个问题:
现有一套数据(x[m],y[n],z[m,n]),其中x[m]和y[n]是参数,z[m,n]是(x[m],y[n])
相对应的值,m从1到M,n从1到N,一共有M*N个数据点.
根据我的观测和初步检验,我认为下面的方程可以很好地拟合(fit)这些数据:
z=f(x,y)=a11*x^2*y^2+a12*x^2*y+a13*x^2+a21*x*y^2+a22*x*y+a23*x+a31*y^2+a32*y
+a33
现在的问题是,如何准确地得到方程中的系数coeffients(a11,a12,a13,a21,a22,a23,
a31,a32,a33)?
据我所知,对于二维数据(x[i],y[i]), 很容易用方程y=f(x)来拟合,有很多工具(
gnuplot,xmgrace)可以做到。但是对于三维数据(x[m],y[n],z[m,n]),请问有什么现
成的工具可以用来拟合?若是没有现成的工具,我需要自己写程序,该如何着手?
我在网上搜索“least-squares fitting”,找到一些说明和子程序,但全都是用y=f(x
)来拟合数据(x,y),没有用z= | l*****d 发帖数: 359 | 2 自己写程序应该不难吧。
求 E(x,y,z)=sigma(square(z-f(x,y))的最小值 , sigma对你所有已知(x[m],y[n],z[m,n])的求和
用拉格朗日法,解 E对每个系数aij求导等于零的方程组。 | r**u 发帖数: 130 | 3 lease square 拟和,自己构造X矩阵。
系数=inv(X'X)X'Y | h****a 发帖数: 580 | 4 我数学基础比较差,尤其是一遇到矩阵就头疼。老大能否给个网页链接,讲这个“拉格
朗日法”的原理的,最好有子程序。
先发个包子。
[m,n])的求和
【在 l*****d 的大作中提到】 : 自己写程序应该不难吧。 : 求 E(x,y,z)=sigma(square(z-f(x,y))的最小值 , sigma对你所有已知(x[m],y[n],z[m,n])的求和 : 用拉格朗日法,解 E对每个系数aij求导等于零的方程组。
| h****a 发帖数: 580 | 5 请问这个X矩阵包含些什么?能否讲详细些?
先发个包子。
【在 r**u 的大作中提到】 : lease square 拟和,自己构造X矩阵。 : 系数=inv(X'X)X'Y
| l*****d 发帖数: 359 | 6 名字不要管他, 道理很简单 ,一次导数为0的点就是最小值。那个E其实就是观测值z
和你拟合的f(x,y)之间的的均方差函数,把所有的x[m],y[n],z[m,n]带入E,你就得到
一个以系数aij为变量的多元函数E(a11,a12,a13,a21,...)。n个系数的话,这个函数就
有n个变量(为了方便,把你的系数aij换成x1,x2...,xn)。 分别对这n个变量求导,你
就得到n个线性方程。这个线性方程组中xi的系数可以写成n阶矩阵。 求xi向量就等于
求解一个矩阵方程。自己编程求矩阵的每个元素吧, 并不难。写个matlab function A
= get_coef(x[m],y[n],z[m,n])。A是N阶方阵,N为你系数个数,在你的case中N=9
【在 h****a 的大作中提到】 : 我数学基础比较差,尤其是一遇到矩阵就头疼。老大能否给个网页链接,讲这个“拉格 : 朗日法”的原理的,最好有子程序。 : 先发个包子。 : : [m,n])的求和
| r**u 发帖数: 130 | 7
k=0;
for i=1:m
for j=1:n
k++;
X(k,1)=x(i)^2*y(j)^2;
X(k,2)=.............;
X(k,3)=.............;
...
y(k)=z(i,j)
X is the matrix, y is the output vector
【在 h****a 的大作中提到】 : 请问这个X矩阵包含些什么?能否讲详细些? : 先发个包子。
| l*****d 发帖数: 359 | 8 为何要写inv(X'X)X'Y?
难道不是 inv(X'X)X'Y = inv(X)inv(X')X'Y = inv(X)Y?
【在 r**u 的大作中提到】 : lease square 拟和,自己构造X矩阵。 : 系数=inv(X'X)X'Y
| D*******a 发帖数: 3688 | 9 x不是方阵
【在 l*****d 的大作中提到】 : 为何要写inv(X'X)X'Y? : 难道不是 inv(X'X)X'Y = inv(X)inv(X')X'Y = inv(X)Y?
| d*****l 发帖数: 8441 | 10 伪逆(pseudo-inverse),因为精确拟合不可能,只能有最小二乘意义上的精确度。
确实是X^(+)=(X'X)^(-1)X'。
【在 l*****d 的大作中提到】 : 为何要写inv(X'X)X'Y? : 难道不是 inv(X'X)X'Y = inv(X)inv(X')X'Y = inv(X)Y?
| d*****l 发帖数: 8441 | 11 不尽然,高维情况下,“一次导数为0的点”也可能是鞍点。
不过,在least-square sense下,应当是成立的啦。
z
A
【在 l*****d 的大作中提到】 : 名字不要管他, 道理很简单 ,一次导数为0的点就是最小值。那个E其实就是观测值z : 和你拟合的f(x,y)之间的的均方差函数,把所有的x[m],y[n],z[m,n]带入E,你就得到 : 一个以系数aij为变量的多元函数E(a11,a12,a13,a21,...)。n个系数的话,这个函数就 : 有n个变量(为了方便,把你的系数aij换成x1,x2...,xn)。 分别对这n个变量求导,你 : 就得到n个线性方程。这个线性方程组中xi的系数可以写成n阶矩阵。 求xi向量就等于 : 求解一个矩阵方程。自己编程求矩阵的每个元素吧, 并不难。写个matlab function A : = get_coef(x[m],y[n],z[m,n])。A是N阶方阵,N为你系数个数,在你的case中N=9
| l*****d 发帖数: 359 | 12 我说的仅限于最小二乘法的误差函数。否则还可能是最大值呢
【在 d*****l 的大作中提到】 : 不尽然,高维情况下,“一次导数为0的点”也可能是鞍点。 : 不过,在least-square sense下,应当是成立的啦。 : : z : A
| l*****d 发帖数: 359 | 13 靠,LS fitting竟然可以这么简单。为啥我以前不知道?
【在 d*****l 的大作中提到】 : 伪逆(pseudo-inverse),因为精确拟合不可能,只能有最小二乘意义上的精确度。 : 确实是X^(+)=(X'X)^(-1)X'。
| d*****l 发帖数: 8441 | 14 矩阵论里面应当讲的。
【在 l*****d 的大作中提到】 : 靠,LS fitting竟然可以这么简单。为啥我以前不知道?
|
|