My c++ code is here:

//最小二乘曲面拟合.cpp

#include

#include

using namespace std;

//x[n] 存放给定数据点的X坐标。

//y[n] 存放给定数据点的Y坐标。

//z[n][m] 存放给定n*m个网点上函数值。

//n X坐标个数。

//m Y坐标个数。

//a[p][q] 返回二元拟合多项式的系数。

//p 拟合多项式中x的最高次为p-1。要求p<=min(n,20)。

//q 拟合多项式中y的最高次为q-1。要求q<=min(n,20)。

//dt[3] 分别返回误差平方和，误差绝对值之和与误差绝对值的最大值。

void pir2(double x[], double y[], double z[], int n, int m,

          double a[], int p, int q, double dt[])

{

    int i,j,k,l,kk;

  double apx[20],apy[20],bx[20],by[20],u[20][20];

  double t[20],t1[20],t2[20],d1,d2,g,g1,g2;

  double x2,dd,y1,x1,*v;

  v = new double[20*m];

  for (i=0; i<=p-1; i++)

  {

          l=i*q;

      for (j=0; j<=q-1; j++) a[l+j]=0.0;

  }

  if (p>n) p=n;

  if (p>20) p=20;

  if (q>m) q=m;

  if (q>20) q=20;

  d1=1.0*n; apx[0]=0.0;

  for (i=0; i<=n-1; i++)  apx[0]=apx[0]+x[i];

  apx[0]=apx[0]/d1;

  for (j=0; j<=m-1; j++)

  {

          v[j]=0.0;

      for (i=0; i<=n-1; i++)  v[j]=v[j]+z[i*m+j];

      v[j]=v[j]/d1;

  }

  if (p>1)

  {

          d2=0.0; apx[1]=0.0;

      for (i=0; i<=n-1; i++)

      {

               g=x[i]-apx[0];

          d2=d2+g*g;

          apx[1]=apx[1]+x[i]*g*g;

      }

      apx[1]=apx[1]/d2;

      bx[1]=d2/d1;

      for (j=0; j<=m-1; j++)

      {

               v[m+j]=0.0;

          for (i=0; i<=n-1; i++)

          {

                     g=x[i]-apx[0];

              v[m+j]=v[m+j]+z[i*m+j]*g;

          }

          v[m+j]=v[m+j]/d2;

      }

      d1=d2;

  }

  for (k=2; k<=p-1; k++)

  {

          d2=0.0; apx[k]=0.0;

      for (j=0; j<=m-1; j++) v[k*m+j]=0.0;

      for (i=0; i<=n-1; i++)

      {

               g1=1.0; g2=x[i]-apx[0];

          for (j=2; j<=k; j++)

          {

                     g=(x[i]-apx[j-1])*g2-bx[j-1]*g1;

              g1=g2; g2=g;

          }

          d2=d2+g*g;

          apx[k]=apx[k]+x[i]*g*g;

          for (j=0; j<=m-1; j++) v[k*m+j]=v[k*m+j]+z[i*m+j]*g;

      }

      for (j=0; j<=m-1; j++) v[k*m+j]=v[k*m+j]/d2;

      apx[k]=apx[k]/d2;

      bx[k]=d2/d1;

      d1=d2;

  }

  d1=m; apy[0]=0.0;

  for (i=0; i<=m-1; i++) apy[0]=apy[0]+y[i];

  apy[0]=apy[0]/d1;

  for (j=0; j<=p-1; j++)

  {

          u[j][0]=0.0;

      for (i=0; i<=m-1; i++)  u[j][0]=u[j][0]+v[j*m+i];

        u[j][0]=u[j][0]/d1;

  }

  if (q>1)

  {

          d2=0.0; apy[1]=0.0;

      for (i=0; i<=m-1; i++)

      {

               g=y[i]-apy[0];

          d2=d2+g*g;

          apy[1]=apy[1]+y[i]*g*g;

      }

      apy[1]=apy[1]/d2;

      by[1]=d2/d1;

      for (j=0; j<=p-1; j++)

          {

               u[j][1]=0.0;

          for (i=0; i<=m-1; i++)

          {

                     g=y[i]-apy[0];

                  u[j][1]=u[j][1]+v[j*m+i]*g;

          }

            u[j][1]=u[j][1]/d2;

      }

      d1=d2;

  }

  for (k=2; k<=q-1; k++)

  {

          d2=0.0; apy[k]=0.0;

        for (j=0; j<=p-1; j++) u[j][k]=0.0;

      for (i=0; i<=m-1; i++)

      {

               g1=1.0;

          g2=y[i]-apy[0];

          for (j=2; j<=k; j++)

          {

                     g=(y[i]-apy[j-1])*g2-by[j-1]*g1;

              g1=g2; g2=g;

          }

          d2=d2+g*g;

          apy[k]=apy[k]+y[i]*g*g;

          for (j=0; j<=p-1; j++) u[j][k]=u[j][k]+v[j*m+i]*g;

      }

      for (j=0; j<=p-1; j++) u[j][k]=u[j][k]/d2;

      apy[k]=apy[k]/d2;

      by[k]=d2/d1;

      d1=d2;

  }

  v[0]=1.0; v[m]=-apy[0]; v[m+1]=1.0;

  for (i=0; i<=p-1; i++)

  for (j=0; j<=q-1; j++) a[i*q+j]=0.0;

  for (i=2; i<=q-1; i++)

  {

          v[i*m+i]=v[(i-1)*m+(i-1)];

      v[i*m+i-1]=-apy[i-1]*v[(i-1)*m+i-1]+v[(i-1)*m+i-2];

      if (i>=3)

      for (k=i-2; k>=1; k--)

        v[i*m+k]=-apy[i-1]*v[(i-1)*m+k]+

                 v[(i-1)*m+k-1]-by[i-1]*v[(i-2)*m+k];

      v[i*m]=-apy[i-1]*v[(i-1)*m]-by[i-1]*v[(i-2)*m];

  }

  for (i=0; i<=p-1; i++)

  {

          if (i==0) { t[0]=1.0; t1[0]=1.0;}

      else

      {

               if (i==1)

          {

                     t[0]=-apx[0]; t[1]=1.0;

              t2[0]=t[0]; t2[1]=t[1];

          }

          else

          {

                     t[i]=t2[i-1];

              t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];

              if (i>=3)

              for (k=i-2; k>=1; k--)

                t[k]=-apx[i-1]*t2[k]+t2[k-1]

                     -bx[i-1]*t1[k];

              t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];

              t2[i]=t[i];

              for (k=i-1; k>=0; k--)

              { t1[k]=t2[k]; t2[k]=t[k];}

          }

      }

      for (j=0; j<=q-1; j++)

      for (k=i; k>=0; k--)

      for (l=j; l>=0; l--)

           a[k*q+l]=a[k*q+l]+u[i][j]*t[k]*v[j*m+l];

  }

  dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;

  for (i=0; i<=n-1; i++)

  {

          x1=x[i];

      for (j=0; j<=m-1; j++)

      {

               y1=y[j];

          x2=1.0; dd=0.0;

          for (k=0; k<=p-1; k++)

          {

                     g=a[k*q+q-1];

              for (kk=q-2; kk>=0; kk--) g=g*y1+a[k*q+kk];

              g=g*x2; dd=dd+g; x2=x2*x1;

          }

          dd=dd-z[i*m+j];

          if (fabs(dd)>dt[2]) dt[2]=fabs(dd);

          dt[0]=dt[0]+dd*dd;

          dt[1]=dt[1]+fabs(dd);

      }

  }

  delete[] v;

  return;

}

//最小二乘曲面拟合例

#include

#include

#include

// #include "最小二乘曲面拟合.cpp"

using namespace std;

int main()

{

    int i,j;

  double x[11],y[21],z[11][21],a[6][5],dt[3];

  for (i=0; i<=10; i++) x[i]=0.2*i;

  for (i=0; i<=20; i++) y[i]=0.1*i;

  for (i=0; i<=10; i++)

  for (j=0; j<=20; j++)

      z[i][j]=exp(x[i]*x[i]-y[j]*y[j]);

  pir2(x,y,&z[0][0],11,21,&a[0][0],6,5,dt);

  cout <<"二元拟合多项式系数矩阵 : " <<endl;

  for (i=0; i<=5; i++)

  {

          for (j=0; j<=4; j++)

          cout <<setw(14) <<a[i][j];

      cout <<endl;

  }

  cout <<"误差平方和 = " <<dt[0] <<endl;

    cout <<"误差绝对值和 = " <<dt[1] <<endl;

    cout <<"误差绝对值最大值 = " <<dt[2] <<endl;

  return 0;

}

结果:

二元拟合多项式系数矩阵 : 0.888696 0.096095 -1.38677 0.903053 -0.171411 8.59716 0.929614 -13.4154 8.73604 -1.65821 -45.7529 -4.94728 71.3951 -46.492 8.82477 85.1474 9.20702 -132.868 86.5229 -16.4231 -62.8532 -6.79634 98.0792 -63.8686 12.1231 16.99 1.83714 -26.5121 17.2645 -3.27702 误差平方和 = 5.80139 误差绝对值和 = 25.3761 误差绝对值最大值 = 0.555972

这个结果应该是正确的.

rust 代码:

#[no_mangle]

extern "C" fn fit_surface_least_squares_m( n:usize, m:usize,mut p:usize,mut q:usize,x:&[f64;100],y:&[f64;100],z:&[[f64;100];100],a:&mut [[f64;20];20],dt:&mut [f64;3])->bool

{

  //let mut l:usize;

  //let mut kk:usize;

let mut apx=[0.0;20];

let mut apy=[0.0;20];

let mut bx=[0.0;20];

let mut by=[0.0;20];

let mut t=[0.0;20];

let mut t1=[0.0;20];

let mut t2=[0.0;20];

//let mut d2:f64;

let mut d=[0.0;2];

let mut g=[0.0;3];

//let mut g1:f64;

//let mut g2:f64;

let mut x2:f64;

let mut dd:f64;

let mut y1:f64;

let mut x1:f64;

let mut v=[[0.0;100];20];

let mut u=[[0.0;20];20];

  //matrix<_Ty> v(20,m), u(20,20);



  //int p = a.GetRowNum();       //拟合多项式中变量x的最高次数加1

  //并要求p<=n且p<=20; 否则在本函数中自动取p=min{n,20}处理



  //int q = a.GetColNum();        //拟合多项式中变量y的最高次数加1

  //并要求q<=m且q<=20; 否则在本函数中自动取q=min{m,20}处理



for i in 0..p

{

  for j in 0..q

  {

        a[i][j]=0.0;

  }

}   



if p>n {p=n;}

if p>20 {p=20;}

if q>m {q=m;}

if q>20 {q=20;}



let mut xx:f64=0.0;





for i in 0..n

{

  xx+=  x[i] / (n as f64);

}

let mut yy:f64=0.0;

for i in 0..m

{

  yy+= y[i] / (m as f64);

}

d[0]=n as f64;

apx[0]=0.0;

for i in 0..n

{

  apx[0]+=x[i] - xx;

}

apx[0] = apx[0] / d[0];

for j in 0..m

{

        v[0][j]=0.0;

    for i in 0..n

    {

        v[0][j]=v[0][j]+z[i][j];

    }  

    v[0][j]=v[0][j]/d[0];

}

if p>1

{

        d[1]=0.0;

        apx[1]=0.0;

    for i in 0..n

    {

             g[0]=x[i]-xx-apx[0];

        d[1]+=g[0]*g[0];

        apx[1]+=(x[i]-xx)*g[0]*g[0];

    }

    apx[1]=apx[1]/d[1];

    bx[1]=d[1]/d[0];

    for j in 0..m

    {

             v[1][j]=0.0;

        for i in 0..n

        {

                   g[0]=x[i]-xx-apx[0];

            v[1][j]=v[1][j]+z[i][j]*g[0];

        }

        v[1][j]/=d[1];

    }

    d[0]=d[1];

}

for k in 2..p

{

        d[1]=0.0;

        apx[k]=0.0;

    for j in 0..m

    {

        v[k][j]=0.0;

    }

    for i in 0..n

    {

             g[1]=1.0;

             g[2]=x[i]-xx-apx[0];



        for j in 2..=k

        {

                   g[0]=(x[i]-xx-apx[j-1])*g[2]-bx[j-1]*g[1];

            g[1]=g[2];

            g[2]=g[0];

        }

        d[1]+=g[0]*g[0];

        apx[k]=apx[k]+(x[i]-xx)*g[0]*g[0];

        for j in 0..m    

        {

              v[k][j]+=z[i][j]*g[0];

        }

    }

    for j in 0..m

    {

        v[k][j]/=d[1];

    }  

    apx[k]=apx[k]/d[1];

    bx[k]=d[1]/d[0];

    d[0]=d[1];

}

d[0]=m as f64;

  apy[0]=0.0;

for i in 0..m

{

  apy[0]+=y[i]-yy;

}

apy[0]=apy[0]/d[0];

for j in 0..p

{

        u[j][0]=0.0;

    for i in 0..m

    {

        u[j][0]+=v[j][i];

    }  

        u[j][0]/=d[0];

}

if q>1

{

        d[1]=0.0;

        apy[1]=0.0;

    for i in 0..m

    {

             g[0]=y[i]-yy-apy[0];

        d[1]+=g[0]*g[0];

        apy[1]+=(y[i]-yy)*g[0]*g[0];

    }

    apy[1]=apy[1]/d[1];

    by[1]=d[1]/d[0];

    for j in 0..p

        {

             u[j][1]=0.0;

        for i in 0..m

        {

                   g[0]=y[i]-yy-apy[0];

                   u[j][1]+=v[j][i]*g[0];

        }

             u[j][1]/=d[1];

    }

    d[0]=d[1];

}

for k in 2..q

{

        d[1]=0.0;

        apy[k]=0.0;

        for j in 0..p

        {

             u[j][k]=0.0;

        }

    for i in 0..m

    {

             g[1]=1.0;

        g[2]=y[i]-yy-apy[0];

        for j in 2..=k

        {

                   g[0]=(y[i]-yy-apy[j-1])*g[2]-by[j-1]*g[1];

            g[1]=g[2];

                   g[2]=g[0];

        }

        d[1]+=g[0]*g[0];

        apy[k]=apy[k]+(y[i]-yy)*g[0]*g[0];

        for j in 0..p

        {

              u[j][k]+=v[j][i]*g[0];

        }    

    }

    for j in 0..p

    {

        u[j][k]/=d[1];

    }  

    apy[k]/=d[1];

    by[k]=d[1]/d[0];

    d[0]=d[1];

}

v[0][0]=1.0;

  v[1][0]=-apy[0];

  v[1][1]=1.0;

for i in 0..p

{

  for j in 0..q

  {

        a[i][j]=0.0;

  }    

}



for i in 2..q

{

        v[i][i]=v[i-1][i-1];

    v[i][i-1]=-apy[i-1]*v[i-1][i-1]+v[i-1][i-2];

    if i>=3

    {

        for k in (1..=(i-2)).rev()

        {

              v[i][k]=-apy[i-1]*v[i-1][k]+v[i-1][k-1]-by[i-1]*v[i-2][k];

        }

       

    }

     

    v[i][0]=-apy[i-1]*v[i-1][0]-by[i-1]*v[i-2][0];

}

for i in 0..p

{

        if i==0

        {

             t[0]=1.0;

             t1[0]=1.0;

        }

    else

    {

             if i==1

        {

                   t[0]=-apx[0];

                   t[1]=1.0;

            t2[0]=t[0];

                   t2[1]=t[1];

        }

        else

        {

                   t[i]=t2[i-1];

            t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];

            if i>=3

            {

              for k in (1..=(i-2)).rev()

              {

                    t[k]=-apx[i-1]*t2[k]+t2[k-1]-bx[i-1]*t1[k];

                    

              }

               

            }

             

            t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];

            t2[i]=t[i];

            for k in (0..=(i-1)).rev()

            {

                         t1[k]=t2[k];

                         t2[k]=t[k];

                   }

        }

    }

    for j in 0..q

    {

        for k in (0..=i).rev()

        {

              for l in (0..=j).rev()

                   {

                         a[k][l]+=u[i][j]*t[k]*v[j][l];

                   }

        }

       

    }

     

}

for i in 0..3

{

  dt[i]=0.0;

}

  //计算dt

//dt[0]存放拟合多项式与数据点误差的平方和

//dt[1]存放拟合多项式与数据点误差的绝对值之和

//dt[2]存放拟合多项式与数据点误差绝对值的最大值

for i in 0..n

{

        x1=x[i]-xx;

    for j in 0..m

    {

             y1=y[j]-yy;

        x2=1.0;

             dd=0.0;

        for k in 0..p

        {

                   g[0]=a[k][q-1];

            for kk in (0..=(q-2)).rev()

            {

              g[0]=g[0]*y1+a[k][kk];

            }

            g[0]=g[0]*x2;

                   dd+=g[0];

                   x2*=x1;

        }

        dd-=z[i][j];

        if dd.abs()>dt[2]

        {

              dt[2]=dd.abs();

        }

        dt[0]=dt[0]+dd*dd;

        dt[1]=dt[1]+dd.abs();

    }

}

true

}

When I compile it to c dll and call it by RubyFFI got the following result:

系数矩阵a=[[12.610828495549994, -6.4031174330927225, -48.421984673330485, 7.374037165505316, 38.79011866159077], [-6.304552141689332, 3.2011209842825186, 24.207682095178615, -3.68651447611775, -19.39240754646785], [-1.5801195738843747, 0.8023018625207973, 6.06720849590951, -0.9239567779291769, -4.860348849570485], [1.3068055530349345, -0.6635273345643454, -5.01776060806868, 0.7641395424233058, 4.019652038542775], [0.043835634614448844, -0.022257436637848578, -0.16831633450534478, 0.025632384020965664, 0.13483566673677283], [-0.12004056416126747, 0.06095030388610434, 0.46092152946540504, -0.07019234159007005, -0.3692372574619461]]

误差平方和 = 14546.87281940942 误差绝对值和 =1105.0604796924251 误差绝对值最大值 =38.4058183240485

The ruby FFI function as following:

module Qin

extend FFI::Library

….

attach_function :fit_surface_least_squares_m,[:int,:int,:int,:int,:pointer,:pointer,:pointer,:pointer,:pointer],:bool

def self.fit_surface_least_squares(p,q,x_coord,y_coord,z_coord)

  #extern "C" fn fit_surface_least_squares_m( n:usize, m:usize,mut p:usize,mut q:usize,x:&[f64;100],y:&[f64;100],z:&[[f64;100];100],a:&mut [[f64;100];100],dt:&mut [f64;3])->bool

 #int p = a.GetRowNum();        //拟合多项式中变量x的最高次数加1

 #int q = a.GetColNum();          //拟合多项式中变量y的最高次数加1



   n=x_coord.length

   m=y_coord.length



   pointer_x=FFI::MemoryPointer.new :double,100

 pointer_y=FFI::MemoryPointer.new :double,100

 pointer_x.put_array_of_double(0,x_coord)

 pointer_y.put_array_of_double(0,y_coord)





pointer_z=FFI::MemoryPointer.new :double,100*100

pointer_a=FFI::MemoryPointer.new :double,20*20

pointer_dt=FFI::MemoryPointer.new :double,3

z_coord.each_with_index do |e,i|

    pointer_z.put_array_of_double(i*100*8,e)  ##8-->f64--float,n=100

end



 if n<=100 and m<=100

    if self.fit_surface_least_squares_m(n,m,p,q,pointer_x,pointer_y,pointer_z,pointer_a,pointer_dt)

       a=[]

       (0..(p-1)).each do |i|

           a[i]=pointer_a.get_array_of_double(i*20*8,q)

       end

       dt=pointer_dt.get_array_of_double(0,3)

       return [a,dt]

    end

 else

    nil

 end

end

…

end #end module

x_coord=(0..10).to_a.map{|e| e*0.2}

y_coord=(0..20).to_a.map{|e| e*0.1}

z_coord=Array.new(21) { Array.new(11) { |i| 0.0} }

for j in 0..20

  for i in 0..10

  z_coord[j][i]=Math.exp(x_coord[i]**2-y_coord[j]**2)

  end

end

#puts z_coord.inspect

a=Qin.fit_surface_least_squares(p=6,q=5,x_coord,y_coord,z_coord)

puts a[0].inspect

puts a[1]

转换似乎正确，结果却错误。这是为什么呢？

n=1000000
m=10
x=runif(n)#n个点的X坐标
y=runif(n)#n个点的Y坐标，因为是随机数，所以与x不同
z=x+y-x*y+x^2-y^2+x*y^2-x^2*y+x^3+y^3+rnorm(1000000)#带误差的
lm(z~1+x+y+I(x*y)+I(x^2)+I(y^2)+I(x^3)+I(x^2*y)+I(x*y^2)+I(y^3))#（还可以继续加但我懒得写了）

Call:
lm(formula = z ~ 1 + x + y + I(x * y) + I(x^2) + I(y^2) + I(x^3) + 
    I(x^2 * y) + I(x * y^2) + I(y^3))

Coefficients:
(Intercept)            x            y     I(x * y)       I(x^2)       I(y^2)  
   0.006309     1.000321     0.942730    -0.862455     0.973187    -0.975042  
     I(x^3)   I(x^2 * y)   I(x * y^2)       I(y^3)  
   1.017611    -1.052587     0.928569     1.018850