Inferencia Bayesiana: MLE y MAP

Cuando tratamos con modelos probabilisticos indirectamente siempre hablamos de MAP y MLE. ¿Pero que son exactamente?
Para ello debemos empezar por el Teorema de Bayes:
$\underbrace{P(\sigma | X)}_{\text{posterior}} = \frac{\overbrace{P(X | \sigma)}^{\text{likelihood } L(\sigma | X)} \cdot \overbrace{P(\sigma)}^{\text{prior}}}{\underbrace{P(X)}_{marginal}}$

La formula nos permite «girar» la probabilidad condicionada. Esto es muy interesante cuando planteamos esta formula en términos de entrada de datos $X$ y parametros $\sigma$ . $P(X | \sigma)$ nos esta diciendo: Que probabilidad tiene esa entrada $X$ , teniendo en cuenta que los parametros son $\sigma$ y $P(\sigma | X)$ : Cuanto de probables son esos parametros $\sigma$ si la entrada es $X$ .
Inicialmente cuando trabajamos con modelos probabilisticos no queremos exactamente esta probabilidad sino que ese será nuestro modelo final. Sino en que el teorema de bayes es la vía que nos permite estimar los parametros $\sigma$ .
Aquí es donde los conceptos MAP y MLE entran en juego, ambos son estimadores y tratan de resolver el siguiente problema:
$\sigma = {arg max}_{\sigma} P(X | \sigma)$ ¿Pero que les diferencia uno del otro? Simple, las restricciones que se aplican.
$P(\sigma | X) \propto P(X | \sigma) \cdot P(\sigma)$ . Podemos decir que $P(\sigma | X)$ es proporcional a $P(X | \sigma) \cdot P(\sigma)$ , porque el denominador será constante si X no varia. Así que podemos reescribir el problema en:
$\sigma = {arg max}_\sigma P(X | \sigma) \cdot P(\sigma)$

MLE considera que $P(\sigma)$ es equiprobable, es decir tiene una distribución uniforme. Se utiliza cuando no sabemos nada sobre $P(\sigma)$ .
$\sigma_{MLE} = {arg max}_\sigma P(X | \sigma)$ $\forall \sigma : P(\sigma) = \frac{1}{|\sigma|}$
MAP considera que $P(\sigma)$ es conocida y por lo tanto podemos aportar más información al sistema y tener un modelo mas preciso.
$\sigma_{MAP} = {arg max}_\sigma P(X | \sigma) \cdot P(\sigma)$

Ambas formulas són extremadamente potente y es que se puede usar tanto para estimar $\sigma$ categóricos como numéricas. Normalmente por Internet las encontráis como MAP vs MLE y eso lleva a confusión ya que no son algoritmos distintos. Si conocéis $P(\sigma)$ siempre sera mejor MAP.
No acostumbra a ser factible usar MAP cuando lo que estamos buscando es generar un modelo, no podemos conocer $P(\sigma)$ .
Vamos a plantear dos ejemplos:

Realizar un naive bayes usando MAP.
Realizar una regresión logistica usando MLE.

Realizar un naive bayes usando MAP:

Vamos a realizar un algoritmo basando en MAP y Naive Bayes que permitirá clasificar si un email es spam o no. Este ejemplo es ampliamente conocido en el sector. Es extremadamente simple y es bastante bueno.
Vamos a crear un modelo probabilistico basado en palabras, las palabras tendrán asociadas tu probabilidad de ser o no spam. Por ejemplo:
$P(\sigma=spam | X_1=te, X_2=vendo, X_3=viagra) =P(X_1=te, X_2=vendo, X_3=viagra | \sigma=spam)\cdot P(\sigma=spam) = 1$
$P(\sigma = spam | X_1=hola, X_2=cariño, X_3=te, X_4=has, X_5=acordado, X_6=de, X_7=coger, X_8=la, X_9=viagra) = \ P(X_1=hola, X_2=cariño, X_3=te, X_4=has, X_5=acordado, X_6=de, X_7=coger, X_8=la, X_9=viagra | \sigma = spam) \cdot P(\sigma=spam)$
El problema que pretenderíamos resolver seria:
$P(\sigma |X_1=palabra_1, X_2=palabra_2, ..., X_n=palabra_n) \propto \ P(X_1=palabra_1, X_2=palabra_2, ..., X_n=palabra_n | \sigma) \cdot P(\sigma)$
Donde $\sigma$ puede ser spam o no-spam.
Hasta aquí ya tenemos un problema gordo… ¡Que narices es $<span style="font-size: 10pt;">P(X_1=palabra_1, X_2=palabra_2, ..., X_n=palabra_n | \sigma)$ !
Nos dice que probabilidad tiene una frase de ser spam o no. ¿Como se materializa?
Como sabréis si tenéis algunos conocimientos de probabilidad, no es tan fácil sacar la probabilidad condicionada respecto $\sigma$ si las palabras son dependientes entre sí.
Aquí es donde entra en juego Naive Bayes, que su traducción seria algo similar a el «Bayes ingenuo». Esté va a asumir que todas las palabras entre ellas son independientes. Esto es una asunción muy fuerte porque sabemos claramente que en la realidad no es así y que la preposición «de» si va delante de «viagra» probablemente la condicione a ser spam.
El problema se debe a que calcular todas las relaciones hacen que el modelo sea infinitamente complejo y si hacemos esa suposición lo hacemos totalmente viable. Pero, ¿Considerando esa suposición podemos encontrar un modelo suficientemente bueno? Cuando trabajamos en problemas de machine learning buscamos un compromiso entre esfuerzo, rendimiento y complejidad. Entonces podemos decir que:
$P(X_1=te, X_2=vendo, X_3=viagra | \sigma=spam) = \\ P(\sigma=spam)P(X_1=te)P(X_2=vendo)P(X_3=viagra)$
$P(X_1=hola, X_2=cariño, X_3=te, \ X_4=has, X_5=acordado, X_6=de, X_7=coger, X_8=la, X_9=viagra |\sigma = spam) = \\ P(\sigma=spam)P(X_1=hola)P(X_2=cariño)P(X_3=te)P(X_4=has)P(X_5=acordado)P(X_6=de)P(X_7=coger)P(X_8=viagra)$
Así que la formula del Naive Bayes queda de la siguiente forma:
$\sigma = argmax_{\sigma \in \text{( spam, no-spam)}} P(\sigma) \prod_i^N P(X_i=palabra_i | \sigma)$
Al ser las palabras independientes entre sí, podemos estimar su distribución obteniendo la probabilidad frecuencial de una base de datos en la que tengamos emails y si el email es spam o no.
$P(palabra_i|\sigma=spam) = \frac{P(palabra_i, \sigma=spam)}{P(\sigma=spam)}= \frac{\text{numero palabra_i en emails spam}}{\text{numero palabras totales en spam}}$
Esta técnica se emplea mucho en machine learning y se conoce como Bag of Words. Al leer todos los emails construimos dos vocabularios (spam y no-spam) con su numero de veces asociado. Con esto conseguimos las frecuencias absolutas de cada palabra para luego convertirlas en probabilidades.

¡Al código!

import pandas as pd
import numpy as np
from math import log
import string
from sklearn.metrics import confusion_matrix, accuracy_score
# dataset:
# https://www.kaggle.com/uciml/sms-spam-collection-dataset
dataset = pd.read_csv('spam.csv')[['v1', 'v2']]
# Distribuciones de los messages
# Probabilidad de ser spam o no
# Calculamos cuantos mensajes son spam y cuantos no
P ={}
P['spam'] = float(len(dataset['v1'] == 'spam')) / len(dataset)
P['ham'] = float(len(dataset['v1'] == 'ham')) / len(dataset)
# Convierte la palabra a ASCII
# La pone en minusculas
# Elimina signos de puntuacion
# Preprocesados mas potentes => Tokenizers
def preprocessWord(word):
    word = unicode(word.lower(), errors='ignore')
    word = str(word.decode('ascii', 'ignore'))
    word = word.translate(None, string.punctuation)
    return word
# Distribuciones de palabras
# Frecuencia de ser una palabra X teniendo en cuenta que es SPAM / HAM
# En la siguiente funcion convertiremos esta frecuencia absoluta en una probabilidad
wordCount = {'ham': {},'spam': {}}
totalWords = {'ham': 0, 'spam': 0}
for messageBox in dataset.iterrows():
    spamOrHam = messageBox[1]['v1']
    message = messageBox[1]['v2']
    for word in message.split(" "):
        word = preprocessWord(word)
        if word not in wordCount[spamOrHam]:
            wordCount[spamOrHam][word] = 1
        else:
            wordCount[spamOrHam][word] += 1
        totalWords[spamOrHam] += 1
numWords = sum(totalWords.values())
# Probabilidad Condicionadas de las palabras
# P(palabra | spam)
# Se añade el concepto de laplace smooth para evitar que palabras desconocidas devuelvan 0. log(0) = big fail
# https://en.wikipedia.org/wiki/Additive_smoothing
L = 0.00001 # laplace smooth
def PConditional(word, class_w):
    return float(wordCount[class_w].get(word, 0) + L)/(totalWords[class_w] + L*numWords)
# Model
# Se usa un modelo MAP, pero se maximiza usando el log-likelihood para evitar problemas numericos
# Eso se debe a que el producto continuado de probabilidades muy pequeñas acaba siendo un valor inifinitamente pequeño
# y en terminos computacionales somos incapaces de representarlo por lo que acaba valiendo 1.
# Para evitar esto se introduce el logaritmo y se convierte el productorio en un sumatorio.
# Using log-maximize, avoid low values product
def is_spam(sentence):
    acc = {'spam': log(P['spam']), 'ham': log(P['ham'])}
    for word in sentence.split(" "):
        word = preprocessWord(word)
        acc['spam'] += log(PConditional(word, 'spam'))
        acc['ham'] += log(PConditional(word, 'ham'))
    return acc['spam'] > acc['ham']
# Obtenemos los resultados
# Accuracy of training set
true = []
pred = []
for mb in dataset.iterrows():
    spamOrHam = mb[1]['v1']
    message = mb[1]['v2']
    prediction_is_spam = is_spam(message)
    true.append(spamOrHam == 'spam')
    pred.append(prediction_is_spam)
print('Accuracy:', accuracy_score(true, pred))
print('Confusion Matrix:\n')
print(confusion_matrix(true, pred))

import pandas as pd

import numpy as np

from math import log

import string

from sklearn.metrics import confusion_matrix, accuracy_score

# dataset:

# https://www.kaggle.com/uciml/sms-spam-collection-dataset

dataset = pd.read_csv('spam.csv')[['v1', 'v2']]

# Distribuciones de los messages

# Probabilidad de ser spam o no

# Calculamos cuantos mensajes son spam y cuantos no

P ={}

P['spam'] = float(len(dataset['v1'] == 'spam')) / len(dataset)

P['ham'] = float(len(dataset['v1'] == 'ham')) / len(dataset)

# Convierte la palabra a ASCII

# La pone en minusculas

# Elimina signos de puntuacion

# Preprocesados mas potentes => Tokenizers

def preprocessWord(word):

word = unicode(word.lower(), errors='ignore')

word = str(word.decode('ascii', 'ignore'))

word = word.translate(None, string.punctuation)

return word

# Distribuciones de palabras

# Frecuencia de ser una palabra X teniendo en cuenta que es SPAM / HAM

# En la siguiente funcion convertiremos esta frecuencia absoluta en una probabilidad

wordCount = {'ham': {},'spam': {}}

totalWords = {'ham': 0, 'spam': 0}

for messageBox in dataset.iterrows():

spamOrHam = messageBox[1]['v1']

message = messageBox[1]['v2']

for word in message.split(" "):

word = preprocessWord(word)

if word not in wordCount[spamOrHam]:

wordCount[spamOrHam][word] = 1

else:

wordCount[spamOrHam][word] += 1

totalWords[spamOrHam] += 1

numWords = sum(totalWords.values())

# Probabilidad Condicionadas de las palabras

# P(palabra | spam)

# Se añade el concepto de laplace smooth para evitar que palabras desconocidas devuelvan 0. log(0) = big fail

# https://en.wikipedia.org/wiki/Additive_smoothing

L = 0.00001 # laplace smooth

def PConditional(word, class_w):

return float(wordCount[class_w].get(word, 0) + L)/(totalWords[class_w] + L*numWords)

# Model

# Se usa un modelo MAP, pero se maximiza usando el log-likelihood para evitar problemas numericos

# Eso se debe a que el producto continuado de probabilidades muy pequeñas acaba siendo un valor inifinitamente pequeño

# y en terminos computacionales somos incapaces de representarlo por lo que acaba valiendo 1.

# Para evitar esto se introduce el logaritmo y se convierte el productorio en un sumatorio.

# Using log-maximize, avoid low values product

def is_spam(sentence):

acc = {'spam': log(P['spam']), 'ham': log(P['ham'])}

for word in sentence.split(" "):

word = preprocessWord(word)

acc['spam'] += log(PConditional(word, 'spam'))

acc['ham'] += log(PConditional(word, 'ham'))

return acc['spam'] > acc['ham']

# Obtenemos los resultados

# Accuracy of training set

true = []

pred = []

for mb in dataset.iterrows():

spamOrHam = mb[1]['v1']

message = mb[1]['v2']

prediction_is_spam = is_spam(message)

true.append(spamOrHam == 'spam')

pred.append(prediction_is_spam)

print('Accuracy:', accuracy_score(true, pred))

print('Confusion Matrix:\n')

print(confusion_matrix(true, pred))

Como resultado obtenemos para todo el dataset:

Verdadero \ Predicción	Ham	Spam
Ham	4809	16
Spam	2	745

Con un accuracy de 99.6%. ¿Nada mal para lo simple que es, no?

Realizar una regresión logística usando MLE:

En este caso usaremos MLE para crear un regresor logístico. En este caso lo usaremos para estimar si un paciente tiene cáncer o no a partir de una base de datos: https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Coimbra
Utilizar el MLE para resolver este problema es similar a cuando usamos una función de coste o perdida, pero en vez de buscar el minimo buscaremos el maximo.
Un regresor logistico es una estimador que se usa para estimar la probabilidad de una variable categorica (pertenecer o no, ser o no). Son regresores muy utiles para datos que siguen una distribución de Bernoulli ( $X \sim B(n, p)$ ).

¿Porque se llama regresor logistico?

Se estima la función logit mediante una función lineal:
$\underbrace{ln \left ( \underbrace{\frac{\overbrace{p}^{probabilidad}}{1-p}} _ {odd} \right )} _ {logit} = W \cdot X + B$ Un logit no es más que el logaritmo de un odd (no importa la base mientras sea superior 1, aunque normalmente la base es $e$ ). Un odd es la proporción entre el exito y fracaso de un evento binario.
En este tipo de modelo tenemos conceptos interesantes:

Hemos relacionado la probabilidad de un evento binaria con una función facil de computar.
El odd convierte la probabilidad en un valor de 0 a infinito.
El logaritmo convierte el odd en un valor de entre -inf a inf. Un rango que podemos aproximar facilmente mediante una función lineal.Podemos aislar $p$ y podremos estimar directamente la probabilidad.
$P(x | w, b) = \frac{1}{1 + e^{- W \cdot X + B}}$ ¡Si of fijais es la función sigmoide!¿Cual es la función que vamos a maximizar?

Vamos a definir la función de likelihood que mediante MLE seremos capaces de encontrar los mejores valores de $W$ y $B$ .
$\mathcal{L}(W,B | X_i) = \left [ sigmoid(X_i \cdot X + B) \right ]^{ y_i }\left [ 1 - sigmoid(X_i \cdot X + B) \right ]^{\left ( 1 - y_i \right ) }$ $W_{MLE} = argmax_{W} \prod_i \underbrace{P(X_i | W, B)} _ {\mathcal{L}(W,B | X_i)}$ Podemos aplicar el log-likelihood que vimos en el anterior ejemplo para simplificar el modelo.
$\mathcal{L}(W,B | X_i) = \left [ sigmoid(X_i \cdot X + B) \right ]^{ y_i }\left [ 1 - sigmoid(X_i \cdot X + B) \right ]^{\left ( 1 - y_i \right ) }$ $W_{MLE} = argmax_{W} ln \left ( \sum_i \mathcal{L}(W,B | X_i) \right )$ $log\mathcal{L}(W,B | X_i) = ln \left ( \left [ sigmoid(X_i \cdot X + B) \right ]^{ y_i }\left [ 1 - sigmoid(X_i \cdot X + B) \right ]^{\left ( 1 - y_i \right ) } \right )$ $W_{MLE} = argmax_{W} \left ( \sum_i ln\mathcal{L}(W,B | X_i) \right )$ Ahora nos centramos en $log\mathcal{L}(W,B|X_i)$ (1 muestra). Para trabajar mejor diremos que $sigmoid = \frac{1}{1+e^{f}}$ y la variable $f = -WX+B$ .
Ahora separaremos log-likelihood en pequeños terminos :
$= ln \left ( \frac{1}{1 + e^{f}} \right ) \cdot y + \underbrace{ln \left (1- \frac{1}{1 + e^{f}} \right ) \cdot (1 - y)} _ {\underbrace{ln \left (\frac{\cancel{1}+e^f \cancel{- 1}}{1 + e^{f}} \right ) \cdot (1 - y)} _ {ln \left (\frac{e^f}{1 + e^{f}} \right ) - ln \left (\frac{e^f}{1 + e^{f}} \right ) \cdot y}}$ $=ln \left ( \frac{1}{1 + e^{f}} \right ) \cdot y + ln \left (\frac{e^f}{1 + e^{f}} \right ) - ln \left (\frac{e^f}{1 + e^{f}} \right ) \cdot y$ Agrupamos por $y$ :
$=\left [ln \left ( \frac{1}{1 + e^{f}} \right ) - ln \left (\frac{e^f}{1 + e^{f}} \right ) \right ]\cdot y + ln \left (\frac{e^f}{1 + e^{f}} \right )$ Simplificamos el termino izquierdo:
$=\left [ln \left ( \frac{1}{\cancel{1 + e^{f}}} \frac{\cancel{1 + e^f}}{e^{f}} \right ) \right ]\cdot y + ln \left (\frac{e^f}{1 + e^{f}} \right )$ $=ln \left ( \frac{1}{e^f} \right ) \cdot y + ln \left (\frac{e^f}{1 + e^{f}} \right )$ $=- f\cdot y + ln \left (\frac{e^f}{1 + e^{f}} \right )$ Simplificamos el termino derecho:
$=- f\cdot y + ln \left (\frac{e^f}{1 + e^{f}} \frac{e^{-f}}{e^{-f}} \right )$ $=- f\cdot y + ln \left (\frac{1}{1 + e^{-f}} \right )$ $=- f\cdot y - ln \left (1 + e^{-f} \right )$ Substituimos $f$ por $-WX+B$ :
$=(WX + B)\cdot y - ln \left ( 1 + e^{WX+B} \right )$ Nuestra función a MLE queda como:
$W_{MLE} = argmax_{W} \sum_i\left ((WX_i + B)\cdot y_i - ln \left (1 + e^{WX_i+B} \right ) \right )$ Esta será la función que debemos maximizar, para no tener que derivar dos veces. Haremos un truco ya usado en estos articulos, expandiremos el vector $X_i$ con un 1 y así $W$ incluirá el bias ( $B$ ).
$W_{MLE} = argmax_{W} \sum_i\left ((W\begin{bmatrix} X_i\\ 1 \end{bmatrix})\cdot y_i + ln \left (1 + e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}} \right ) \right )$ Derivemos con respecto $W$ :
$\nabla_W = \sum_i \begin{bmatrix} X_i\\ 1 \end{bmatrix} y_i - \frac{1}{1 + e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}} \left (1 + e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}} \right ) \begin{bmatrix} X_i\\ 1 \end{bmatrix}$ Agrupamos con respecto $y_i$ :
$\nabla_W = \sum_i \begin{bmatrix} X_i\\ 1 \end{bmatrix} \left [ y_i - \frac{1}{1 + e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}}\left (1 + e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}} \right ) \right ]$ Simplificamos:
$\nabla_W = \sum_i \begin{bmatrix} X_i\\ 1 \end{bmatrix} \left [ y_i - \frac{\cancel{e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}}}{1 + e^{W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}} \frac{\cancel{e^{-W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}}}{e^{-W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}} \right ]$ $\nabla_W = \sum_i \begin{bmatrix} X_i\\ 1 \end{bmatrix} \left [ y_i -\underbrace{\frac{1}{1 + e^{-W\begin{bmatrix} X_i\\ 1 \end{bmatrix}}}}_{sigmoid \left ( W\begin{bmatrix} X_i\\ 1 \end{bmatrix} \right )} \right ]$ $\nabla_W =\sum_i \begin{bmatrix} X_i\\ 1 \end{bmatrix} \left [ y_i - sigmoid \left ( W\begin{bmatrix} X_i\\ 1 \end{bmatrix} \right ) \right ]$ Como sucede habitualmente en problemas de ML, no existe una solución explicita (aislar $W$ ), de modo que tendremos que aproximarla. Podemos:

Utilizar un descenso de gradiente (explicado anteriormente).
Utilizar Newton-Ransom sobre la primera derivada.

Aplicando Newton-Raphson:

Para una función $y= f(x)$ , el método Newton-Raphson es un método iterativo que aproxima el valor $x$ más cercano al $y=0$ en cada iteración. ¿Como lo hace? Es extremadamente simple se basa en la recta tangente a esa función.

Newton-Raphson https://commons.wikimedia.org/wiki/File:Newton_iteration.svg

La recta tangente era $y = \nabla f(x)\cdot (p - x) + f(x)$ . Donde $x$ es el punto de la función $f(x)$ donde se pretende encontrar la recta tangente y $p$ es la variable independiente de la recta.
Entonces en Newton-Raphson solamene deberemos poner $y=0$ y aislar $p$ .
$0 =\nabla f(x)\cdot (p - x) + f(x)$ $0 = \nabla f(x)p - \nabla f(x)x + f(x)$ $\nabla f(x) p = \nabla f(x)x - f(x)$ $p = \frac{\nabla f(x)x - f(x)}{\nabla f(x)}$ $p = x - \frac{f(x)}{\nabla f(x)}$

¿Como aplicamos Newton-Raphson en nuestro problema?

Como veis el algoritmo nos permite llegar al punto donde $y=0$ , aunque no lo podemos aplicar directamente a nuestro problema… Pero si en la primera derivada ( $\nabla_W$ ). Porque buscamos el 0 en la primera derivada. Quedando de la siguiente forma:
$W_{t+1} = W_{t} - \frac{\nabla W_{t}}{\nabla^2 W_{t}}$ Volvamos a derivar respecto $W$ …
$\nabla^2_W =\sum_i \cancel{\begin{bmatrix} X_i\\ 1 \end{bmatrix} y_i} - \nabla_W \left [ \begin{bmatrix} X_i\\ 1 \end{bmatrix} sigmoid \left ( W\begin{bmatrix} X_i\\ 1 \end{bmatrix} \right ) \right ]$ $\nabla_W \left [ \begin{bmatrix} X_i\\ 1 \end{bmatrix} sigmoid \left ( W\begin{bmatrix} X_i\\ 1 \end{bmatrix} \right ) \right ] = sigmoid \left ( W\begin{bmatrix} X_i\\ 1 \end{bmatrix} \right ) \cdot \left (1 - sigmoid \left ( W\begin{bmatrix} X_i\\ 1 \end{bmatrix} \right ) \right ) \cdot \begin{bmatrix} X_i\\ 1 \end{bmatrix}$ $\nabla^2_W = \sum_i sigmoid \left ( W\begin{bmatrix} X_i \\ 1 \end{bmatrix} \right ) \cdot \left (1 - sigmoid \left ( W\begin{bmatrix} X_i \\ 1 \end{bmatrix} \right ) \right ) \cdot \begin{bmatrix} X_i \\ 1 \end{bmatrix}$ Ya tenemos todas las piezas, veamos como funciona todo esto. ¡Al código!

import pandas as pd
import numpy as np
from sklearn.metrics import confusion_matrix, accuracy_score
# dataset
#https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Coimbra
dataset = pd.read_csv('dataR2.csv')
attributes = list(dataset.columns.values)
dataset = np.array(dataset)
x = dataset[:, :-1]
y = (dataset[:, -1] == 2).astype(x.dtype)
x = np.hstack((x, np.ones(shape=(x.shape[0], 1))))
w_newton = np.zeros(shape=(len(attributes), )) # (2*np.random.rand(*(len(attributes), )) - 1)/(1000**2)
w_descend = np.zeros(shape=(len(attributes), )) #(2*np.random.rand(*(len(attributes), )) - 1)/(1000**2)
logistic_regresor = lambda x, w: 1.0/(1.0 + np.exp(-np.dot(x, w)))
def dd(x, w):
    return - np.sum(np.sum(x*x, axis=1)*(logistic_regresor(x, w)*(1 - logistic_regresor(x, w))), axis=0)
def d(y, x, w):
    return np.sum(np.tile(np.expand_dims(y - logistic_regresor(x, w), axis=-1), [1, x.shape[-1]])*x, axis=0)
#print(w)
for i in range(100000):
    w_newton = w_newton - (d(y, x, w_newton) + 0.00000001)/(dd(x, w_newton) + 0.00000001)
    w_descend = w_descend + 0.001*d(y, x, w_descend)
# newton
print('Newton-Raphson')
true = y
pred = np.round(logistic_regresor(x, w_newton))
print(confusion_matrix(true, pred))
print(accuracy_score(true, pred))
print('Gradient Descend:')
true = y
pred = np.round(logistic_regresor(x, w_descend))
print(confusion_matrix(true, pred))
print(accuracy_score(true, pred))

import pandas as pd

import numpy as np

from sklearn.metrics import confusion_matrix, accuracy_score

# dataset

#https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Coimbra

dataset = pd.read_csv('dataR2.csv')

attributes = list(dataset.columns.values)

dataset = np.array(dataset)

x = dataset[:, :-1]

y = (dataset[:, -1] == 2).astype(x.dtype)

x = np.hstack((x, np.ones(shape=(x.shape[0], 1))))

w_newton = np.zeros(shape=(len(attributes), )) # (2*np.random.rand(*(len(attributes), )) - 1)/(1000**2)

w_descend = np.zeros(shape=(len(attributes), )) #(2*np.random.rand(*(len(attributes), )) - 1)/(1000**2)

logistic_regresor = lambda x, w: 1.0/(1.0 + np.exp(-np.dot(x, w)))

def dd(x, w):

return - np.sum(np.sum(x*x, axis=1)*(logistic_regresor(x, w)*(1 - logistic_regresor(x, w))), axis=0)

def d(y, x, w):

return np.sum(np.tile(np.expand_dims(y - logistic_regresor(x, w), axis=-1), [1, x.shape[-1]])*x, axis=0)

#print(w)

for i in range(100000):

w_newton = w_newton - (d(y, x, w_newton) + 0.00000001)/(dd(x, w_newton) + 0.00000001)

w_descend = w_descend + 0.001*d(y, x, w_descend)

# newton

print('Newton-Raphson')

true = y

pred = np.round(logistic_regresor(x, w_newton))

print(confusion_matrix(true, pred))

print(accuracy_score(true, pred))

print('Gradient Descend:')

true = y

pred = np.round(logistic_regresor(x, w_descend))

print(confusion_matrix(true, pred))

print(accuracy_score(true, pred))

Usando Newton-Raphson obtenemos un accuracy de 77.5%.
Usando Gradient Descend obtenemos un accuracy de 62.9%.
¿Que diferencia hay entre el Gradient Descend y el Newton-Raphson?
Si os fijáis normalmente cuando hacemos gradient descend puro y duro, añadimos el termino que conocemos como learning rate. Podríamos decir que Newton-Raphson estima el mejor learning rate haciendo uso de la segunda derivada $lr = - \frac{1}{\nabla^2 W}$ . Entonces os preguntareis, ¿Porque no se usa Newton-Raphson en redes neuronales? Pues muy simple, no siempre es factible encontrar la segunda derivada tanto por su dificultad como por su coste computacional.
Espero que os haya gustado este post, nos vemos en los próximos.

Referencias:

Tags:

Inferencia Bayesiana: MLE y MAP

Realizar un naive bayes usando MAP:

¡Al código!

Realizar una regresión logística usando MLE:

¿Porque se llama regresor logistico?

Aplicando Newton-Raphson:

¿Como aplicamos Newton-Raphson en nuestro problema?

Referencias:

Tags:

dreamer

Lo último en el blog…

Categorías

Previous PostNeural network from Scratch: Convolution Neural Network: Convolución (Parte 9)

Next PostNeural network from Scratch: Multilayer perceptron 2 (Parte 8)

Leave a Reply Cancel Reply