In [40]:
import math
import random
import numpy as np
from matplotlib import mlab
from matplotlib import pylab as plt
%matplotlib inline
# Построение графика композиции гауссианов с заданными параметрами и определение функции вычисления её значения:
def gaussian(omg,sgm,omg0):
return (1/(sgm*math.sqrt(2*math.pi)))*math.exp(-(omg-omg0)**2/(2*sgm**2))
def gaussComp(omg, sgm, omg0):
r = 0
for i in range(len(sgm)):
r = r + gaussian(omg, sgm[i], omg0[i])
return r
sigma = [1, 1.5]
omega_0 = [0, 10]
dOmega = 0.01
omegaMin = -15
omegaMax = 25
omegaList = mlab.frange (omegaMin, omegaMax, dOmega)
gaussList = [gaussComp(omega, sigma, omega_0) for omega in omegaList]
plt.plot (omegaList, gaussList,'green')
Out[40]:
In [41]:
# выборка из "непрерывной" композиции гауссианов с заданными параметрами тренировочного набора данных:
dOmegaTrain = 0.2
omegaTrain = mlab.frange(omegaMin, omegaMax, dOmegaTrain)
gaussTrain = [gaussComp(omega, sigma, omega_0) for omega in omegaTrain]
plt.plot (omegaTrain, gaussTrain,'x')
#Учебный набор в виде списка списков:
trainData = [[omega] for omega in omegaTrain]
m = len(omegaTrain)
for i in range(m):
trainData[i].append(gaussTrain[i])
#Пусть так же известно и колическо гауссианов в комбинации, так уже легко можно будет построить гипотезу:
gaussN = len(sigma)
#Запишем в отдельную переменную длину учебного набора:
m = len(omegaTrain)
#Параметр определяющий количество последовательных узлов для формирования mini-batch набора
L = 7
In [42]:
def J(sgm, omg0): #Весовая функция, в виде суммы квадратов разности значений гипотезы и заданной функции в узлах. Пакетная.
cost = 0
for i in range(m):
omg = trainData[i][0]
y = trainData[i][1]
cost = cost + (1/(2*m))*(gaussComp(omg, sgm, omg0)-y)**2
return cost
def JSmg_Der(k, j, sgm, omg0): # Частная производная функции J по sigma:
der = 0
for i in range(L):
omg = trainData[k][0]
y = trainData[k][1]
der = der + (gaussComp(omg, sgm, omg0)-y)*((-1/(math.sqrt(2*math.pi)*sgm[j]**2))*math.exp(-(omg-omg0[j])**2/(2*sgm[j]**2)) + gaussian(omg, sgm[j], omg0[j])*((omg-omg0[j])**2/sgm[j]**3))
return der
def Jomg0_Der(k, j, sgm, omg0): # Частная производная функции J по omega0:
der = 0
for i in range(L):
omg = trainData[k][0]
y = trainData[k][1]
der = (gaussComp(omg, sgm, omg0)-y)*gaussian(omg, sgm[j], omg0[j])*(omg-omg0[j])*(1/(sgm[j]**2))
return der
In [43]:
alpha = 0.5 #скорость обучения
# Стартовый набор коэффицентов проинициализируем единицами(нулями нельзя, т.к. сигма встречается в знаменателе, и получим ошибку)
# А так же костыль для временного набора параметров, нулями
#prm = []
#prmTemp0 = []
#for i in range(gaussN):
# prm.append([1, 1])
# prmTemp0.append([0, 0])
prm = [[0.5, -5], [2, 15]]
prmTemp0 = [[0,0],[0,0]]
prmTemp = prmTemp0
for i in range(500000): # основной цикл реализации градиентного спуска
k = random.randint(0,m-1)
tSigma = []
tOmega_0 = []
for v in range(gaussN):
tSigma.append(prm[v][0])
tOmega_0.append(prm[v][1])
for j in range(gaussN):
prmTemp[j][0] = prm[j][0] - alpha * JSmg_Der(k, j, tSigma, tOmega_0)
prmTemp[j][1] = prm[j][1] - alpha * Jomg0_Der(k, j, tSigma, tOmega_0)
prm = prmTemp
prmTemp = prmTemp0
#Разделим лист параметров на отдельные листы для сигма и омега_0.
apprSigma = []
apprOmega_0 = []
for i in range(gaussN):
apprSigma.append(prm[i][0])
apprOmega_0.append(prm[i][1])
In [44]:
apprGauss = [gaussComp(omega, apprSigma, apprOmega_0) for omega in omegaList]
plt.plot(omegaList, apprGauss)
Out[44]:
In [45]:
J(apprSigma, apprOmega_0)
prm
Out[45]:
In [46]:
error = []
for i in range(len(omegaList)):
error.append(gaussList[i] - apprGauss[i])
plt.plot(omegaList, error)
Out[46]: