Kleinste Quadrate Minimierung Komplexe Zahlen

Question

Ich habe meine Matlab verwendet, aber es ist meine Vision, schließlich zu all meinen Analysen in Python zu wechseln, da es eine tatsächliche Programmiersprache und ein paar andere Gründe ist.

Das jüngste Problem, das ich anzugehen versucht habe, ist die Minimierung der kleinsten Quadrate komplexer Daten. Ich bin ein Ingenieur und wir beschäftigen uns ziemlich oft mit komplexer Impedanz, und ich versuche Kurvenanpassung zu verwenden, um ein einfaches Schaltungsmodell an gemessene Daten anzupassen.

Die Impedanzgleichung ist wie folgt:

Z (w) = 1/(1/R + j * B * C) + j * w * L

Ich versuche dann zu finden, die Werte von R, C und L, so dass die Kurve der kleinsten Quadrate gefunden wird.

Ich habe versucht, mit dem Optimierungspaket wie optimize.curve_fit oder optimize.leastsq, aber sie arbeiten nicht mit komplexen Zahlen.

Ich versuchte dann, meine Restfunktion die Größe der komplexen Daten zurückgeben, aber das hat auch nicht funktioniert.

Answer 1

Letztlich ist das Ziel, den absoluten Wert der Summe der quadrierten Unterschiede zwischen dem Modell zu reduzieren und die beobachtete Z:

abs(((model(w, *params) - Z)**2).sum()) 

Mein original answer vorgeschlagen leastsq auf eine residuals Funktion der Anwendung, die eine skalare zurückgibt die Summe der Quadrate der realen und imaginären Unterschiede darstellen:

def residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    return diff.real**2 + diff.imag**2 

Mike Sulzer suggested mit einer Restfunktion, die einen Vektor von Floats zurückgibt. Hier

ist ein Vergleich des Ergebnisses dieser restlichen Funktionen:

from __future__ import print_function 
import random 
import numpy as np 
import scipy.optimize as optimize 
j = 1j 

def model1(w, R, C, L): 
    Z = 1.0/(1.0/R + j*w*C) + j*w*L 
    return Z 

def model2(w, R, C, L): 
    Z = 1.0/(1.0/R + j*w*C) + j*w*L 
    # make Z non-contiguous and of a different complex dtype 
    Z = np.repeat(Z, 2) 
    Z = Z[::2] 
    Z = Z.astype(np.complex64) 
    return Z 

def make_data(R, C, L): 
    N = 10000 
    w = np.linspace(0.1, 2, N) 
    Z = model(w, R, C, L) + 0.1*(np.random.random(N) + j*np.random.random(N)) 
    return w, Z 

def residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    return diff.real**2 + diff.imag**2 

def MS_residuals(params, w, Z): 
    """ 
    https://stackoverflow.com/a/20104454/190597 (Mike Sulzer) 
    """ 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    z1d = np.zeros(Z.size*2, dtype=np.float64) 
    z1d[0:z1d.size:2] = diff.real 
    z1d[1:z1d.size:2] = diff.imag 
    return z1d 

def alt_residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    return diff.astype(np.complex128).view(np.float64) 

def compare(*funcs): 
    fmt = '{:15} | {:37} | {:17} | {:6}' 
    header = fmt.format('name', 'params', 'sum(residuals**2)', 'ncalls') 
    print('{}\n{}'.format(header, '-'*len(header))) 
    fmt = '{:15} | {:37} | {:17.2f} | {:6}' 
    for resfunc in funcs: 
     # params, cov = optimize.leastsq(resfunc, p_guess, args=(w, Z)) 
     params, cov, infodict, mesg, ier = optimize.leastsq(
      resfunc, p_guess, args=(w, Z), 
      full_output=True) 
     ssr = abs(((model(w, *params) - Z)**2).sum()) 
     print(fmt.format(resfunc.__name__, params, ssr, infodict['nfev'])) 
    print(end='\n') 

R, C, L = 3, 2, 4 
p_guess = 1, 1, 1 
seed = 2013 

model = model1 
np.random.seed(seed) 
w, Z = make_data(R, C, L) 
assert np.allclose(model1(w, R, C, L), model2(w, R, C, L)) 

print('Using model1') 
compare(residuals, MS_residuals, alt_residuals) 

model = model2 
print('Using model2') 
compare(residuals, MS_residuals, alt_residuals) 

ergibt

Using model1 
name   | params        | sum(residuals**2) | ncalls 
------------------------------------------------------------------------------------ 
residuals  | [ 2.86950167 1.94245378 4.04362841] |    9.41 |  89 
MS_residuals | [ 2.85311972 1.94525477 4.04363883] |    9.26 |  29 
alt_residuals | [ 2.85311972 1.94525477 4.04363883] |    9.26 |  29 

Using model2 
name   | params        | sum(residuals**2) | ncalls 
------------------------------------------------------------------------------------ 
residuals  | [ 2.86590332 1.9326829 4.0450271 ] |    7.81 | 483 
MS_residuals | [ 2.85422448 1.94853383 4.04333851] |    9.78 | 754 
alt_residuals | [ 2.85422448 1.94853383 4.04333851] |    9.78 | 754 

So scheint es, die Restfunktion zu verwenden, kann auf der Modellfunktion abhängen. Ich bin zu einem Verlust, um den Unterschied im Ergebnis angesichts der Ähnlichkeit von model1 und model2 zu erklären.

Answer 2

, gibt es keine Notwendigkeit, die verfügbaren Informationen zu reduzieren, indem die Größe im Quadrat in Funktion Residuen, weil leastsq ist egal, ob die Zahlen real oder komplex sind, aber nur, dass sie als 1D-Array ausgedrückt werden, beibehalten die Integrität der funktionalen Beziehung. Hier

ist ein Ersatz Residuen Funktion:

def residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    z1d = np.zeros(Z.size*2, dtype = np.float64) 
    z1d[0:z1d.size:2] = diff.real 
    z1d[1:z1d.size:2] = diff.imag 
    return z1d 

Dies ist die einzige Änderung, die vorgenommen werden müssen. Die Antworten für den Samen 2013 sind: [2.96564781, 1.99929516, 4.00106534].

Die Fehler bezogen auf to unutbu answer's sind deutlich um mehr als eine Größenordnung reduziert.