外挿 - awkの私は、次のヘルプが必要

Question

をもと：私は（「\トン」表形式で区切られた列）のデータファイルを持っているこのdata.dat

# y1 y2  y3  y4 
    17.1685 21.6875 20.2393 26.3158 

のようなこれらは、線形のために4点の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)と作成することができます。

今

私はこれらのポイントパリの3に線形近似をしたいと思い、(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).

私は線形近似との点のこれらの3を持っている場合、私は外挿のx値の値を知りたいのですがポイントはこれらの3つのフィッティングポイントから外れています。

私はこれまでのところ、この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 

私の現在の出力は以下の通りです：

lsq機能の中核計算を想定し

$> ./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

OKです（それは正しいと思われますが、私はそれを精査していません）。入力データセット（パラメータx、y、n）のベストフィットの二乗和の最小値を求めます。私は関数の末尾を理解しているかどうかはわかりません。

「3つのポイントを取って4番目に計算する」問題については、最も簡単な方法は4つのサブセットを生成することです（論理的には、4つの呼び出しごとに4つのセットから1つのポイントを削除します）。

lsqからラインデータ（スロープ、インターセプト）を受け取り、別のy値で値を補間（外挿）する別の関数を呼び出す必要があります。これは簡単な計算（x = m * y + c）ですが、あなたが渡した3つのセットのうちどれが欠けているかを判断する必要があります。

このスキームを「最適化」することができます「二乗和」と「合計」と「積の合計」の値から一度に、傾きを再計算し、切片を再計算し、欠けている点を再び計算する。

（通常は、固定値0、200、400、600のx座標で、y座標は読み取られた値になりますが、したがって、それは重要ではありません。）

---

ここには少なくとも偶発的に動作するコードがあります。 awkは空白で自動的に分割されるため、タブを分割する必要はありません。読み取りループはこれを考慮に入れます。

コードには深刻なリファクタリングが必要です。そこには繰り返しのトンがあります - しかし、私はまた私がやらなければならない仕事があります。それはおそらく、このタスクのための最良の言語ではない、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 
    } 
} 

関数から複数の値をバック通過する簡単な方法を提供しません。また、（時には連想）配列以外の構造を提供しません。その一方で、それは仕事をすることができます。 Least Squaresの計算を、勾配と切片を含む配列を返す関数にバンドルして、それを使用することができます。オプションを探そう。

$ 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)を乗じたスクリプトlsq.awkと示された入力ファイルlsq.data考える

は、私が示した出力が得られます。修正された出力の値は、あなたが期待しているものと一致しているようです。残りの問題は、計算的ではなく、表現的である。