В своей предыдущей заметке на тему обработки данных лабораторных работ я написал об использовании пакета gnuplot – простого и мощного инструмента для решения подобных задач и графического представления результатов. Однако довольно распространённым является мнение, что студенты, которым я советовал использовать gnuplot, вероятно, изучают программирование и способы визуализации данных, и что для них более естественным и полезным будет практическое применение уже полученных навыков в этой сфере. В этом коротком тексте мы рассмотрим применение python с использованием библиотек scipy для обработки данных и matplotlib для представления результатов.
Возьмём тот же пример, что и в предыдущем тексте (там нужно посмотреть формулировку задачи). С моделированием данных python, разумеется, легко справляется: определяем переменные, теоретическую функцию, записываем результаты в текстовый файл.
from random import normalvariate
from math import pi, cos
GtoR, Uo, To, D, B = pi/180., 0.7, 56.0, 1.0, 0.02
def U(T):
return Uo * cos(GtoR*(T + normalvariate(0.0, D) - To))**2 + B
with open('pyexp.dat', 'w') as fp:
for T in range(0,360,10):
fp.write('{:5.1f} {:5.3f}\n'.format(T, U(T)))Ограничимся минимально достаточным набором импортируемых в python модулей. В частности, не будем использовать пакеты csv или pandas для записи/чтения табличных данных. Также воздержимся от прямого импорта numpy, несмотря на то, что numpy и так входит в библиотеку matplotlib. Получившийся в итоге python-сценарий для обработки данных выглядит так:
from math import pi, cos, sin # для вычислений
from scipy.optimize import curve_fit # алгоритм подгонки
from scipy.special import gammainc # чтобы определить p-value
from matplotlib import pyplot as plt # для рисования графиков
GtoR = pi/180.
angle, experiment = [], []
with open('pyexp.dat','r') as fp: data = fp.readlines()
for row in data:
a, u = row.strip().split()
angle.append(GtoR*float(a))
experiment.append(float(u))
def U(X, Uo, To, B): # Теоретическая функция
Y = []
for x in X: Y.append(Uo * cos(x - To)**2 + B )
return Y
def V(X, Uo, To, D): # Её производная
Y = []
for x in X: Y.append(abs(D*(Uo * sin(2*(To - x)))))
return Y
# Подгонка без учёта ошибок измерений
popt, pcov = curve_fit(U, angle, experiment, p0 = [0.7, 1.0, 0.01], bounds=(0.0, [1., 180., 0.5]))
# Переводим ошибки измерений углов (1 градус) в ошибки напряжений
errors = V(angle, popt[0], popt[1], 1.0*GtoR)
# Подгонка с учётом ошибок
popt, pcov = curve_fit(U, angle, experiment, p0 = popt, sigma=errors, absolute_sigma=True)
# "Вытаскиваем" результаты подгонки
theory = U(angle, *popt)
Chi2, NDF = 0.0, len(angle)-len(popt)
for p in range(len(angle)): Chi2 += ((theory[p]-experiment[p])/errors[p])**2
pvalue = 1 - gammainc(0.5*NDF, 0.5*Chi2) # gammainc - regularized upper incomplete gamma function
To, dTo = popt[1]/GtoR, pcov[1][1]**0.5/GtoR
tlabel = 'теория: '
tlabel += r'$\chi^2/NDF$ = {:.1f}/{:2d}'.format(Chi2,NDF)
tlabel += '\nP-value = {:.2f}%'.format(100*pvalue)
theta = [GtoR*i for i in range(0,360)] # чтобы рисовать функцию без изломов
plt.rc('grid', color='#316931', linewidth=0.5, linestyle='-.')
fig = plt.figure(figsize=(8, 8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection='polar', facecolor='#fffff0')
ax.set_title('Поворот оси поляроида: ' + r'$\theta_0={:.2f}^\circ\pm{:.2f}^\circ$'.format(To, dTo)+'\n')
ax.plot(angle, experiment, 'ro', label='эксперимент')
ax.plot(theta, U(theta, *popt), 'b-', label = tlabel)
ax.legend()
plt.show()В результате выполнения программы matplotlib нарисует нам график результатов моделирования (или измерений) и теоретической функции в полярных координатах:
Выводы
Итак, для решения «старой» задачи мы выбрали другой «набор инструментов» –
«python + scipy + matplotlib» вместо gnuplot. И да, решение получено. Преимущества нового метода отчасти сформулированы во введении, кроме того, использование python расширяет границы «дозволенного» почти до бесконечности. Недостатки, по сравнению с gnuplot-версией, на мой взгляд, таковы:
Неподготовленному человеку понять смысл сценария труднее, чем это было в случае с gnuplot.
Библиотека scipy предоставляет необходимые алгоритмы для подгонки теоретической функции к экспериментальным данным. Однако, для доступа к результатам пришлось повозиться: ошибки определения параметров подгонки извлечь из ковариационной матрицы,
– посчитать вручную, для определения p-value использовать регуляризованную гамма-функцию. Библиотека matplotlib весьма лаконична, синтаксис – на любителя, и ещё – не смогла изобразить ошибки (усы) данных на графике в полярных координатах. Последнее обстоятельство, конечно же, можно преодолеть написанием дополнительного кода.
Библиотека NumPy
Благодаря замечательным рецептам, полученным в комментарии @masai к вышеприведенному тексту, посмотрим, как изменится код если всё-таки использовать возможности библиотеки numpy.
import numpy as np # numeric python
from matplotlib import pyplot as plt # для рисования графиков
from scipy.optimize import curve_fit # алгоритм подгонки
from scipy.special import gammainc # чтобы определить p-value
data = np.loadtxt("pyexp.dat")
angle = np.radians(data[:, 0]) # Первый столбик, переводим в радианы
experiment = data[:, 1] # Второй столбик
# Теоретическая функция
def U(X, Uo, To, B): return Uo * np.cos(X - To)**2 + B
# Её производная
def V(X, Uo, To, D): return np.abs(D*(Uo * np.sin(2*(To - X))))
# Приблизительные значения параметров "на глазок"
[Uo, To, B] = 0.7, 1.0, 0.01
# Подгонка без учёта ошибок измерений
[Uo, To, B], pcov = curve_fit(U, angle, experiment, p0 = [Uo, To, B])
# Переводим ошибки измерений углов (1 градус) в ошибки напряжений
errors = V(angle, Uo, To, np.radians(1.0))
# Подгонка с учётом ошибок
[Uo, To, B], pcov = curve_fit(U, angle, experiment, p0 = [Uo, To, B], sigma=errors, absolute_sigma=True)
# "Вытаскиваем" результаты подгонки
theory = U(angle, Uo, To, B)
Chi2,NDF = np.sum(((theory-experiment)/errors)**2), len(angle)-len([Uo, To, B])
pvalue = 1 - gammainc(0.5*NDF, 0.5*Chi2) # gammainc - regularized upper incomplete gamma function
tlabel = f'теория: $\chi^2/NDF$ = {Chi2:.1f}/{NDF:2d}\nP-value = {100*pvalue:.2f}%'
theta = np.linspace(0, 2 * np.pi, 1000) # чтобы нарисовать функцию без изломов
gTo, dTo = np.degrees(To), np.degrees(pcov[1][1]**0.5)
# Рисуем график в полярных координатах
plt.rc('grid', color='#316931', linewidth=0.5, linestyle='-.')
fig = plt.figure(figsize=(8, 8))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection='polar', facecolor='#fffff0')
ax.set_title(f'Поворот оси поляроида: $\\theta_0={gTo:.2f}^\circ\pm{dTo:.2f}^\circ$\n')
ax.plot(angle, experiment, 'ro', label='эксперимент')
ax.plot(theta, U(theta, Uo, To, B), 'b-', label = tlabel)
ax.legend()
plt.show()Скрипт стал выглядеть лаконичнее, спору нет. В своей практике я что-то давно забросил numpy, так как в основном пользуюсь библиотеками ROOT, а смешивать разные идеологии внутри кода не очень удобно...
