Interpolazione lineare - Python

Question

Sono abbastanza nuovo in programmazione e ho pensato di provare a scrivere una funzione di interpolazione lineare.

dire che sono forniti dati come segue: x = [1, 2.5, 3.4, 5.8, 6] y = [2, 4, 5.8, 4.3, 4]

voglio progettare una funzione che interpolerà linearmente tra 1 e 2.5, da 2.5 a 3.4 e così via usando Python.

Ho provato a cercare attraverso http://docs.python.org/tutorial/ ma non riesco ancora a capirlo.

Answer 1

Come ho capito la tua domanda, vuoi scrivere una qualche funzione y = interpolate(x_values, y_values, x), che ti darà il valore a qualche x? L'idea di base quindi segue questi passaggi:

1. Trova gli indici dei valori in x_values che definiscono un intervallo contenente x. Ad esempio, per x=3 con gli elenchi di esempio, l'intervallo contenente sarebbe [x1,x2]=[2.5,3.4], e gli indici sarebbe i1=1, i2=2
2. calcolare la pendenza in questo intervallo (y_values[i2]-y_values[i1])/(x_values[i2]-x_values[i1]) (cioè dy/dx).
3. Il valore a x è ora il valore a x1 più la pendenza moltiplicata per la distanza da x1.

Sarà inoltre necessario decidere che cosa succede se x è al di fuori dell'intervallo di x_values, o si tratta di un errore, o si potrebbe interpolare "a ritroso", supponendo che il pendio è lo stesso del primo/ultimo intervallo.

Ti è stato d'aiuto o hai bisogno di un consiglio più specifico?

Answer 2

ho pensato una soluzione piuttosto elegante (IMHO), quindi non posso resistere a pubblicarlo:

from bisect import bisect_left 

class Interpolate(object): 
    def __init__(self, x_list, y_list): 
     if any(y - x <= 0 for x, y in zip(x_list, x_list[1:])): 
      raise ValueError("x_list must be in strictly ascending order!") 
     x_list = self.x_list = map(float, x_list) 
     y_list = self.y_list = map(float, y_list) 
     intervals = zip(x_list, x_list[1:], y_list, y_list[1:]) 
     self.slopes = [(y2 - y1)/(x2 - x1) for x1, x2, y1, y2 in intervals] 

    def __getitem__(self, x): 
     i = bisect_left(self.x_list, x) - 1 
     return self.y_list[i] + self.slopes[i] * (x - self.x_list[i]) 

I mappa per float in modo che la divisione intera (pitone < = 2.7), non entreranno in gioco e rovina le cose se x1, x2, y1 e y2 sono tutti numeri interi per alcuni iterval.

In __getitem__ sto approfittando del fatto che self.x_list è ordinato in ordine crescente utilizzando bisect_left a (molto) di trovare rapidamente l'indice del più grande elemento più piccolo di x in self.x_list.

Utilizzare la classe in questo modo:

i = Interpolate([1, 2.5, 3.4, 5.8, 6], [2, 4, 5.8, 4.3, 4]) 
# Get the interpolated value at x = 4: 
y = i[4] 

non ho affrontato le condizioni di confine a tutti qui, per semplicità. Così com'è, i[x] per x < 1 funzionerà come se la linea da (2,5, 4) a (1, 2) fosse stata estesa a meno infinito, mentre i[x] per x == 1 o x > 6 aumenterà un IndexError. Meglio sarebbe sollevare un IndexError in tutti i casi, ma questo è lasciato come esercizio per il lettore.:)

Answer 3

import scipy.interpolate 
y_interp = scipy.interpolate.interp1d(x, y) 
print y_interp(5.0) 

scipy.interpolate.interp1d esegue l'interpolazione lineare e può essere personalizzata per gestire le condizioni di errore.

Answer 4

La soluzione non ha funzionato in Python 2.7. Si è verificato un errore durante il controllo dell'ordine degli elementi x. Ho dovuto cambiare a codice a questo per farlo funzionare:

from bisect import bisect_left 
class Interpolate(object): 
    def __init__(self, x_list, y_list): 
     if any([y - x <= 0 for x, y in zip(x_list, x_list[1:])]): 
      raise ValueError("x_list must be in strictly ascending order!") 
     x_list = self.x_list = map(float, x_list) 
     y_list = self.y_list = map(float, y_list) 
     intervals = zip(x_list, x_list[1:], y_list, y_list[1:]) 
     self.slopes = [(y2 - y1)/(x2 - x1) for x1, x2, y1, y2 in intervals] 
    def __getitem__(self, x): 
     i = bisect_left(self.x_list, x) - 1 
     return self.y_list[i] + self.slopes[i] * (x - self.x_list[i]) 

Answer 5

Invece di estrapolare le estremità, è possibile restituire l'estensione della y_list. La maggior parte delle volte l'applicazione è ben educata e lo Interpolate[x] sarà nello x_list. Gli effetti (presumibilmente) lineari di estrapolazione dei fini possono indurre in errore a credere che i dati siano ben educati.

• Tornando risultato non lineare (delimitata dai contenuti di x_list e y_list) il comportamento del vostro programma può avvisare di un problema per i valori notevolmente al di fuori x_list. (Comportamento lineare va banane quando somministrato ingressi non lineari!)

• Tornando l'estensione della y_list per Interpolate[x] al di fuori del x_list significa anche che si conosce la portata del valore di uscita. Se estrapoli in base a x molto, molto meno di x_list[0] o x molto, molto maggiore di x_list[-1], il tuo risultato di ritorno potrebbe essere al di fuori dell'intervallo di valori che ti aspettavi.

  def __getitem__(self, x): 
      if x <= self.x_list[0]: 
       return self.y_list[0] 
      elif x >= self.x_list[-1]: 
       return self.y_list[-1] 
      else: 
       i = bisect_left(self.x_list, x) - 1 
       return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])