Python curve_fit с измеренными точками данных

Я измерил точки данных, я хочу сопоставить формулу для определения двух объектов. Однако я получаю сообщение об ошибке:

TypeError: ufunc 'bitwise_xor' не поддерживается для входных типов, и входные данные не могут быть безопасно приведены к каким-либо поддерживаемым типам в соответствии с правилом приведения ''safe''

создается следующим кодом Python (я использую версию 3):

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

def func(T, fun, Tc):
    return fun*np.sqrt(np.cos(np.pi/2*(T/Tc)^2))

xdata=(4.61, 4.89, 4.92, 4.95, 5.06, 5.10, 5.21, 5.38, 5.41, 5.57, 5.80, 6.14, 6.61, 7.27, 7.66, 7.90, 8.91, 8.29, 8.23, 7.30, 7.86,
       8.30, 8.89, 8.99, 9.24, 9.35, 9.50, 8.77, 8.27, 8.37, 7.72, 7.57, 7.99, 8.13) # these are temperature values <-> T

ydata=(2.85, 2.84, 2.83, 2.825, 2.82, 2.81, 2.80, 2.765, 2.76, 2.74, 2.695, 2.62, 2.50, 2.265, 2.105, 1.975, 1.23, 1.75, 1.81, 2.26,
       2.005, 1.75, 1.31, 1.14, 1.015, 1.045, 1.06, 1.40, 1.75, 1.69, 2.075, 2.15, 1.93, 1.855) # these are energy values <-> func

popt, pcov = curve_fit(func, xdata, ydata)
popt #display these optimized values

Вот вышеприведенная ошибка!!!

Я видел способ сделать это, если у вас есть фиксированная формула и вы добавляете немного шума, но я измерил точки данных (они должны быть дискретными).

Благодарю вас! Карстен


python data-analysis measurement curve-fitting
person Carsten    schedule 25.04.2018    source источник
comment
Оператор мощности ** не ^... ^ означает побитовое исключающее или   -  person DavidG    schedule 25.04.2018
comment
@DavidG, это хороший улов, молодец.   -  person James Phillips    schedule 25.04.2018
comment
Когда я использую **, я получаю другую ошибку после popt: C:\Users\B\Anaconda3\lib\site-packages\ipykernel_main_.py:2: RuntimeWarning: недопустимое значение обнаружено в sqrt из ipykernel импортировать kernelapp как приложение И позже после той же ячейки: RuntimeError: Оптимальные параметры не найдены: количество вызовов функции достигло maxfev = 600.   -  person Carsten    schedule 25.04.2018
comment
Это другая проблема. Вам нужно задать другой вопрос (ответ на который вам, вероятно, нужно предоставить начальное предположение curve_fit с использованием аргумента p0=)   -  person DavidG    schedule 25.04.2018
comment
Так что я не должен редактировать вопрос, но оставить его и задать новый вопрос? (я новичок в stackoverflow)   -  person Carsten    schedule 25.04.2018
comment
И вы думаете, что ядру требуется слишком много времени, чтобы вычислить это соответствие, чтобы оно остановилось после 600 попыток, поэтому я должен дать python некоторый начальный раунд правильных значений, которые он может оптимизировать?   -  person Carsten    schedule 25.04.2018
comment
Пожалуйста, попробуйте эти начальные параметры Tc = 9,7 и fun = 2,9.   -  person James Phillips    schedule 25.04.2018

Ответы (1)


arrow_upward
0
arrow_downward

Я думаю, что большая часть ошибки, которую вы видите, заключается в том, что вы не указали начальные значения для подходящих переменных fun и Tc. К сожалению, curve_fit() scipy позволяет это и молча присваивает значения равным 1.0, что поощряет плохую практику и является очень плохой «функцией». Не используйте его.

Разрешите мне рекомендовать lmfit (https://lmfit.github.io/lmfit-py) который обеспечивает интерфейс более высокого уровня для подбора кривой, который проще в использовании, лучше предотвращает плохое поведение и имеет множество полезных функций, недоступных с curve_fit().

С lmfit ваша проблема соответствия будет выглядеть так (для ясности я изменил несколько имен переменных):

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

def func(t, scale, tc):
    return scale*np.sqrt(np.cos(np.pi/2*(t/tc)**2))

tdata = np.array([4.61, 4.89, 4.92, 4.95, 5.06, 5.10, 5.21, 5.38, 5.41, 5.57, 5.80, 6.14, 6.61, 7.27, 7.66, 7.90, 8.91, 8.29, 8.23, 7.30, 7.86,
                  8.30, 8.89, 8.99, 9.24, 9.35, 9.50, 8.77, 8.27, 8.37, 7.72, 7.57, 7.99, 8.13]) # these are temperature values <-> T

energy = np.array([2.85, 2.84, 2.83, 2.825, 2.82, 2.81, 2.80, 2.765, 2.76, 2.74, 2.695, 2.62, 2.50, 2.265, 2.105, 1.975, 1.23, 1.75, 1.81, 2.26,
                   2.005, 1.75, 1.31, 1.14, 1.015, 1.045, 1.06, 1.40, 1.75, 1.69, 2.075, 2.15, 1.93, 1.855]) # these are energy values <-> func

fmodel = Model(func)

params = fmodel.make_params(scale=5, tc=10)

result = fmodel.fit(energy, params, t=tdata)
print(result.fit_report())

plt.plot(tdata, energy, 'o', label='data')
plt.plot(tdata, result.best_fit, '+', label='fit')

plt.legend()
plt.show()

который распечатает отчет

[[Model]]
    Model(func)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 24
    # data points      = 34
    # variables        = 2
    chi-square         = 0.34988407
    reduced chi-square = 0.01093388
    Akaike info crit   = -151.601474
    Bayesian info crit = -148.548753
[[Variables]]
    scale:  2.87776739 +/- 0.02737439 (0.95%) (init = 5)
    tc:     9.68051725 +/- 0.03597889 (0.37%) (init = 10)
[[Correlations]] (unreported correlations are < 0.100)
    C(scale, tc) = -0.506

и создайте график, подобный

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

person M Newville    schedule 25.04.2018
comment
Как вы получили масштаб = 5 (ширина запрещенной зоны сверхпроводника при T=0K) и tc=10 (критическая температура сверхпроводника - ниобия)? Вы только что попробовали некоторые значения и проверили, что дало результат? (На самом деле эти значения исходят из физической оценки?) - person Carsten; 25.04.2018
comment
Физические значения верны по сравнению с литературой (приблизительно), и я смог воспроизвести их после установки «lmfit». Я не знаю, как это обычно видят в сообществе python, но установка новых пакетов всегда сопряжена с риском того, что где-то есть ошибка, не так ли? - Так что я отмечу как правильный. Тем не менее, ответ без нового пакета (или других для тех, кто предпочитает другие) по-прежнему очень приветствуется - большое спасибо, мистер Ньювиэль, за вашу помощь. - person Carsten; 25.04.2018
comment
Чтобы угадать начальные значения, я построил данные и увидел: а) максимальное значение энергии (y) ~ 2,5 и б) максимальное значение t (x) ~ 10. Это были бы лучшие предположения (!), но я использовал их, чтобы предположить что scale, вероятно, было между 1 и 10 (то есть 5), а tc, вероятно, было между 1 и 100 (то есть 10). Визуализация данных может быть чрезвычайно полезной! При правильных начальных значениях вы должны получить одинаковые результаты с curve_fit или leastsq. lmfit использует scipy.optimize.leastsq() по умолчанию. Любое программное обеспечение имеет риск ошибок. В данном случае было продемонстрировано, что curve_fit хуже, чем lmfit ;). - person M Newville; 25.04.2018