Простая программа на Python для гиперболической аппроксимации статистических данных

Зачем это нужно


Законы Зипфа оописывают закономерности частотного распределения слов в тексте на любом естественном языке[1]. Эти законы кроме лингвистики применяться также в экономике [2]. Для аппроксимации статистических данных для объектов, которые подчиниться Законам Зипфа используется гиперболическая функция вида:

(1)

где: a.b – постоянные коэффициенты: x – статистические данные аргумента функции (в виде списка): y- приближение значений функции к реальным данным полученным методом наименьших квадратов[3].

Обычно для аппроксимации гиперболической функцией методом логарифмирования её приводят к линейной, а затем определяют коэффициенты a,b и делают обратное преобразование [4]. Прямое и обратное преобразование приводит к дополнительной погрешности аппроксимации. Поэтому привожу простую программу на Python, для классической реализации метода наименьших квадратов.

Алгоритм



Исходные данные задаются двумя списками где n ─ количество данных в списках.

Получим функцию для определения коэффициентов

Коэффициенты a,b найдём из следующей системы уравнений:(3)

Решение такой системы не представляет особого труда:

(4),
(5).

Средняя ошибка аппроксимации по формуле: (6)

Код Python


#!/usr/bin/python
# -*- coding: utf-8 -*
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.family'] = 'fantasy'
def mnkGP(x,y): # функция которую можно использзовать в програме
              n=len(x) # количество элементов в списках
              s=sum(y) # сумма значений y
              s1=sum([1/x[i] for i in  range(0,n)]) #  сумма 1/x
              s2=sum([(1/x[i])**2 for i in  range(0,n)]) #  сумма (1/x)**2
              s3=sum([y[i]/x[i]  for i in range(0,n)])  # сумма y/x                   
              a= round((s*s2-s1*s3)/(n*s2-s1**2),3) # коэфициент а с тремя дробными цифрами
              b=round((n*s3-s1*s)/(n*s2-s1**2),3)# коэфициент b с тремя дробными цифрами
              s4=[a+b/x[i] for i in range(0,n)] # список значений гиперболической функции              
              so=round(sum([abs(y[i] -s4[i]) for i in range(0,n)])/(n*sum(y))*100,3)   # средняя ошибка аппроксимации
              plt.title('Аппроксимация гиперболой Y='+str(a)+'+'+str(b)+'/x\n Средняя ошибка--'+str(so)+'%',size=14)
              plt.xlabel('Координата X', size=14)
              plt.ylabel('Координата Y', size=14)
              plt.plot(x, y, color='r', linestyle=' ', marker='o', label='Data(x,y)')
              plt.plot(x, s4, color='g', linewidth=2, label='Data(x,f(x)=a+b/x')
              plt.legend(loc='best')
              plt.grid(True)
              plt.show()
x=[10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86] 
y=[0.1, 0.0714, 0.0556, 0.0455, 0.0385, 0.0333, 0.0294, 0.0263, 0.0238, 0.0217,
   0.02, 0.0185, 0.0172, 0.0161, 0.0152, 0.0143, 0.0135, 0.0128, 0.0122, 0.0116] # данные для проверки по функции y=1/x
mnkGP(x,y)

Результат




Мы брали данные из функции равносторонней гиперболы, поэтому и получили a=0,b=10 и абсолютную погрешность в 0,004%. Значить функция mnkGP(x,y) работает правильно и её можно вставлять в прикладную программу

Аппроксимация для степенных функций


Для этого в Python есть модуль scipy, но он не поддерживает отрицательную степень d полинома. Рассмотрим код реализации аппроксимации данных полиномом.

#!/usr/bin/python
# coding: utf8
import scipy as sp
import matplotlib.pyplot as plt
def mnkGP(x,y):
              d=2 # степень полинома
              fp, residuals, rank, sv, rcond = sp.polyfit(x, y, d, full=True) # Модель
              f = sp.poly1d(fp) # аппроксимирующая функция
              print('Коэффициент -- a %s  '%round(fp[0],4))
              print('Коэффициент-- b %s  '%round(fp[1],4))
              print('Коэффициент -- c %s  '%round(fp[2],4))
              y1=[fp[0]*x[i]**2+fp[1]*x[i]+fp[2] for i in range(0,len(x))] # значения функции a*x**2+b*x+c
              so=round(sum([abs(y[i]-y1[i]) for i in range(0,len(x))])/(len(x)*sum(y))*100,4) # средняя ошибка
              print('Average quadratic deviation '+str(so)) 
              fx = sp.linspace(x[0], x[-1] + 1, len(x)) # можно установить вместо len(x) большее число для интерполяции
              plt.plot(x, y, 'o', label='Original data', markersize=10)
              plt.plot(fx, f(fx), linewidth=2)
              plt.grid(True)
              plt.show()

x=[10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86] 
y=[0.1, 0.0714, 0.0556, 0.0455, 0.0385, 0.0333, 0.0294, 0.0263, 0.0238, 0.0217,
   0.02, 0.0185, 0.0172, 0.0161, 0.0152, 0.0143, 0.0135, 0.0128, 0.0122,
   0.0116] # данные для проверки по функции y=1/x
mnkGP(x,y)

Результат




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

Полученные функции будут применены для анализа Законов Зипфа (Ципфа), но это будет сделано уже в следующей статье.

Ссылки:

1. Законы Зипфа (Ципфа) tpl-it.wikispaces.com/Законы+Зипфа+%28Ципфа%29
2. Закон Ципфа это. dic.academic.ru/dic.nsf/ruwiki/24105
3. Законы Зипфа. wiki.webimho.ru/законы-зипфа
4. Лекция 5. Аппроксимация функций по методу наименьших квадратов. mvm-math.narod.ru/Lec_PM5.pdf
5. Средняя ошибка аппроксимации. math.semestr.ru/corel/zadacha.php

Similar posts

Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 9

    0
    Спасибо, долго искал решение задачи, и получил грамотно составленный ответ на него! Удачи в работе
      0
      Что-то мне подсказывает, что эта задача метода наименьших квадратов решается для обобщенного полинома более-менее легко и непринужденно.
        +4
        В питоне принято вместо
        s1=sum([1/x[i] for i in  range(0,n)])
        

        писать
        s1 = sum(1 / value for value in x)
        

        Но ещё лучше использовать numpy:
        import numpy as np
        x, y = np.array(x), np.array(y)
        s1 = np.sum(1 / x)
        s2 = np.sum(1 / y)
        s3 = np.sum(y / x)
        ...
        

          +3
          Для этого в Python есть модуль scipy, но он не поддерживает отрицательную степень d полинома. Рассмотрим код реализации аппроксимации данных полиномом.

          Не обязательно использовать линейный метод наименьших квадратов. Можно и нелинейный использовать, который, правда, уже так просто не реализуешь. Но функция же есть! :)


          from scipy.optimize import curve_fit
          
          def f(x, a, b):
              return a + b/x
          
          x=[10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86]
          
          y=[0.1, 0.0714, 0.0556, 0.0455, 0.0385, 0.0333, 0.0294, 0.0263, 0.0238, 0.0217, 0.02, 0.0185, 0.0172, 0.0161, 0.0152, 0.0143, 0.0135, 0.0128, 0.0122, 0.0116]
          
          popt, _ = curve_fit(f, x, y)
          
          fit_y = [f(xi, popt[0], popt[1]) for xi in x]
          
          import matplotlib.pyplot as plt
          plt.plot(x, y, 'o', x, fit_y, '-')
          

          Скрытый текст
            0
            А как питон в плане производительности при больших объёмах данных?
              0
              Такой же как и любой интерпретируемый язык — если удается сводить обработку к редким (относительно) вызовам сишных функций, то все хорошо, в противном случае — не очень
                +2
                Вполне себе отлично, даже на вполне себе риалтайме уровня 15 000 000 событий в секунду, доставляемые Kafka в предсказательную модель на временных рядах.
                Но кто-то и Hadoop разворачивает для поиска буковок в 15 Гб данных, где отлично справляется grep.
                  0
                  Ну у меня следующая проблема: порой нужны фиты для набора в 100к точек и более. Ваять код на С было не очень удобно, я решил воспользоваться гнуплотом. Работает отлично, но функционал ограничен. Вот я и думаю: если питон будет это делать так же быстро (порядка нескольких секунд), то быть может стоит попробовать?
                    0

                    numpy и scipy написаны на C и Fortran, всё работает довольно быстро, поищите бенчмарки в сети.
                    Если данные обрабатывать на чистом питоне, можно накрутить JIT в виде numba-декораторов. Смотрите инфу по пакету numba.

              Only users with full accounts can post comments. Log in, please.