Уточнение данных полиномом n-степени в Python (scipy, polyfitw)

Задача - произвести подгонку данных полиномом n-степени. Известно, что каждое, экспериментально измеренное значение, представляет собой два числа: собственно значение и погрешность измерения. Таким образом, желательно произвести подгонку с учётом весовых коэффициентов каждой точки. Ведь может оказаться так, что в эксперименте есть точки, измеренные с меньшей точностью и степенью достоверности, но которые не желательно просто так выбрасывать, но и нельзя учитывать наравне с другими, более достоверными.

В качестве инструмента выбираем Python, благодаря расширяемости этого языка.

Модули в Python

Scipy

В Python для численных операций имеются модули SciPy и NumPy. Однако, несмотря на всю мощь этих модулей, я не нашёл в них возможности производить уточнение именно с учетом весовых коэффициентов. В scipy содержится функция polyfit, которая принимает только значения аргумента, функции и степени полинома, но не учитывает погрешности.

Модуль CARSMath

На сайте The Consortium for Advanced Radiation Sources есть разработанный Марком Риверсом модуль CARSMath. Он содержит функцию polyfitw, которая помимо аргумента, функции и степени полинома принимает также и весовые коэффициенты (именно весовые коэффициенты, а не погрешность, то есть 1/погрешность).

Установка модуля CARSMath

Скачайте модуль отсюда. Бросьте тарбол, куда-нибудь, например, в ~/python/python_epics.

Распакуйте:

tar xvf python_epics.tar

Нас интересует файл CARSMath.py.

Теперь, чтобы использовать этот модуль, надо его бросать в рабочую директорию с вашим скриптом, что очень неудобно. Поэтому надо воспользоваться переменной окружения $PYTHONPATH. Для этого отредактируйте файл ~/.bashrc. Необходимо добавить следующую строчку:

export PYTHONPATH="$HOME/distrib/python/lib"

Естественно, директория должна быть указана та, где лежит скачанный модуль CARSMath (я его перекинул в ~/distrib/python/lib, там же у меня лежат и другие модули, не распространяемые с дистрибутивом).

После этого инициализируем оболочку заново:

source ~/.bashrc

Если вы запускаете иксы не через startx, а через менеджер входа в систему, то вам надо узнать каким образом в этом случае инициализировать переменные окружения. Ключевое слово .XClients, если я ничего не путаю.

Теперь можно, находясь в любой директории, импортировать модуль CARSMath стандартным образом:

import CARSMath

Сравнение двух способов подгонки

Сравним два способа подгонки экспериментальных данных полиномом n-й степени:

• с помощью polyfit из scipy
• с помощью polyfitw из CARSMath

Для этого создадим массив точек, сделаем несколько «выбросов» с большими погрешностями и аппроксимируем полиномом (в одном случае без учёта погрешностей, а во втором - с учётом).

Ниже проведён скрипт, реализующий эти действия.

#!/usr/bin/env python
#-*-coding: utf-8 -*-

import scipy
from numpy import arange
import pylab as pl
from CARSMath import polyfitw

x = arange(0, 2, 0.01)
y = 0.5 * x**3 - 2 * x**2 + x
error = 0.005 * scipy.rand(len(x))

# генерим "выпавшие" точки с большими погрешностями
for i in range(30, 60):
    y[i] = y[i] - 2* abs(scipy.rand(1))
    error[i] = 1 * abs(scipy.rand(1))

# находим коэффициенты полинома без учёта погрешностей
fit_coeffs = scipy.polyfit(x, y, 3)
# структура коэффициентов: (x**(max), x**(max-1), x**(max-2), ...)

def fit_func(x, fit_coeffs):
    """
    фит без учёта весов

    """
    fit_func = 0
    # максимальная степень полинома.
    # В данном случае 3 = 4 - 1
    max_degree = len(fit_coeffs)-1

    for i in range(max_degree):
        fit_func += fit_coeffs[i]* x**(max_degree-i)
    return fit_func


# коэффициенты полинома с учетом весов 1/error
# i-й коэффициент соответствует i-й степени полинома
# здесь полная степень полинома = 3
weight_coeffs = polyfitw(x, y, 1/error, 3)

def weight_func(x, weight_coeffs):
    """
    фит с учётом весов

    """
    weight_func = 0

    # обратите внимание на другую структуру вектора weight_coeffs
    # (x**0, x**1, x**2, x**3)
    for i in range(len(weight_coeffs)):
        weight_func += weight_coeffs[i] * x**i
    return weight_func


pl.errorbar(x, y, error, fmt= 'bo')
# график уточнения без "весов"
pl.plot(x, fit_func(x, fit_coeffs), 'r-', \
        linewidth=2, label='no weights')
# график уточнения с "весами"
pl.plot(x, weight_func(x, weight_coeffs), 'g-', \
        linewidth=2, label='with weights')

pl.title('Comparison between two fitting methods')
pl.legend()
pl.show()

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

В результате получим следующую картину:

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

Копируете статью - поставьте ссылку на оригинал!