Комплексные числа минимизации наименьших квадратов

Question

Я использую свой Matlab, но мое видение заключается в том, чтобы в конечном итоге переключиться на выполнение всего моего анализа в Python, поскольку это настоящий язык программирования и несколько других причин.

Недавняя проблема, которую я пытался решить, заключается в минимизации сложных данных методом наименьших квадратов. Я инженер, и мы довольно часто имеем дело со сложным импедансом, и я пытаюсь использовать подгонку кривой для подгонки простой модели схемы к измеренным данным.

Уравнение импеданса выглядит следующим образом:

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

Затем я пытаюсь найти значения R, C и L так, чтобы была найдена кривая наименьших квадратов.

Я пытался использовать пакет оптимизации, такой как optimize.curve_fit или optimize.leastsq, но они не работают со сложными числами.

Затем я попытался заставить мою функцию невязки возвращать величину сложных данных, но это тоже не сработало.

Answer 1

Обращаясь к ответам unutbu, нет необходимости уменьшать доступную информацию, принимая квадрат в квадрате от остатков функций, потому что leastsq не имеет значения, являются ли числа действительными или сложными, а только то, что они выражены в виде одномерного массива, сохраняя целостность функциональных отношений.

Вот функция замены остатков:

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

Это единственное изменение, которое нужно сделать. Ответы на семена 2013 года: [2.96564781, 1.99929516, 4.00106534].

Ошибки относительно ответов unutbu уменьшены значительно более чем на порядок.

Answer 2

В конечном счете, цель состоит в том, чтобы уменьшить абсолютную величину суммы квадратов различий между моделью и наблюдаемой Z:

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

Мой оригинальный ответ предложил применить leastsq к residuals функция, которая возвращает скаляр, представляющий сумму квадратов действительных и мнимых разностей:

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

Майк Зульцер предложил использовать остаточную функцию, которая возвращает вектор с плавающей точкой.

Вот сравнение результата с использованием этих остаточных функций:

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://stackru.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)

доходность

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

Таким образом, представляется, что остаточная функция для использования может зависеть от функции модели. Я затрудняюсь объяснить разницу в результате, учитывая сходство model1 а такжеmodel2,