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竟然可以这么简单。为啥我以前不知道？
|
|
