Tutorial of linear regresion in www.bigneuralvision.com

tutorial_linear_regresion.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D

# numero iteraciones
ITERS = 10000

# learning rate
LR = 0.0000000001

# Create GIF
GIF = True

# Generamos datos sobre el salario y alquiler.
SAMPLES = 1000
FACTOR = 0.6 + np.random.random_sample(size=[SAMPLES])*0.2
BIAS = np.random.random_sample(size=[SAMPLES])*100

# Espacio de soluciones
SOL_POINTS_SHOW_A = 50
SOL_POINTS_SHOW_B = 100
SOL_RANGE_A = [-6, 6]
SOL_RANGE_B = [-1000, 1000]

salario = np.random.randint(500, 3000, size=[SAMPLES])
alquiler = np.round(FACTOR*salario + BIAS).astype('float64')

# Modelo
x = salario.astype('float64')
y = alquiler.astype('float64')
AB = np.random.random_sample(size=[2]).astype('float64')

# Gradient descent (MUY LENTO)
"""for iter in range(ITERS):
	dAB = np.zeros(shape=AB.shape[0], dtype=AB.dtype)
	for i in range(x.shape[0]):
		x_ampl = np.append(x[i], 1.0)
		pred = np.sum(AB * x_ampl)
		dAB += (pred - y[i]) * x_ampl
	dAB /= x.shape[0]
	AB -= LR*dAB
print(AB)"""

# Gradient descent (MAS RAPIDO)
if GIF:
	AB_HIST = np.zeros(shape=(ITERS, AB.shape[0]), dtype=AB.dtype)
for iter in range(ITERS):
	dAB = np.zeros(shape=AB.shape[0], dtype=AB.dtype)
	pred = (AB[0] * x) + AB[1]
	diff = (pred - y)
	dAB = np.array([np.sum(diff*x), np.sum(diff*1)])
	dAB /= x.shape[0]
	AB -= LR*dAB

	if GIF:
		AB_HIST[iter,:] = AB
print(AB)

# PLOT ESPACIO DE SOLUCIONES (recortado)
fig = plt.figure()
def showSolutionSpace():
	# Generamos los rangos de A y B.
	A = np.linspace(SOL_RANGE_A[0], SOL_RANGE_A[1], num=SOL_POINTS_SHOW_A)
	B = np.linspace(SOL_RANGE_B[0], SOL_RANGE_B[1], num=SOL_POINTS_SHOW_B)
	A, B = np.meshgrid(A, B, sparse=False)
	# Cambiamos el tamaño para que se puedan aplicar la formula de manera matricial
	Aaux = np.tile(np.reshape(A, [SOL_POINTS_SHOW_A, SOL_POINTS_SHOW_B, 1]), [1, 1, SAMPLES])
	Baux = np.tile(np.reshape(B, [SOL_POINTS_SHOW_A, SOL_POINTS_SHOW_B, 1]), [1, 1, SAMPLES])
	X = np.tile(np.reshape(x, [1, 1, SAMPLES]), [SOL_POINTS_SHOW_A, SOL_POINTS_SHOW_B, 1])
	Y = np.tile(np.reshape(y, [1, 1, SAMPLES]), [SOL_POINTS_SHOW_A, SOL_POINTS_SHOW_B, 1])
	pred = Aaux*X + Baux
	# Se calcula el error
	error = np.sum((pred - Y)**2, axis=-1)
	# Se muestra el espacio de soluciones
	ax = Axes3D(fig)
	ax.scatter(A, B, error)
	ax.set_xlabel("A")
	ax.set_ylabel("B")
	ax.set_zlabel("Error")
showSolutionSpace()
plt.show()

# PLOT RESULTADO FINAL
def plotFrame(AB):
	plt.clf()
	plt.plot(x, y, 'bo')
	rango = np.arange(np.min(x), np.max(x))
	plt.plot(rango, AB[0]*rango + AB[1], 'r')
fig = plt.figure()
plotFrame(AB)
plt.show()


# GIF PROCESO
if GIF:
	fig = plt.figure()
	def gifFrame(i):
	    label = 'salario\niteracion {0}'.format(i)
	    plotFrame(AB_HIST[i])
	    plt.ylabel("alquiler")
	    plt.xlabel(label)
	
	anim = FuncAnimation(fig, gifFrame, frames=np.arange(0, ITERS, ITERS//10), interval=200)
	anim.save('regresion_lineal.gif', dpi=80, writer='imagemagick')