В своей предыдущей заметке на тему обработки данных лабораторных работ я написал об использовании пакета 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, а смешивать разные идеологии внутри кода не очень удобно...
