外推——基于 awk

Question

我需要以下方面的帮助：我有一个像 data.dat 一样的数据文件（由“\t”表格分​​隔的列）

    # y1    y2      y3      y4
    17.1685 21.6875 20.2393 26.3158 

这些是线性拟合的 4 个点的 x 值。四个 y 值是常数：0, 200, 400, 600。

我可以创建点对 (x,y): (x1,y1)=(17.1685,0), (x2,y2)=(21.6875,200), (x3,y3)=(20.2393,400), (x4,y4)=(26.3158,600) 的线性拟合。

现在我想对其中三个巴黎点进行线性拟合，(x1,y1), (x2,y2), (x3,y3) and (x2,y2), (x3,y3), (x4,y4) and (x1,y1), (x3,y3), (x4,y4) and (x1,y1), (x2,y2), (x4,y4).

如果我有这三个线性拟合的点，我想知道外推点的 x 值在这三个拟合点之外的值。

到目前为止，我有这个 awk 代码：

#!/usr/bin/awk -f

BEGIN{
  z[1] = 0;
  z[2] = 200;
  z[3] = 400;
  z[4] = 600;
}

{
  split($0,str,"\t");
  n = 0.0;

  for(i=1; i<=NF; i++)
  {
    centr[i] = str[i];
    n += 1.0;
  # printf("%d\t%f\t%.1f\t",i,centr[i],z[i]);
  }
  # print "";

  if (n > 2)
  {
    lsq(n,z,centr);
  }
}

function lsq(n,x,y)
{
  sx  = 0.0
  sy  = 0.0
  sxx = 0.0
  syy = 0.0
  sxy = 0.0
  eps = 0.0

  for (i=1;i<=n;i++)
  {
    sx  += x[i]
    sy  += y[i]
    sxx += x[i]*x[i]
    sxy += x[i]*y[i]
    syy += y[i]*y[i]
  }

  if ( (n==0) || ((n*sxx-sx*sx)==0) )
  {
    next;
  }
#   print "number of data points = " n;
  a = (sxx*sy-sxy*sx)/(n*sxx-sx*sx)
  b = (n*sxy-sx*sy)/(n*sxx-sx*sx)

  for(i=1;i<=n;i++)
  {
    ycalc[i] = a+b*x[i]
    dy[i]    = y[i]-ycalc[i]
    eps     += dy[i]*dy[i]
  }

  print "# Intercept =\t"a"
  print "# Slope     =\t"b"

  for (i=1;i<=n;i++)
  {
    printf("%8g %8g %8g \n",x[i],y[i],ycalc[i])
  }

} # function lsq()

所以，

    If we extrapolate to the place of 4th 
    0   17.1685   <--(x1,y1)
    200 21.6875   <--(x2,y2)
    400 20.2393   <--(x3,y3)
    600 22.7692 <<< (x4 = 600,y1 = 22.7692)

    If we extrapolate to the place of 3th
    0   17.1685   <--(x1,y1)
    200 21.6875   <--(x2,y2)
    400 23.6867 <<< (x3 = 400,y3 = 23.6867)
    600 26.3158   <--(x4,y4)

    0   17.1685
    200 19.35266 <<<
    400 20.2393
    600 26.3158

    0   18.1192 <<<
    200 21.6875
    400 20.2393
    600 26.3158

我目前的输出如下：

$> ./prog.awk data.dat
# Intercept =   17.4537
# Slope     =   0.0129968   
       0  17.1685  17.4537 
     200  21.6875  20.0531 
     400  20.2393  22.6525 
     600  26.3158  25.2518 

Answer 1

假设lsq 函数中的核心计算没问题（看起来差不多，但我没有仔细检查过），那么这将为您提供最适合的最小平方和线的斜率和截距输入数据集（参数 x、y、n）。我不确定我是否理解函数的结尾。

对于您的“取三个点并计算第四个”问题，最简单的方法是生成 4 个子集（从逻辑上讲，通过在四个调用中的每个调用中从一组四个中删除一个点），然后重新计算。

您需要调用另一个函数，该函数从lsq 获取线数据（斜率、截距）并在另一个 y 值处内插（外推）该值。这是一个直接的计算 (x = m * y + c)，但您需要确定您传入的 3 个集合中缺少哪个 y 值。

您可以通过从“平方和”、“总和”和“乘积之和”值中一次删除一个值，重新计算斜率、截距，然后再次计算缺失点。

（我还会观察到通常它是具有固定值 0、200、400、600 的 x 坐标，而 y 坐标是读取的值。但是，这只是方向问题，所以这并不重要。）

---

这至少是可行的代码。由于awk 自动在空白处拆分，因此您无需专门在选项卡上拆分；读取循环会考虑到这一点。

代码需要认真重构；里面有很多重复 - 然而，我也有我应该做的工作。

#!/usr/bin/awk -f
BEGIN{
  z[1] = 0;
  z[2] = 200;
  z[3] = 400;
  z[4] = 600;
}

{
  for (i = 1; i <= NF; i++)
  {
    centr[i] = $i
  }

  if (NF > 2)
  {
    lsq(NF, z, centr);
  }
}

function lsq(n, x, y)
{
  if (n == 0) return

  sx  = 0.0
  sy  = 0.0
  sxx = 0.0
  syy = 0.0
  sxy = 0.0

  for (i = 1; i <= n; i++)
  {
    print "x[" i "] = " x[i] ", y[" i "] = " y[i]
    sx  += x[i]
    sy  += y[i]
    sxx += x[i]*x[i]
    sxy += x[i]*y[i]
    syy += y[i]*y[i]
  }

  if ((n*sxx - sx*sx) == 0) return

#   print "number of data points = " n;
  a = (sxx*sy-sxy*sx)/(n*sxx-sx*sx)
  b = (n*sxy-sx*sy)/(n*sxx-sx*sx)

  for (i = 1; i <= n; i++)
  {
    ycalc[i] = a+b*x[i]
  }

  print "# Intercept = " a
  print "# Slope     = " b
  print "Line: x = " a " + " b " * y"

  for (i = 1; i <= n; i++)
  {
    printf("x = %8g, yo = %8g, yc = %8g\n", x[i], y[i], ycalc[i])
  }

  print ""
  print "Different subsets\n"

  for (drop = 1; drop <= n; drop++)
  {
    print "Subset " drop
    sx = sy = sxx = sxy = syy = 0
    j = 1
    for (i = 1; i <= n; i++)
    {
      if (i == drop) continue
      print "x[" j "] = " x[i] ", y[" j "] = " y[i]
      sx  += x[i]
      sy  += y[i]
      sxx += x[i]*x[i]
      sxy += x[i]*y[i]
      syy += y[i]*y[i]
      j++
    }
    if (((n-1)*sxx - sx*sx) == 0) continue
    a = (sxx*sy-sxy*sx)/((n-1)*sxx-sx*sx)
    b = ((n-1)*sxy-sx*sy)/((n-1)*sxx-sx*sx)
    print "Line: x = " a " + " b " * y"

    xt = x[drop]
    yt = a + b * xt;
    print "Interpolate: x = " xt ", y = " yt
  }
}

由于awk 没有提供从函数传回多个值的简单方法，也没有提供数组以外的结构（有时是关联的），因此它可能不是完成这项任务的最佳语言。另一方面，它可以完成这项工作。您可以将最小二乘计算捆绑在一个函数中，该函数返回一个包含斜率和截距的数组，然后使用它。轮到你探索选项了。

鉴于脚本lsq.awk 和输入文件lsq.data 显示，我得到显示的输出：

$ cat lsq.data
17.1685 21.6875 20.2393 26.3158
$ awk -f lsq.awk lsq.data
x[1] = 0, y[1] = 17.1685
x[2] = 200, y[2] = 21.6875
x[3] = 400, y[3] = 20.2393
x[4] = 600, y[4] = 26.3158
# Intercept = 17.4537
# Slope     = 0.0129968
Line: x = 17.4537 + 0.0129968 * y
x =        0, yo =  17.1685, yc =  17.4537
x =      200, yo =  21.6875, yc =  20.0531
x =      400, yo =  20.2393, yc =  22.6525
x =      600, yo =  26.3158, yc =  25.2518

Different subsets

Subset 1
x[1] = 200, y[1] = 21.6875
x[2] = 400, y[2] = 20.2393
x[3] = 600, y[3] = 26.3158
Line: x = 18.1192 + 0.0115708 * y
Interpolate: x = 0, y = 18.1192
Subset 2
x[1] = 0, y[1] = 17.1685
x[2] = 400, y[2] = 20.2393
x[3] = 600, y[3] = 26.3158
Line: x = 16.5198 + 0.0141643 * y
Interpolate: x = 200, y = 19.3526
Subset 3
x[1] = 0, y[1] = 17.1685
x[2] = 200, y[2] = 21.6875
x[3] = 600, y[3] = 26.3158
Line: x = 17.7985 + 0.0147205 * y
Interpolate: x = 400, y = 23.6867
Subset 4
x[1] = 0, y[1] = 17.1685
x[2] = 200, y[2] = 21.6875
x[3] = 400, y[3] = 20.2393
Line: x = 18.163 + 0.007677 * y
Interpolate: x = 600, y = 22.7692
$

编辑：在之前版本的答案中，子集乘以 n 而不是 (n-1)。修改后的输出中的值似乎与您的期望一致。剩余问题是表象的，而不是计算的。