Introducción:

En este trabajo resolveremos un problema de clasificación de imagenes mediante redes neuronales bayesianas, para esto, explicaremos de forma resumida en que consiste una red neuronal. ### Redes neuronales.

un modelo de red neuronal es un modelo que busca resolver problemas de aprendizaje de maquinas escogiendo funciones de base fijas pero flexible. Estas son funciones de base de base parametrizadas y no lineales cuyos parametros pueden ser ajustados a través de la base de datos durante el entrenamoiento del modelo. También podemos crear una jerarquía dentro de estas funciones de base. Cada nivel de jerarquía es llamada capa de la red. una red consiste en un número fijo de funciones de base, que toman las salidas de la capa anterior como valores de entrada.

De esta forma, Una red neuronal de 1 capa tendría esta forma:

Red neuronal de 1 capas
Red neuronal de 1 capas

Cada uno de los nodos en la imagen representa una función de base calculada de la siguiente forma:

\[a_{i,k}=\sigma(W'a_{i-1}+b_i)\]

Donde \(a_{i,k}\) es el nodo k de la capa i, \(\sigma\) representa la función de base, en este caso una sigmoide, \(a_{i-1}\) es el vector de valores de la capa anterior, \(W\) son los pesos asociados a esociados a esos valores, mientras que \(b_i\) son los sesgos en cada capa.

Para encontrar la red neuronal usualmente lo que necesitamos para resolver nuestro problema es encontrar los pesos adecuados para minimizar la funcion pérdida de nuestra red neuronal, esto puede ser representado de la siguiente forma:

suponiendo que tenemos una muestra de tamaño \(n\) con valores de entrada y salida \(x_i,y_i\) ;\(i=1,..,n\), la probabilidad de que un valor \(x_i\) este bien clasificado será \(\hat{y_i}=f(x_i,w)\). entonces su función de perdida será:

\[L(W)=\sum_{i=1}^n(f(x_i,w)-y(x_i))^2\]

Ahora, si contamos con que existe error en nuestro modelo, entonces: \[\hat{y}_i=f(x_i,w)+\varepsilon_i~~;~~\varepsilon_i\sim(0,\sigma^2)\] Donde\(f(x_i,w)\) es una red neuronal de una capa,de valor: \[f(x_i,w)=\sigma(\sum_{i=1}^H\sigma(\sum_{j=1}^Kx_{ij}w_{ij})w_i)=\sigma(\sum_{i=1}^Ha_iw_i)\] Con \(H\) la cantidad de nodos, \(K\) la cantidad de variables, y \(w\) los pesos

dado que y puede tomar los valores 0 y 1, asumiremos que tiene la siguiente distribución: \[\mathbb{P}(y_i|x_i,w,\sigma^2)=f(x_i,w)^{y_i}(1-f(x_i,w))^{1-y_i}\]

Y su verosimilitud será:

\[L(w|X,Y)=\prod_{i=1}^{n} f(x_i,w)^{y_i}(1-f(x_i,w))^{1-y_i}\]

Su log-verosimilitud será:

\[l(w|X,Y,)=\sum_{i=1}^n y_1log(f(x_i,w))+(1-y_i)log(1-f(x_i,w))\] Ahora tenemos la ventaja de que máximizando la verosimilitud podemos tener predicciones parecidas a las que tendríamos usando la función de pérdida cuadrática, pero para eso necesitamos encontrar valores para sus pesos.

Si llamamos \(alpha_{ij}\) a los pesos \(w_{ij}\) y \(\beta_i\) a los pesos \(w_i\), y les damos distribuciones a priori \(\alpha\sim N(0,I_{H*K}),\beta_i\sim N(0,I_H)\) independientes entre si, entonces la distribución a posteriori \(\pi(\alpha,\beta|X,Y)\) será:

\[\pi(\alpha,\beta|X,Y)\propto L(w|X,Y)\pi(\alpha)\pi(\beta)\] A partir de esta posteriori podemos calcular la siguiente predicción a posteriori:

\[\mathbb{P}(\hat{y}|X,Y)=\int_{\alpha,\beta} \mathbb{P}(\hat{y}|w)\mathbb{P}(w|X,Y)dw=\mathbb{E}_{w|X,Y}(\mathbb{P(\hat{y}|w)})\]

Ejempleficaremos lo anterior mediante el siguiente código:

library(brnn)
## Loading required package: Formula
## Loading required package: truncnorm
datos1 <- read.delim("C:/Users/ernes/OneDrive/Escritorio/bayesiana/pep2bayes/datos1.txt")
xt<-as.matrix(datos1[1:2])
x<-xt[1:40,]
x2<-xt[41:50,]
y<-(datos1[3])
yt<-as.vector(y)
y<-yt$y[1:40]
y<-as.vector(y)
y2<-yt$y[41:50]


out<-brnn(x,y,neurons = 4,Monte_Carlo = T,samples = 2)
## Number of parameters (weights and biases) to estimate: 16 
## Nguyen-Widrow method
## Scaling factor= 0.7246854 
## gamma= 3.8404     alpha= 0.5551   beta= 0.8218
brnn::coef.brnn(out)
##          w[1]          b[1]     beta[1,1]     beta[1,2]          w[2] 
##  0.2076431841  0.0171540664  0.0846015413  0.1649159426  0.0254790585 
##          b[2]     beta[2,1]     beta[2,2]          w[3]          b[3] 
##  0.0006850829  0.0076559618  0.0147951739 -1.2219527314 -0.7385557009 
##     beta[3,1]     beta[3,2]          w[4]          b[4]     beta[4,1] 
## -0.8865401417 -0.6948968571  0.1667803323  0.0153066854  0.0970168380 
##     beta[4,2] 
##  0.1882609756
plot(y,predict(out))

predict(out)
##  [1] 0.2198083 0.7843291 0.8402229 0.1415335 0.8727823 0.9681185 0.9611676
##  [8] 0.5765287 0.9578284 0.4162267 0.8483713 1.0058242 0.6729093 0.7286126
## [15] 0.9591380 1.0193565 0.5089443 0.8822698 0.8880871 1.0838190 0.8512371
## [22] 0.6529083 0.6808755 0.9716940 0.9911848 0.4960066 0.7003748 0.4758619
## [29] 0.8638742 0.5810219 0.9572682 0.7104914 0.4135783 0.7044981 0.6604631
## [36] 0.7288418 0.4720780 0.4687967 0.7429960 0.5077889
cor(out$y,predict(out))
## [1] 0.5598802
sens<-ifelse(predict(out)<=0.50,0,1)
table(y,sens)
##    sens
## y    0  1
##   0  6  5
##   1  2 27
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
par(pty="s")
roc(y,predict(out),plot=T,print.auc=T)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

## 
## Call:
## roc.default(response = y, predictor = predict(out), plot = T,     print.auc = T)
## 
## Data: predict(out) in 11 controls (y 0) < 29 cases (y 1).
## Area under the curve: 0.8245

```