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