Python: установка нелинейной функции на поверхность

У меня есть функция Гаусса ошибок erfc(x), которую мне нужно вписать в мои данные о поверхности. Полное уравнение:

 Z = Z_0 * erfc(x / 2*sqrt(D*t))

Я знаю из данных Z, Z_0, x, t ... единственный параметр, который я ищу, это D. Использование curve_fit подходит для отдельных линий, но мне нужно найти только один постоянный параметр D для всей поверхности. Поверхность выглядит так: surface]. Есть идеи, пожалуйста? Спасибо


person Tripo    schedule 28.04.2018    source источник
comment
Являются ли обе переменные x и t независимыми?   -  person MPA    schedule 28.04.2018
comment
Да, это положение и время.   -  person Tripo    schedule 28.04.2018
comment
Примечание: аргумент erfc должен быть x / (2*sqrt(D*t))?   -  person MPA    schedule 28.04.2018


Ответы (1)


Я создал пример, демонстрирующий многомерность curve_fit. Обратите внимание, что я использовал объект класса для хранения параметра Z0, но есть и другие способы сделать это (см. этот вопрос) .

import numpy as np
from scipy.optimize import curve_fit
from scipy.special import erfc

# Class to contain model and parameters
class fitClass:

    def __init__(self):
        pass

    # Model with unknown parameter D
    def func(self, p, D):
        x, t = p
        Z = self.Z0 * erfc(x / 2*np.sqrt(D*t))
        return Z

# Instantiate class and define parameters
inst = fitClass()
inst.Z0 = 1.0
D = 10.0
Nx = int(1e2)
Nt = int(1e1)

# Independent variables
x = np.linspace(-1.0, 1.0, Nx)
t = np.linspace(1.0, 5.0, Nt)

X, T = np.meshgrid(x, t)

# Merge independent variables
xdata = np.vstack([X.reshape(-1), T.reshape(-1)])

# Synthetic ydata (noisy measurement)
noise = 0.5*(np.random.rand(Nx*Nt)-0.5)
Z = inst.func(xdata, D)
Z_noisy = Z + noise

# Fit model to data
popt, pcov = curve_fit(inst.func, xdata, Z_noisy)
D_fit = popt[0]
print(D_fit)
person MPA    schedule 28.04.2018
comment
Спасибо, это выглядит полезным. Однако мои данные Z представляют собой 2D-массив. Когда вы посмотрите на прилагаемый рисунок. Одна ось — это позиция x, вторая горизонталь — время t, а третья вертикаль — Z. - person Tripo; 28.04.2018
comment
Сюжет в вашем вопросе - это визуализация вашей функции на сетке, но ваши данные - это дискретные измерения? т.е. для каждого значения x есть одно соответствующее t и Z? Возможно, вы можете дать немного больше контекста, в котором можно интерпретировать ваш вопрос. - person MPA; 28.04.2018
comment
Простите за это. Это данные сканирующей зондовой микроскопии, и данные дискретны, поэтому одна горизонтальная ось — это время, скажем (0, 2, 4, ... секунды). Вторая горизонтальная ось — это положение, а вертикальная линия — это Z=f(x,t). В зависимости от этого каждая комбинация x и t имеет определенное значение Z. Если я подойду достаточно хорошо, я смогу найти значение параметра D, которое должно быть постоянным для всей поверхности, а Z_0 является максимальным значением Z. Я сделал это яснее? - person Tripo; 28.04.2018
comment
Порядок данных не должен иметь значения, поэтому xdata может быть представлен двумя одномерными массивами вместо двух структурированных матриц. Я обновил свое решение, в котором X и T представляют данные сетки, которые вы могли бы получить в результате анализа SPM. - person MPA; 28.04.2018