In [3]:
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 = [2, 0.5, 2]
omega_0 = [0, 10, 20]
dOmega = 0.01
omegaMin = -15
omegaMax = 30
omegaList = mlab.frange (omegaMin, omegaMax, dOmega)
gaussList = [gaussComp(omega, sigma, omega_0) for omega in omegaList]
plt.plot (omegaList, gaussList,'green')
Out[3]:
In [199]:
# выборка из "непрерывной" композиции гауссианов с заданными параметрами тренировочного набора данных:
dOmegaTrain = 0.5
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)
In [200]:
#def J(k, sgm, omg0): #Весовая функция, в виде суммы квадратов разности значений гипотезы и заданной функции в узлах. Для одного узла:
# omg = trainData[k][0]
# y = trainData[k][1]
# cost = (gaussComp(omg, sgm, omg0)-y)**2
# return cost
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:
omg = trainData[k][0]
y = trainData[k][1]
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:
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 [ ]:
In [201]:
alpha = 2 #скорость обучения
# Стартовый набор коэффицентов проинициализируем единицами(нулями нельзя, т.к. сигма встречается в знаменателе, и получим ошибку)
# А так же костыль для временного набора параметров, нулями
#prm = []
#prmTemp0 = []
#for i in range(gaussN):
# prm.append([1, 1])
# prmTemp0.append([0, 0])
prm = [[3, -5], [1, 15], [10, 24]]
prmTemp0 = [[0,0],[0,0],[0,0]]
prmTemp = prmTemp0
for i in range(150000): # основной цикл реализации градиентного спуска
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 [202]:
apprGauss = [gaussComp(omega, apprSigma, apprOmega_0) for omega in omegaList]
plt.plot(omegaList, apprGauss)
Out[202]:
In [203]:
J(apprSigma, apprOmega_0)
prm
Out[203]:
In [204]:
error = []
for i in range(len(omegaList)):
error.append(gaussList[i] - apprGauss[i])
plt.plot(omegaList, error)
Out[204]:
In [7]:
A = mlab.frange (-5, 15, dOmega)
B = [gaussComp(omega, [1, 1.5], [1, 5]) for omega in A]
plt.plot(A,B)
Out[7]: