Студенты, лабы и gnuplot: обработка данных

Я преподаю курс экспериментальной оптики в одном из российских университетов. Складывается ощущение, что развитие вычислительной техники и программного обеспечения не оказывает влияния на ожидаемый прогресс в удобстве проведения измерений и качестве представления результатов лабораторных работ. Думаю, причина состоит в том, что офисные приложения для работы с электронными таблицами воспринимаются студентами как единственный подходящий инструмент для решения подобных задач. Чтобы развеять это сверхпопулярное заблуждение, я на простом примере расскажу об использовании пакета gnuplot.

Установка gnuplot

Gnuplot – высокоуровневый язык команд и сценариев, предназначенный для построения графиков математических функций и работы с данными, активно развивающийся с 1986 года. Исходный код gnuplot защищён авторским правом, но распространяется как свободное программное обеспечение и способен работать на Linux, OS/2, MS Windows, OSX, VMS и многих других платформах. Для установки gnuplot на компьютер с MS Windows наиболее удачным решением будет обратиться к официальному сайту gnuplot.info, перейти по ссылке «Download» и загрузить актуальную версию с ресурса sourceforge.net. Размер скачиваемого файла – около 30 Мб, после установки пакет занимает примерно 100 Мб дискового пространства.

Установщик задаст несколько простых вопросов о конфигурации (на английском языке), одно из окошек будет называться «Select Additional Tasks», в нем я рекомендую выбрать тип терминала wxt, а также поставить галочку в последнем пункте «Add application directory to your PATH environment variable». После установки в меню запуска программ появятся два новых пункта, кажется они называются «gnuplot» и «gnuplot console version». Если выбрать второй вариант, появится черное окно с командной строкой, как показано на рисунке:

Если ввести командуpwd, вы увидите директорию, в которой запускается gnuplot. Команда
plot sin(x)откроет графический терминал wxt и нарисует в нем график функции

$\sin(x)$

Если запускать просто gnuplot (не консоль), окно для ввода команд будет белое с черным текстом, при этом в нем будет отображаться верхнее меню с многочисленными командами. В gnuplot есть команда help для быстрого получения справки по любой из команд. Для облегчения работы с gnuplot пользователям MS Windows настоятельно рекомендую обзавестись нормальным текстовым редактором (я не очень в курсе текущей ситуации с Notepad).

Эксперимент

В качестве примера для настоящей заметки выберем простой опыт, относящийся к теме «поляризация света». Пусть у нас есть направленный источник линейно-поляризованного света, луч от которого проходит через линейный поляризатор и попадает на фотодетектор.

Измеряем напряжение на фотодетекторе, которое пропорционально мощности прошедшего через поляризатор излучения. Поляризатор представляет из себя тонкую пластинку, плоскость которой перпендикулярна лучу света. Поляризатор установлен в оправу, позволяющую вращать его вокруг направления луча, изменяя угол между плоскостью поляризации света и оптической осью поляризатора. Угол поворота оси поляризатора считывается со шкалы на оправе.

Теоретическая зависимость измеряемого напряжения

от угла поворота поляризатора

$\theta$

описывается законом Малюса:

$U(\theta) = U_0\cos^2(\theta-\theta_0),$

где

– максимальное значение напряжения,

$\theta_0$

– угол между плоскостью поляризации света и оптической осью поляризатора.

Моделирование

Применим gnuplot для моделирования возможных результатов эксперимента, в котором измеряется зависимость

от угла поворота

$\theta$

. Стоит отметить, что в реальных экспериментах моделирование изучаемого явления широко используется для создания и тестирования средств анализа результатов измерений. Отличие наблюдаемых величин от модельных данных заставляет исследователей проверять и совершенствовать модель, добавляя в нее уже известные науке эффекты. Ситуация, в которой результаты эксперимента невозможно описать известными явлениями называется открытием. Итак, создадим текстовый файл model.gpl со сценарием из gnuplot-команд и комментариев к ним:

В окне gnuplot набираем командуload 'model.gpl'.Подразумевается, что созданный нами файл находится в рабочей директории, в противном случае директорию можно сменить, введя команду cd 'your_path'. После выполнения скрипта мы получим текстовый файл cрезультатами моделированияexp.dat, который выглядит примерно так:

Анализ данных

Начнем обработку данных с построения графика данных моделирования (или измерений). Также используем метод наименьших квадратов для сравнения полученных результатов с теорией. Редуцированное значение взвешенного параметра хи-квадрат определяется так:

$\chi_{\nu}^2 = \frac{\chi^2}{\nu} = \frac{1}{\nu}\sum_i \frac{(Y_i - f(X_i))^2}{\sigma_i^2}.$

$\nu$
называется числом степеней свободы алгоритма поиска минимального значения
$\chi^2$
и определяется как число измерений (точек на графике) минус количество свободных параметров функции
,
– точки на оси x, в которых проведены измерения,
– результат измерения в точке
,
– значение теоретической функции в точке
,
$\sigma_i$
– среднеквадратическое отклонение измеренного значения величины
от ее среднего значения
$\langle Y_i \rangle$
.

В gnuplot отношение

$\chi^2/\nu$

из последней формулы обозначается как WSSR / NDF, что является аббревиатурой от фраз «Weighted Sum of Squared Residuals» и «Number of Degrees of Freedom».

Рассматриваемый здесь пример лабораторной работы является частным случаем обширного класса задач, в которых есть набор данных измерений/моделирования и теоретическая функция, которая «должна» описывать изучаемое явление. Для проверки соответствия между теорией и экспериментом мы будем варьировать свободные параметры функции так, чтобы минимизировать величину

$\chi^2$

. Для этого в gnuplot есть команда fit, изпользующая алгоритм Левенберга — Марквардта. На русский язык фраза «fit function to data» переводится как «подогнать функцию к данным».

Создадим текстовый файлlook.gpl со сценарием для gnuplot.

Символ «\» используется в сценариях для переноса длинных строк и должень быть последним символом в своей строке. В окне gnuplot вводим команду load 'look.gpl', в результате выполнения которой у нас появится файл look.png с построенными графиками данных и теоретической функции:

Итак, мы видим, что «на глазок» данные хорошо совпадают с теорией. Нам даже удалось правильно измерить угол

$\theta_0$

с точностью 0.2 градуса. Однако, мы имеем среднеквадратическое отклонение теории от моделирования

$\sigma \simeq 9$

милливольт. Что бы это могло означать?

Если посмотреть на использование команды fit (строка №11 в файле look.gpl), можно заметить, что мы не передаём процедуре подгонки никаких данных об ошибках измерений, т. е. алгоритм ничего «не знает» о величинах

$\sigma_i$

из последней формулы. Что же он (алгоритм) делает в таком случае? А вот что: все

$\sigma_i$

считаются равными безразмерной единице и мы получаем невзвешенное значение редуцированного параметра хи-квадрат, квадратный корень из которого является размерной величиной, описывающей среднеквадратическое отклонение функции от экспериментальных точек – те самые 8.88 милливольт, приведённые на графике.

Если в процессе измерений текущее наблюдаемое значение напряжения в каждой точке ведет себя достаточно стабильно, можно предположить, что отличие теории и эксперимента может быть связано с неточностью отсчета угла по шкале оправы с вращающимся поляризатором. Чтобы преобразовать ошибки отсчета угла в эквивалентные ошибки измерения напряжения, продифференцируем функцию

$U(\theta)$

и назовем получившуюся функцию

$V(\theta)$

$V(\theta) \equiv \frac{\partial U(\theta)}{\partial \theta} = U_0 \sin(2(\theta_0-\theta)) .$

В нашем моделировании (см. файлmodel.gplв начале этого текста) мы «учли» погрешность измерения угла, добавив его к закону Малюса как случайную величину с нормальным распределением и среднеквадратическим разбросом в 1 градус. Чтобы учесть погрешности измерений, создадим усовершенствованный скрипт для обработки данных. Заодно, нарисуем графики в полярных координатах (удобных для визуализации поляризационных явлений), а также сохраним картинку в формате pdf.

Следующее изображение - скриншот получившегося файла closelook.pdf.

В результате применения алгоритма подгонки мы получили редуцированное значение взвешенного параметра хи-квадрат

$\chi^2/\text{NDF}=33.9/33$

. Используя распределение хи-квадрат с соответствующим числом степеней свободы gnuplot находит значение вероятности (p-value) получения наблюдаемого результата в предположении, что исходная теоретическая гипотеза верна. В нашем случае эта вероятность довольно высока - 42%.

Заключение

Изучение пакета gnuplot, как и многих других инструментов, целесообразно начинать на примере решения какой-нибудь простой задачи. Существование готовых сценариев позволяет в дальнейшем легко адаптировать и видоизменять их под свои нужды. Открытость и мультиплатформенность программы и языка сценариев обеспечивает легкий обмен знаниями, данными, навыками и результатами работы.