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

1. Неподготовленному человеку понять смысл сценария труднее, чем это было в случае с gnuplot.

2. Библиотека scipy предоставляет необходимые алгоритмы для подгонки теоретической функции к экспериментальным данным. Однако, для доступа к результатам пришлось повозиться: ошибки определения параметров подгонки извлечь из ковариационной матрицы, – посчитать вручную, для определения p-value использовать регуляризованную гамма-функцию.

3. Библиотека 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, а смешивать разные идеологии внутри кода не очень удобно...