Gnuplot. МНК

Одно из замечательных свойств программы Gnuplot, это возможность использовать нелинейный метод наименьших квадратов (МНК) для анализа экспериментальных зависимостей. Чтобы ею воспользоваться, нужно указать файл с данными, вид искомой зависимости и ввести команду fit.

Данные в файле могут быть представлены в виде "x y" или "x y s", где x - независимая переменная, y - значение в точке x, s - стандартное отклонение, используемое программой для расчёта "веса" в виде 1/(s**2). Для функции двух переменных воспользуйтесь форматом "x y z" или "x y z s" соответственно.

Искомая зависимость может быть выражена через любые известные Gnuplot функции. В том числе, возможно использование выражений, содержащих условие перехода if.

Поиск коэффициентов запускается выражением

fit f(x) 'file.dat' via 'coeff.par' или 

fit f(x) 'file.dat' via a, b, c ...

где f(x) - заданная функция, file.dat - файл с данными, coeff.par - файл с начальными значениями параметров функции f(x), а a, b, c - искомые параметры. Если значения параметров не заданы, они принимаются равными 1. Однако, часто от близости начальных значений к реальным параметрам зависит адекватность решения, поэтому к их выбору стоит отнестись серьёзно.

Если вы разбираетесь в алгоритме Маркарда-Левенберга, Gnuplot готов предоставить вам возможность настроить параметры расчёта. Однако, большинству пользователей будет достаточно двух: FIT_LIMIT - требуемая точность, FIT_MAXITER - максимальное число итераций.

Рассмотрим пример. Попробуем найти коэффициенты квадратичной зависимости d1 и d2 для массива данных points.dat.

Зададим файл начальных значений коэффициентов d.par.

По каким-то причинам программа отказалась считывать запись "d2 = -0.5", поэтому пришлось ввести её в "урезанном" виде. Закомментированная первая строка отображает искомую зависимость. Теперь нужно запустить Gnuplot, записать функцию и воспользоваться командой fit.

    g(x) = d1*x + d2*x**2

    fit g(x) 'points.dat' via 'd.par'

Можно было обойтись и без файла d.par, просто задав

    fit g(x) 'points.dat' via d1, d2

Каждая итерация поиска решения будет отображена программой. В конце концов, мы получим сообщение

Таким образом, d1 = 1.22, d2 = -0.39, сумма квадратов отклонений теоретической зависимости от экспериментальной 0.0025. Построим график:

    set xrange [0:3]

    plot 'points.dat', g(x)