Не могу применить сетку напрямую для создания Surfaceplot

Question

Рассмотрим следующий пример:

from numpy import array, exp, linspace, meshgrid
import matplotlib.pyplot as plt

x = array([1, 2, 3, 4, 5, 6])
y = array([10, 5.49, 0.89, -0.14, -1.07, 0.84])
a1, a2 = 1, 4

def model(x, a1, a2):
    return a1 * exp(-a2 * x)

def residuen(x, y, a1, a2):
    return modell(x, a1, a2) - y

def S(x, y, a1, a2): 
    return 0.5 * sum(residuen(x, y, a1, a2) ** 2)

Входные значения x и y являются фиксированными. Я хочу вычислить значения функции S для различных значений a1, a2 и отобразить результат в виде графика 3D-поверхности. В учебнике matplotlib я нашел пример, в котором используется функция meshgrid. Я пытался применить эту функцию для моей проблемы:

a1 = linspace(-100, 100, 1000)
a2 = linspace(-100, 100, 1000)
A1, A2 = meshgrid(a1, a2)
Z = S(x, y, A1, A2) 

Когда я запускаю этот код, я получаю ошибку

ValueError: operands could not be broadcast together with shapes (1000,1000) (6,)

Ясно, почему приведенный выше код не работает, но Я не знаю, как исправить / переопределить функцию S, чтобы я мог применить функцию meshgrid, чтобы получить желаемые значения Z для графика поверхности. У кого-нибудь есть идеи, как решить эту проблему?

Answer 1

Как показывает ваша ошибка, ваша проблема не имеет ничего общего с matplotlib. Вы пытаетесь вычислить остаток между вашей моделью и вашими 6 точками данных в каждом из 1e6 местоположений. Поскольку остаток представляет собой сумму, вы можете представить массив 1000x1000x6 в некоторой точке и суммировать его по последнему измерению. Вы также можете создать массив 6x1000x1000 и суммировать по первому измерению, но я покажу первый случай.

Метод, который вы хотите использовать, называется broadcasting . Это означает выравнивание форм, начиная с последних размеров, и заполнение размеров единиц до соответствия. Самое простое изменение, которое необходимо внести для введения вещания, - это model:

a1[..., None] * exp(-a2[..., None] * x)

. Многоточие (...) в [..., None] означает «взять все существующие измерения». В этом случае это сокращение от 1011 *, но столько раз, сколько необходимо для конкретного массива. None эквивалентно np.newaxis и означает «добавить единицу измерения». В итоге, индексное выражение - это идиоматический c способ добавить измерение к вашей фигуре. Другим способом было бы сделать что-то вроде a1.reshape(a1.shape + (1,)).

Теперь у вас есть следующие операции:

1. Умножение формы (1000, 1000, 1) (расширение a2) на форму (6,) (x) -> Форма результата: (1000, 1000, 6). Последнее измерение a2 - , транслируемое . x предоставляется два неявных начальных размера размером 1, которые также транслируются .

   Здесь произошла ваша ошибка. Вы не можете умело умножать массивы формы (1000, 1000) и (6,) поэлементно, но вы можете с новыми единицами измерения.

2. Экспонента формы (1000, 1000, 6) массив сохраняет форму.
3. Умножьте результат на форму (1000, 1000, 1) (расширен a1). Последнее измерение - , транслируемое .

Теперь остаток в residuen в порядке: modell(x, a1, a2) - y находит разницу между массивом (1000, 1000, 6) и массивом (6,) , Окончательные размеры идеально выстраиваются, а вещание позаботится обо всем остальном.

Единственное другое отличие состоит в том, что вы должны указать sum, на какую ось работать, поскольку ваши массивы не длиннее плоское:

0.5 * sum(residuen(x, y, a1, a2)**2, axis=1)

Это не совсем правильное вычисление RMS. Я бы ожидал увидеть

sum(residuen(x, y, a1, a2)**2, axis=-1)**0.5

Более простой способ найти среднеквадратичное значение массива - использовать np.linalg.norm:

np.linalg.norm(residuen(x, y, a1, a2), axis=-1)

TL; DR

Ваша программа должна использовать трансляцию:

x = array([1, 2, 3, 4, 5, 6])
y = array([10, 5.49, 0.89, -0.14, -1.07, 0.84])
a1 = np.linspace(-100, 100, 1000)[None, :]
a2 = np.linspace(-100, 100, 1000)[:, None]

def model(x, a1, a2):
    return a1[..., None] * exp(-a2[..., None] * x)

def residuen(x, y, a1, a2):
    return model(x, a1, a2) - y

def S(x, y, a1, a2): 
    return np.linalg.norm(residuen(x, y, a1, a2), axis=-1)

Z = S(x, y, A1, A2)