In [121]:

    

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.5, 1.5]
omega_0 = [0, 35]

dOmega = 0.01
omegaMin = -10
omegaMax = 45

omegaList = mlab.frange (omegaMin, omegaMax, dOmega)
gaussList = [gaussComp(omega, sigma, omega_0) for omega in omegaList]

plt.plot (omegaList, gaussList,'green')


    

    
Out[121]:





[<matplotlib.lines.Line2D at 0x88df710>]






    

In [122]:

    

# выборка из "непрерывной" композиции гауссианов с заданными параметрами тренировочного набора данных:
dOmegaTrain = 0.4

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 [123]:

    

#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, sgm, omg0): # Частная производная функции J по sigma в точке:
    omg = trainData[k][0]
    y = trainData[k][1]
    der = (gaussian(omg, sgm, omg0)-y)*((-1/(math.sqrt(2*math.pi)*sgm**2))*math.exp(-(omg-omg0)**2/(2*sgm**2)) + gaussian(omg, sgm, omg0)*((omg-omg0)**2/sgm**3))
    return der

def Jomg0_Der(k, sgm, omg0): # Частная производная функции J по omega0 в точке:
    omg = trainData[k][0]
    y = trainData[k][1]
    der = (gaussian(omg, sgm, omg0)-y)*gaussian(omg, sgm, omg0)*(omg-omg0)*(1/(sgm**2))
    return der


    

In [124]:

    

len(trainData)


    

    
Out[124]:





139

In [125]:

    

alpha = 0.5 #скорость обучения

# Стартовый набор коэффицентов проинициализируем единицами(нулями нельзя, т.к. сигма встречается в знаменателе, и получим ошибку)
# А так же костыль для временного набора параметров, нулями
#prm = []
#prmTemp0 = []
#for i in range(gaussN):
#    prm.append([1, 1])
#    prmTemp0.append([0, 0])

prm = [[1, 2], [3, 35]]

prmTemp0 = [[0,0],[0,0]]    
    


prmTemp = prmTemp0
for i in range(100000): # основной цикл реализации градиентного спуска
    k = random.randint(0,m-1)
    for j in range(gaussN):  
        prmTemp[j][0] = prm[j][0] - alpha * JSmg_Der(k, prm[j][0], prm[j][1])
        prmTemp[j][1] = prm[j][1] - alpha * Jomg0_Der(k, prm[j][0], prm[j][1])
    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 [126]:

    

apprGauss = [gaussComp(omega, apprSigma, apprOmega_0) for omega in omegaList]
plt.plot(omegaList, apprGauss)


    

    
Out[126]:





[<matplotlib.lines.Line2D at 0x890c970>]






    

In [127]:

    

J(apprSigma, apprOmega_0)
prm


    

    
Out[127]:





[[1.4999999999999996, 1.3878995883212396e-16],
 [1.5000000000000528, 35.000000000000611]]

In [128]:

    

error = []
for i in range(len(omegaList)):
    error.append(gaussList[i] - apprGauss[i])
plt.plot(omegaList, error)


    

    
Out[128]:





[<matplotlib.lines.Line2D at 0x87a7910>]






    

In [ ]: