У меня есть 4 нелинейных уравнения с тремя неизвестными X, Y и Z, для которых я хочу решить. Уравнения имеют вид:
F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ
... где a, b и c - это константы, зависящие от каждого значения F в четырех уравнениях.
Каков наилучший способ решить это?
Есть два способа сделать это.
Итак, насколько я понимаю ваш вопрос, вы знаете F, a, b и c в 4 разных точках, и вы хотите инвертировать параметры модели X, Y и Z. У нас есть 3 неизвестных и 4 наблюдаемых данных точек, поэтому проблема переопределена. Поэтому мы будем решать в наименьшем квадрате смысл.
В этом случае более распространено использование противоположной терминологии, поэтому позвольте перевернуть свое уравнение вокруг. Вместо:
F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ
Пусть пишут:
F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)
Где мы знаем F, X, Y и Z в 4 разных точках (например, F_0, F_1, ... F_i).
Мы просто меняем имена переменных, а не само уравнение. (Это больше для моей легкости мышления, чем что-либо еще.)
Фактически можно линеаризовать это уравнение. Вы можете легко решить для a^2, b^2, a b cos(c) и a b sin(c). Чтобы сделать это немного проще, давайте снова переделаем вещи:
d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)
Теперь уравнение намного проще: F_i = d + e X_i + f Y_i + g Z_i. Легко сделать линейную инверсию наименьших квадратов для d, e, F и g. Тогда мы можем получить a, b и c из:
a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)
Хорошо, напишите это в матричной форме. Мы собираемся перевести 4 наблюдения (код, который мы напишем, будет занимать любое количество наблюдений, но пусть он будет оставаться конкретным в данный момент):
F_i = d + e X_i + f Y_i + g Z_i
В:
|F_0| |1, X_0, Y_0, Z_0| |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2| |1, X_2, Y_2, Z_2| |f|
|F_3| |1, X_3, Y_3, Z_3| |g|
Или: F = G * m (я геофизик, поэтому мы используем g для "зеленых функций" и m для "параметров модели". Обычно мы использовали d для "данных" вместо F.)
В python это будет выглядеть так:
def invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
Вы также можете решить эту проблему, используя scipy.optimize, как предположил @Joe. Наиболее доступной функцией в scipy.optimize является scipy.optimize.curve_fit, которая по умолчанию использует метод Levenberg-Marquardt.
Levenberg-Marquardt - это алгоритм "восхождения на холм" (ну, в данном случае он идет вниз, но этот термин используется в любом случае). В некотором смысле вы делаете первоначальное предположение о параметрах модели (все по умолчанию в scipy.optimize) и следуете по склону observed - predicted в вашем пространстве параметров вниз донизу.
Предостережение: Выбор правильного метода нелинейной инверсии, первоначальная догадка и настройка параметров метода - это очень "темное искусство". Вы только учитесь этому, делая это, и есть много ситуаций, когда вещи не будут работать должным образом. Levenberg-Marquardt - хороший общий метод, если пространство параметров довольно гладкое (это должно быть). Есть много других (включая генетические алгоритмы, нейронные сети и т.д. В дополнение к более распространенным методам, таким как имитированный отжиг), которые лучше в других ситуациях. Я не собираюсь вникать в эту часть здесь.
Существует одна общая информация о том, что некоторые инструменты оптимизации пытаются исправить, что scipy.optimize не пытается обрабатывать. Если ваши параметры модели имеют разные величины (например, a=1, b=1000, c=1e-8), вам нужно будет перемасштабировать вещи так, чтобы они были одинаковыми по величине. Иначе scipy.optimize "алгоритмы подъема холма" (например, LM) не будут точно рассчитать оценку локального градиента и дадут дико неточные результаты. На данный момент я предполагаю, что a, b и c имеют относительно близкие значения. Кроме того, имейте в виду, что по существу все нелинейные методы требуют от вас первоначального предположения и чувствительны к этой догадки. Я оставляю его ниже (просто передайте его как p0 kwarg на curve_fit), потому что по умолчанию a, b, c = 1, 1, 1 является довольно точным предположением для a, b, c = 3, 2, 1.
С учетом предостережений curve_fit ожидает, что будет передана функция, набор точек, в которых были сделаны наблюдения (как один массив ndim x npoints), и наблюдаемые значения.
Итак, если мы напишем такую функцию:
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
Нам нужно обернуть его, чтобы принять несколько разные аргументы, прежде чем передавать его на curve_fit.
В двух словах:
def nonlinear_invert(f, x, y, z):
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
Чтобы дать вам полную реализацию, вот пример, который
import numpy as np
import scipy.optimize as opt
def main():
nobservations = 4
a, b, c = 3.0, 2.0, 1.0
f, x, y, z = generate_data(nobservations, a, b, c)
print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
print linear_invert(f, x, y, z)
print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
print nonlinear_invert(f, x, y, z)
def generate_data(nobservations, a, b, c, noise_level=0.01):
x, y, z = np.random.random((3, nobservations))
noise = noise_level * np.random.normal(0, noise_level, nobservations)
f = func(x, y, z, a, b, c) + noise
return f, x, y, z
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
def linear_invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
def nonlinear_invert(f, x, y, z):
# "curve_fit" expects the function to take a slightly different form...
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
main()
Вероятно, вы хотите использовать scipy нелинейные решатели, они очень легкие: http://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html