estrapolazione - awk basa

Question

ho bisogno di aiuto nella seguente: Ho un file di dati (colonne separate da "\ t" tabellare) come questo data.dat

# y1 y2  y3  y4 
    17.1685 21.6875 20.2393 26.3158 

Questi sono valori x di 4 punti per lineare in forma. I quattro valori y sono costanti: 0, 200, 400, 600.

Posso creare un adattamento lineare delle coppie di punti (x,y): (x1,y1)=(17.1685,0), (x2,y2)=(21.6875,200), (x3,y3)=(20.2393,400), (x4,y4)=(26.3158,600).

Ora vorrei fare un fit lineare su tre di questi punti di Parigi, (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).

Se ho questi tre punti con un adattamento lineare vorrei conoscere il valore del valore x del estrapolato punto di essere fuori di questi tre punti montati.

ho finora questo codice 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() 

Quindi,

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 

La mia uscita di corrente è il seguente:

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

Assumendo il calcolo di base in funzione lsq va bene (sembra giusto, ma non l'ho scrutato), poi questo ti dà la pendenza e io ntercept per la linea della somma minima dei quadrati che si adatta meglio al set di dati di input (parametri x, y, n). Non sono sicuro di capire la parte finale della funzione.

Per il "prendere tre punti e calcolare il quarto" problema, il modo più semplice è generare i 4 sottoinsiemi (in modo logico, eliminando un punto dall'insieme di quattro su ciascuna delle quattro chiamate) e ripetere il calcolo.

È necessario chiamare un'altra funzione che prende i dati di linea (pendenza, intercetta) da lsq e interpola (estrapola) il valore su un altro valore y. Questo è un calcolo straight-forward (x = m * y + c), ma è necessario determinare quale y valore è mancante dal set di 3 si passa.

Si potrebbe 'ottimizzare' (che significa 'complicato') questo schema facendo cadere un valore alla volta dai valori "somma di quadrati" e "somma" e "somma di prodotti", ricalcolo della pendenza, intercetta e quindi calcolando nuovamente il punto mancante.

(osserverò anche che normalmente sarebbero le coordinate x con i valori fissi 0, 200, 400, 600 e le coordinate y sarebbero i valori letti.Tuttavia, è solo una questione di orientamento, quindi non è fondamentale.)

---

Ecco il codice almeno plausibilmente lavoro. Poiché awk si divide automaticamente su uno spazio bianco, non è necessario dividerlo in modo specifico nelle schede; il ciclo di lettura tiene conto di questo.

Il codice richiede un serio refactoring; c'è un sacco di ripetizioni in esso - tuttavia, ho anche un lavoro che dovrei fare.

#!/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 
    } 
} 

Poiché awk non fornisce un modo semplice per passare indietro valori multipli da una funzione, né fornisce strutture diverse array (talvolta associativi), non è forse il miglior linguaggio per questo compito. D'altra parte, può essere fatto per fare il lavoro. Potresti essere in grado di raggruppare il calcolo dei minimi quadrati in una funzione che restituisce una matrice contenente la pendenza e l'intercetta, e quindi usarla. Il tuo turno per esplorare le opzioni.

Data la sceneggiatura lsq.awk e il file di input lsq.data mostrato, ottengo l'output mostrato:

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

Edit: Nella versione precedente della risposta, i sottoinsiemi si moltiplicavano da n invece di (n-1). I valori nell'output rivisto sembrano essere in accordo con ciò che ti aspetti. I problemi residui sono presentazionali, non computazionali.