Python: подгонка кривой к замаскированным данным с помощью scipy curve_fit

Я пытаюсь написать скрипт с python/numpy/scipy для обработки данных, подгонки и построения графиков измерений магнитосопротивления, зависящих от угла. Я новичок в Python, получил код фрейма от своего научного руководителя и сумел добавить в фрейм несколько сотен строк кода. Через некоторое время я заметил, что в некоторых измерениях было несколько грубых ошибок, и, поскольку скрипт должен выполнять все манипуляции автоматически, я попытался замаскировать эти точки и подогнать кривую к незамаскированным точкам (кривая представляет собой квадрат синуса, наложенный на линейную функцию, так что numpy.ma.polyfit на самом деле не вариант). Однако после маскирования координат x и y проблемных точек подгонка по-прежнему будет учитывать их, даже если они не будут отображаться на графике. Пример упрощен, но происходит то же самое;

import numpy.ma as ma
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit



def Funk(x, k, y0):
 return k*x + y0   

fig,ax= plt.subplots()

x=ma.masked_array([1,2,3,4,5,6,7,8,9,10],mask=[0,0,0,0,0,0,1,1,1,1])
y=ma.masked_array([1,2,3,4,5,30,35,40,45,50], mask=[0,0,0,0,0,1,1,1,1,1])


fitParamsFunk, fitCovariancesFunk = curve_fit(Funk, x, y)

ax.plot(x, Funk(x, fitParamsFunk[0], fitParamsFunk[1]))
ax.errorbar(x, y, yerr = None, ms=3, fmt='-o')
plt.show()

Вторая половина точек маскируется и не отображается на графике, но все же учитывается.< /а>

Во время написания поста я понял, что я могу сделать это:

def Funk(x, k, y0):
    return k*x + y0   

fig,ax= plt.subplots()

x=np.array([1,2,3,4,5,6,7,8,9,10])
y=np.array([1,2,3,4,5,30,35,40,45,50])
mask=np.array([0,0,0,0,0,1,1,1,1,1])

fitParamsFunk, fitCovariancesFunk = curve_fit(Funk, x[mask], y[mask])

ax.plot(x, Funk(x, fitParamsFunk[0], fitParamsFunk[1]))
ax.errorbar(x, y, yerr = None, ms=3, fmt='-o')
plt.show()

Чего я на самом деле хотел

Я предполагаю, что scipy curve_fit не предназначен для работы с замаскированными массивами, но я все же хотел бы знать, есть ли какой-либо обходной путь для этого (мне нужно работать с замаскированными массивами, потому что количество точек данных> 10e6, но я' m только строит 100 сразу, поэтому мне нужно будет взять маску части массива, которую я хочу построить, и назначить ее другому массиву, копируя значения массива в другой или устанавливая исходную маску на False) ? Спасибо за любые предложения


person krdzi    schedule 26.03.2020    source источник


Ответы (3)


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

x = ma.masked_array([1,2,3,4,5,6,7,8,9,10], mask=[0,0,0,0,0,1,1,1,1,1])  # changed mask
y = ma.masked_array([1,2,3,4,5,30,35,40,45,50], mask=[0,0,0,0,0,1,1,1,1,1])

fitParamsFunk, fitCovariancesFunk = curve_fit(Funk, x[~x.mask], y[~y.mask])

PS: обратите внимание, что оба массива должны иметь одинаковое количество допустимых записей.

person joni    schedule 27.03.2020

Использование маски в численном исчислении эквивалентно использованию ступенчатой ​​функции Хевисайда в аналитическом исчислении. Например, это становится очень простым при применении кусочно-линейной регрессии:

введите здесь описание изображения

Это несколько примеров кусочно-линейной регрессии в документе: https://fr.scribd.com/document/380941024/Regression-par-morceaux-Piecewise-Regression-pdf

Используя метод, показанный в этой статье, очень простое исчисление ниже приводит к ожидаемой форме результата:

введите здесь описание изображения

Примечание: В случае большого количества точек, если в области перехода имеется несколько точек с немного отличающимися абсциссами, более точным будет применить случай, рассмотренный на страницах 29-31 упомянутой выше статьи.

person JJacquelin    schedule 28.03.2020

Я думаю, что вы хотите определить маску, в которой перечислены индексы «хороших точек данных», а затем использовать ее в качестве точек для подгонки (и/или для построения).

Как ведущий автор lmfit, я бы порекомендовал использовать эту библиотеку для подгонки кривых: она имеет много полезных функций по сравнению с curve_fit. При этом ваш пример может выглядеть так:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

def Funk(x, k, y0, good_points=None):  # note: add keyword argument
    f = k*x + y0
    if good_points is not None:
        f = f[good_points]       # apply mask of good data points
    return f

x = np.array([1,2,3,4,5, 6,7,8.,9,10])
y = np.array([1,2,3,4,5,30,35.,40,45,50]) 
y += np.random.normal(size=len(x), scale=0.19) # add some noise to make it fun

# make an array of the indices of the "good data points"
# does not need to be contiguous.
good_points=np.array([0,1,2,3,4])

# turn your model function Funk into an lmfit Model
mymodel = Model(Funk)

# create parameters, giving initial values. Note that parameters are
# named using the names of your function's argument and that keyword 
# arguments with non-numeric defaults like 'good points' are seen to
#  *not* be parameters. Like the independent variable `x`, you'll 
# need to pass that in when you do the fit.
# also: parameters can be fixed, or given `min` and `max` attributes

params = mymodel.make_params(k=1.4,  y0=0.2)
params['k'].min = 0

# do the fit to the 'good data', passing in the parameters, the 
# independent variable `x` and the `good_points` mask.
result  = mymodel.fit(y[good_points], params, x=x, good_points=good_points)

# print out a report of best fit values, uncertainties, correlations, etc.
print(result.fit_report())

# plot the results, again using the good_points array as needed.
plt.plot(x, y, 'o', label='all data')
plt.plot(x[good_points], result.best_fit[good_points], label='fit to good data')
plt.legend()
plt.show()

Это распечатает

[[Model]]
    Model(Funk)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 7
    # data points      = 5
    # variables        = 2
    chi-square         = 0.02302999
    reduced chi-square = 0.00767666
    Akaike info crit   = -22.9019787
    Bayesian info crit = -23.6831029
[[Variables]]
    k:   1.02460577 +/- 0.02770680 (2.70%) (init = 1.4)
    y0: -0.04135096 +/- 0.09189305 (222.23%) (init = 0.2)
[[Correlations]] (unreported correlations are < 0.100)
    C(k, y0) = -0.905

и создайте график введите здесь описание изображения

Надеюсь, это поможет вам начать.

person M Newville    schedule 01.04.2020