Numeri minimi di minimizzazione dei quadrati minimi

Question

Utilizzo il mio Matlab, ma alla fine ho deciso di passare a eseguire tutte le mie analisi in Python poiché si tratta di un vero linguaggio di programmazione e di altri motivi.

Il problema recente che ho cercato di affrontare è la minimizzazione dei minimi dei dati complessi. Sono un ingegnere e ci occupiamo di impedenze complesse abbastanza spesso, e sto cercando di utilizzare il fitting della curva per adattare un semplice modello di circuito ai dati misurati.

L'equazione impedenza è la seguente:

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

Sto quindi cercando di trovare i valori di R, C e L tali da trovare la curva dei minimi quadrati.

Ho provato a utilizzare il pacchetto di ottimizzazione come optimize.curve_fit o optimize.leastsq, ma non funzionano con numeri complessi.

Ho quindi provato a rendere la mia funzione residua restituire la grandezza dei dati complessi, ma anche questo non ha funzionato.

Answer 1

In ultima analisi, l'obiettivo è quello di ridurre il valore assoluto della somma dei quadrati differenze tra il modello e l'osservato Z:

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

mio original answer ha suggerito l'applicazione di leastsq a una funzione residuals che restituisce uno scalare che rappresenta la somma dei quadrati delle differenze reali e immaginarie:

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 utilizzando una funzione residua che restituisce un vettore di float.

Ecco un confronto tra il risultato usando queste funzioni residue:

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) 

cede

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 

Così sembra che la funzione residua da utilizzare può dipendere la funzione di modello. Sono in perdita per spiegare la differenza di risultato data la somiglianza di model1 e model2.

Answer 2

Con riferimento a to unutbu answer's, non è necessario ridurre le informazioni disponibili prendendo la magnitudine al quadrato nei residui di funzione poiché leastsq non si cura se i numeri sono reali o complessi, ma solo che sono espressi come un array 1D, preservando l'integrità della relazione funzionale.

Ecco una funzione residui pezzi:

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 

Questa è l'unica modifica che essere apportate. Le risposte per il seme 2013 sono: [2.96564781, 1.99929516, 4.00106534].

Gli errori relativi a to unutbu answer's vengono ridotti di molto più di un ordine di grandezza.