Экстраполяция - на основе awk - PullRequest
6 голосов
/ 18 января 2012

Мне нужна помощь в следующем: у меня есть файл данных (столбцы, разделенные табличкой "\ t"), как этот data.dat

    # y1    y2      y3      y4
    17.1685 21.6875 20.2393 26.3158 

Это значения x 4 точек для линейной подгонки. Четыре значения 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 

1 Ответ

4 голосов
/ 18 января 2012

Если предположить, что вычисление ядра в функции lsq в порядке (оно выглядит правильно, но я его не изучал), то это дает вам наклон и перехват для линии наименьшей суммы квадратов, наиболее подходящей длянабор входных данных (параметры x, y, n).Я не уверен, что понимаю конечный конец функции.

Для вашей задачи «возьмите три точки и вычислите четвертую» задачу, самый простой способ - создать 4 подмножества (логически, удалив одну точку изнабор из четырех на каждый из четырех вызовов) и повторите расчет.

Вам необходимо вызвать другую функцию, которая берет данные линии (наклон, перехват) из lsq и интерполирует (экстраполирует) значение вдругое значение у.Это простой расчет (x = m * y + c), но вам нужно определить, какое значение y отсутствует в наборе 3, который вы передаете.

Вы можете «оптимизировать» (что означает «усложнить»)эта схема, отбрасывая одно значение за раз из значений «суммы квадратов» и «суммы» и «суммы произведений», пересчитывая наклон, перехват, а затем снова вычисляя недостающую точку.

(IТакже отметим, что обычно это будут координаты x с фиксированными значениями 0, 200, 400, 600, а координаты 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).Значения в пересмотренном выводе, кажется, согласуются с тем, что вы ожидаете.Остальные проблемы носят презентационный, а не вычислительный характер.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...