Bayesian Neural Network¶
Autor do notebook: Felipe Maia Polo
Nesta implementação seguiremos em parte a proposta apresentada no livro: Bishop, C. M. (2006). Pattern recognition and machine learning. Springer Science+ Business Media.. Assim, quando citar o autor do livro, estarei me referindo ao Bishop.
Abrindo as bibliotecas necessárias
library(MASS)
library(e1071)
library(ggplot2)
library(sigmoid)
Gerando dados artificiais
n0 <- 1200
n1 <- 1200
X0 <- rbind(mvrnorm(n = n0/3, mu=c(1,2), matrix(c(2,-1,-1,1),2,2)),
mvrnorm(n = n0/3, mu=c(-3,-1), matrix(c(2,-1,-1,1),2,2)),
mvrnorm(n = n0/3, mu=c(-2,2), matrix(c(2,1,1,1),2,2)))
X1 <- mvrnorm(n = n1, mu=c(-0.5,0), matrix(c(1.5,0,0,1),2,2))
y0 <- rep(0,n0)
y1 <- rep(1,n1)
X <- rbind(X0,X1)
y <- c(y0,y1)
Plotando os dados:
data <- as.data.frame(cbind(X, y))
names(data) <- c('x1','x2','y')
qplot(x1, x2, colour = y, data = data)
No modelo proposto pelo autor, temos a seguinte função de verossimilhança:
\begin{align} p(\mathbf{y}|\mathbf{x},\mathbf{w})=\prod_{i=1}^{n}p(y_i|x_i,\mathbf{w})=\prod_{i=1}^{n}h(x_i,\mathbf{w})^{y_i}(1-h(x_i,\mathbf{w}))^{(1-y_i)} \end{align}Sendo que $h(.)$ é uma rede neural com uma camada oculta da seguinte forma:
\begin{align} h(\mathbf{x_i},\mathbf{w})=\sigma \Big (\mathbf{z_{2i}}^\top \mathbf{w_2}\Big)=\sigma_i \end{align}\begin{align} \mathbf{z_{2i}}=\begin{bmatrix}1\\ f\big(\mathbf{z_{1i}}^\top \mathbf{W_1} \big)^\top \end{bmatrix} \end{align}\begin{align} \mathbf{z_{1i}} = \begin{bmatrix}1\\ \mathbf{x_i}\end{bmatrix} \end{align}Onde $f(.)$ é uma função não linear aplicada coordenada a coordenada. Nesse caso, vamos assumir que $f(.)=\text{max}(.,0)=\text{ReLU}(.)$, Veja que internalizamos o intercepto (bias) nos próprios vetores, o que facilitará as contas daqui para frente. Acima temos que as equações foram escritas para o indivÃduo $i$, no entanto, codificaremos as funções já pensando na base inteira, o que deixará nosso algoritmo mais eficiente - pensaremos em $X$ como uma matriz em que cada linha diz respeito a um indivÃduo. Definindo as funções:
sig <- function(z){return(1/(1+exp(-z)))}
f <- function(z){return(apply(z,2,relu))}
itc <- function(Z){return(cbind(rep(1,NROW(Z)),Z))}
h <- function(Z,V1,v2){return(sig(itc(f(itc(Z)%*%V1))%*%v2))}
Definimos então a distribuição a priori dos pesos da rede neural da seguinte maneira
\begin{align} \mathbf{W} |\alpha\sim N(0,\alpha^{-1} I_k) \end{align}Sendo $\alpha$ um hiperparâmetro de precisão. Como queremos fixar uma prior pouco informativa, vamos assumir que $\alpha=\frac{1}{100}$. Abaixo vamos simular previsões feitas pela rede se nenhum aprendizado:
#Tamanho da amostra
n <- NROW(X)
#Número de variáveis
d <- NCOL(X)
#Número de neurônios na camada oculta
k1=5
#Número de neurônios na camada de saÃda
k2=1
#Alpha
alpha=1
Predizendo a probabilidade para os 5 primeiros indivÃduos
iter <- 1000
simul <- c()
mu1 <- rep(0,(d+1)*k1)
Sigma1 <- alpha*diag((d+1)*k1)
mu2 <- rep(0,(k1+1)*k2)
Sigma2 <- alpha*diag((k1+1)*k2)
for(i in 1:iter){
#Amostrando W1
W1 <- matrix(mvrnorm(1, mu1, Sigma1), nrow = (d+1), byrow = F) #cada coluna diz respeito a uma variável de entrada
#Amostrando w2
w2 <- mvrnorm(1, mu2, Sigma2)
#Rbind
simul <- cbind(simul, h(X,W1,w2))
}
Visualizando $\hat{P}(Y_1=1|\mathbf{x_1},\mathbf{w},\alpha)$, sendo $\mathbf{w}$ realizado de acordo com sua priori:
dim(simul)
hist(simul[1,],main='Histograma da probabilidade de sucesso para o indivÃduo 1')
Voltando à teoria, temos a seguinte proporcionalidade pelo teorema de Bayes:
\begin{align} p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha) \propto p(\mathbf{y}|\mathbf{x},\mathbf{w},\alpha)p(\mathbf{w}|\alpha) \end{align}Tirando o logaritmo da posteriori, temos que:
\begin{align} \text{ln} p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha)&=\text{ln}p(\mathbf{y}|\mathbf{x},\mathbf{w},\alpha)+\text{ln}p(\mathbf{w}|\alpha) +\text{cte} \\ &=\text{ln}\prod_{i=1}^{n}p(y_i|x_i,\mathbf{w})+\text{ln}\prod_{j}p(\mathbf{w}_j|\alpha) +\text{cte} \\ &=\text{ln}\prod_{i=1}^{n}h(x_i,\mathbf{w})^{y_i}(1-h(x_i,\mathbf{w}))^{(1-y_i)}+\text{ln}\prod_{j}\sqrt{\frac{\alpha}{2 \pi}} \text{exp}\Big(-\frac{\alpha(w_j^2)}{2}\Big) +\text{cte} \\ &=\sum_{i=1}^{n}{y_i}\text{ln}h(x_i,\mathbf{w})+(1-y_i)\text{ln}(1-h(x_i,\mathbf{w}))-\frac{\alpha}{2}\left \| \mathbf{w} \right \|^2 +\text{cte} \\ \end{align}Sendo que $\circ$ é o produto elemento a elemento de matrizes. Como $p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha)$ assume uma forma muito complicada, aproximaremos essa distribuição da maneira que segue. Primeiramente, aproximaremos o logaritmo da posteriori utilizando o polinômio de Taylor de segunda ordem centrado em $\mathbf{w}_{\text{MAP}}=\underset{\mathbf{w}}{\text{argmax}}~p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha)=\underset{\mathbf{w}}{\text{argmax}}~\text{ln}p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha)$:
\begin{align} \text{ln}p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha) &\approx \text{ln}p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)+ (\mathbf{w}-\mathbf{w}_{\text{MAP}})^\top\nabla \text{ln}p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)-\frac{1}{2}(\mathbf{w}-\mathbf{w}_{\text{MAP}})^\top A(\mathbf{w}-\mathbf{w}_{\text{MAP}}) \\ &= p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)-\frac{1}{2}(\mathbf{w}-\mathbf{w}_{\text{MAP}})^\top A(\mathbf{w}-\mathbf{w}_{\text{MAP}}) \end{align}Sendo que $A=-\frac{\partial ^2 }{\partial \mathbf{w}^2}\text{ln}p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)$ é o negativo da matriz Hessiana do log da posteriori no ponto $\mathbf{w}=\mathbf{w}_{\text{MAP}}$. Tirando o exponencial dos dois lados, temos:
\begin{align} p(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha) &\approx p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)*\text{exp}\Big(-\frac{1}{2}(\mathbf{w}-\mathbf{w}_{\text{MAP}})^\top A(\mathbf{w}-\mathbf{w}_{\text{MAP}})\Big) \\ &= c(\mathbf{w}_{\text{MAP}}) *\text{exp}\Big(-\frac{1}{2}(\mathbf{w}-\mathbf{w}_{\text{MAP}})^\top A(\mathbf{w}-\mathbf{w}_{\text{MAP}})\Big) \end{align}É possÃvel visualizar o núcleo da distirbuição normal na fórmula acima e, assim, é natural pensar que uma possÃvel aproximação para a posteriori é a seguinte distribuição:
\begin{align} q(\mathbf{w}|\mathbf{y},\mathbf{x},\alpha)=N(\mathbf{w}_{\text{MAP}},\mathbf{A}^{-1}) \end{align}Primeiramente, realizaremos o cálculo para $\mathbf{w}_{\text{MAP}}$ e depois para a matriz $A$. Realizaremos o cálculo de $\mathbf{w}_{\text{MAP}}$ achando o argumento que maximiza
\begin{align} g(\mathbf{w})=\sum_{i=1}^{n}{y_i}\text{ln}h(x_i,\mathbf{w})+(1-y_i)\text{ln}(1-h(x_i,\mathbf{w}))-\frac{\alpha}{2}\left \| \mathbf{w} \right \|^2 \end{align}Para tal feito, teremos que recorrer ao algoritmo do Stochastic Grandient Descent e backpropagation. Veja que $\mathbf{w}$ pode ser particionado em $\mathbf{w_1}$ e em $\mathbf{w_2}$ e que $\nabla_{\mathbf{w}_1} g(\mathbf{w})$ depende de $\mathbf{w_2}$ enquanto $\nabla_{\mathbf{w}_2} g(\mathbf{w})$ não depende de $\mathbf{w}_1$. Calcularemos então o segundo gradiente em primeiro lugar:
\begin{align} \nabla_{\mathbf{w}_2} g(\mathbf{w})&= \sum_{i=1}^{n}\Bigg[{y_i}\nabla_{\mathbf{w}_2}\text{ln}h(x_i,\mathbf{w})+(1-y_i)\nabla_{\mathbf{w}_2}\text{ln}(1-h(x_i,\mathbf{w}))\Bigg]-\frac{\alpha}{2}\nabla_{\mathbf{w}_2}\left \| \mathbf{w} \right \|^2 \\ &= \sum_{i=1}^{n}\Bigg[{y_i}\nabla_{\mathbf{w}_2}\text{ln}\sigma \big ( \mathbf{z_{2i}}^\top \mathbf{w_2} \big)+(1-y_i)\nabla_{\mathbf{w}_2}\text{ln}\Big(1-\sigma \big (\mathbf{z_{2i}}^\top \mathbf{w_2} \big)\Big)\Bigg]-\alpha \mathbf{w}_2, \text{ adotando } \sigma_i=\sigma \big (\mathbf{z_{2i}}^\top \mathbf{w_2} \big) \\ &= \sum_{i=1}^{n}\Bigg[y_i(1-\sigma_i)\mathbf{z}_{2i}+(1-y_i)\sigma_i\mathbf{z}_{2i}\Bigg]-\alpha \mathbf{w}_2 \\ &= \sum_{i=1}^{n}\Bigg[(y_i-\sigma_i)\mathbf{z}_{2i}\Bigg]-\alpha \mathbf{w}_2 \\ &=\begin{bmatrix} \mathbf{z}_{2,1} & ...& \mathbf{z}_{2,n}\end{bmatrix}(\mathbf{y}-\boldsymbol{\sigma})-\alpha \mathbf{w}_2\\ &=\mathbf{Z}_2(\mathbf{y}-\boldsymbol{\sigma})-\alpha \mathbf{w}_2 \\ \end{align}Antes de partirmos para o cálculo em si, é necessário frizar que $W_1$ é uma matriz bidimensional (de ordem $d+1$ X $k_1$). Para facilitar o cálculo vamos ver a matriz como uma concatenação de vetores da seguinte forma $\mathbf{W}_1=\begin{bmatrix}\mathbf{w}_{1,1} & ... & \mathbf{w}_{1,k_1}\end{bmatrix}$, sendo que representaremos o jacobiano da seguinte maneira $\frac{\partial }{\partial \mathbf{W}_1} g(\mathbf{w})=\begin{bmatrix}\frac{\partial }{\partial \mathbf{w}_{1,1}} g(\mathbf{w}) & ... & \frac{\partial }{\partial \mathbf{w}_{1,k_1}} g(\mathbf{w})\end{bmatrix}$. Botando a mão na massa:
\begin{align} \nabla_{\mathbf{w}_{1,t}} g(\mathbf{w})&= \nabla_{\mathbf{w}_{1,t}} g(\mathbf{w})= \nabla_{\mathbf{w}_{1,t}} \Bigg[\sum_{i=1}^{n}{y_i}\text{ln}h(x_i,\mathbf{w})+(1-y_i)\text{ln}(1-h(x_i,\mathbf{w}))-\frac{\alpha}{2}\left \| \mathbf{w} \right \|^2\Bigg] \\ &= \nabla_{\mathbf{w}_{1,t}} \Bigg[\sum_{i=1}^{n}{y_i}\text{ln}h(x_i,\mathbf{w})+(1-y_i)\text{ln}(1-h(x_i,\mathbf{w}))\Bigg]-\frac{\alpha}{2}\nabla_{\mathbf{w}_{1,t}} \left \| \mathbf{w} \right \|^2 \\ &= \sum_{i=1}^{n}\Bigg[{y_i} \nabla_{\mathbf{w}_{1,t}} \text{ln}h(x_i,\mathbf{w})+(1-y_i)\nabla_{\mathbf{w}_{1,t}}\text{ln}(1-h(x_i,\mathbf{w}))\Bigg] -\alpha \mathbf{w}_{1,t} \\ &=\sum_{i=1}^{n}\Bigg[{y_i} \frac{\partial \text{ln}h(x_i,\mathbf{w})}{\partial \sigma_i}\frac{\partial \sigma_i}{\partial \mathbf{w}_{1,t}} +(1-y_i) \frac{\partial \text{ln}\big(1-h(x_i,\mathbf{w})\big)}{\partial \sigma_i}\frac{\partial \sigma_i}{\partial \mathbf{w}_{1,t}}\Bigg] -\alpha \mathbf{w}_{1,t} , \text{ adotando } \sigma_i=\sigma \big (\mathbf{w_2}^T \mathbf{z_{2i}} \big) \\ &=\sum_{i=1}^{n}\Bigg[{y_i} (1-\sigma_i) \frac{\partial \mathbf{z}_{2i}^\top\mathbf{w}_2}{\partial \mathbf{w}_{1,t}} -(1-y_i) \sigma_i \frac{\partial \mathbf{z}_{2i}^\top\mathbf{w}_2}{\partial \mathbf{w}_{1,t}}\Bigg] -\alpha \mathbf{w}_{1,t} \\ &=\sum_{i=1}^{n}\Bigg[(y_i-\sigma_i) \frac{\partial \mathbf{z}_{2i}^\top\mathbf{w}_2}{\partial \mathbf{w}_{1,t}} \Bigg] -\alpha \mathbf{w}_{1,t}=\bigstar \\ \end{align}Vamos calcular a derivada que sobrou antes de continuarmos a resolução
\begin{align} \frac{\partial \mathbf{z}_{2i}^\top\mathbf{w}_2}{\partial \mathbf{w}_{1,t}}=\frac{\partial }{\partial \mathbf{w}_{1,t}} \Big[w_{2,0} + \sum^{k_1}_{j=1}w_{2,j}f(\mathbf{z}_{1i}^\top \mathbf{w}_{1,j}) \Big]=w_{2,t}\frac{\partial }{\partial \mathbf{w}_{1,t}}f(\mathbf{z}_{1i}^\top \mathbf{w}_{1,t})=w_{2,t}\delta_{it}\mathbf{z}_{1i} \end{align}Sendo que $\delta_{it}=\mathbb{I}(\mathbf{z}_{1i}^\top \mathbf{w}_{1,t}>0)$ é uma indicadora se o argumento é positivo. Voltando:
\begin{align} \bigstar&=w_{2,t} \Bigg[\sum_{i=1}^{n}(y_i-\sigma_i) \delta_{it} \mathbf{z}_{1i}\Bigg] -\alpha \mathbf{w}_{1,t} \\ &=w_{2,t}\begin{bmatrix} \mathbf{z}_{1,1} & ...& \mathbf{z}_{1,n}\end{bmatrix}\big[(\mathbf{y}-\boldsymbol{\sigma}) \circ \boldsymbol{\delta_t}\big]-\alpha \mathbf{w}_{1,t}\\ &=w_{2,t}\mathbf{Z}_1\big[(\mathbf{y}-\boldsymbol{\sigma}) \circ \boldsymbol{\delta_t}\big]-\alpha \mathbf{w}_{1,t} \\ \end{align}Logo,
\begin{align} \frac{\partial }{\partial \mathbf{W}_1} g(\mathbf{w})&=\begin{bmatrix}w_{2,1}\mathbf{Z}_1\big[(\mathbf{y}-\boldsymbol{\sigma}) \circ \boldsymbol{\delta_1}\big]-\alpha \mathbf{w}_{1,1} & ... & w_{2,k_1}\mathbf{Z}_1\big[(\mathbf{y}-\boldsymbol{\sigma}) \circ \boldsymbol{\delta_{k_1}}\big]-\alpha \mathbf{w}_{1,k_1}\end{bmatrix} \\ &=\mathbf{Z}_1 \Big[\begin{bmatrix}w_{2,1}(\mathbf{y}-\boldsymbol{\sigma})& ... & w_{2,k_1}(\mathbf{y}-\boldsymbol{\sigma})\end{bmatrix}\circ \boldsymbol{\delta} \Big]-\alpha \mathbf{W}_{1} \\ &=\mathbf{Z}_1 \Big[(\mathbf{y}-\boldsymbol{\sigma})\mathbf{w}^\top_2\circ \boldsymbol{\delta} \Big]-\alpha \mathbf{W}_{1}, ~ \boldsymbol{\delta}=(\delta_{it}) \text{ para } i \in \{1,...,n\} \text{ e } t \in \{1,...,k_1\}\\ \end{align}Definindo função que queremos maximizar:
g <- function(w,Z,V1,v2,alpha){return(sum(log(cbind(1-h(Z,V1,v2),h(Z,V1,v2)))*cbind(0^w,w))-alpha*(sum(v2^2)+sum(V1^2)))}
Definindo os nossos gradientes:
grad2<-function(y,X,W1,w2,alpha){return(t(itc(f(itc(X)%*%W1)))%*%(y-h(X,W1,w2))-alpha*w2)}
grad1<-function(y,X,W1,w2,alpha){
delta=1*(itc(X)%*%W1>0)
return(t(itc(X))%*%((y-h(X,W1,w2))%*%t(w2[2:(k1+1)])*delta)-alpha*W1)
}
Apesar de acharmos que as contas que fizemos estão certas, é sempre bom fazer uma comparação de um resultado numérico e do nosso resultado analÃtico. Fixando $W_1$ e $w_2$, vamo fazer as razões dos gradientes numéricos e analÃticos. Primeiramente vamos definir as funções que calculam os gradientes numéricos:
grad2_check<-function(){
d=10^(-10)
grad2_num=rep(0,length(w2))
for(i in 1:length(w2)){
u=rep(0,length(w2))
u[i]=d
grad2_num[i]=(g(y,X,W1,w2+u,alpha)-g(y,X,W1,w2-u,alpha))/(2*d)
}
return(grad2_num)
}
grad1_check<-function(){
d=10^(-10)
grad1_num=matrix(rep(0,prod(dim(W1))),NROW(W1),NCOL(W1))
for(i in 1:NROW(W1)){
for(j in 1:NCOL(W1)){
u=matrix(rep(0,prod(dim(W1))),NROW(W1),NCOL(W1))
u[i,j]=d
grad1_num[i,j]=(g(y,X,W1+u,w2,alpha)-g(y,X,W1-u,w2,alpha))/(2*d)
}
}
return(grad1_num)
}
Vamos efetuar as razões agora:
grad2_check()/grad2(y,X,W1,w2,alpha)
grad1_check()/grad1(y,X,W1,w2,alpha)
É possÃvel ver que as razões todas se aproximam de 1, o que nos dá indÃcios de que nossas contas estão corretas. Vamos avaliar agora a função que queremos maximizar pesos iniciais aleatórios:
scale<-.1
#inicializando W1
mu1 <- rep(0,(d+1)*k1)
Sigma1 <- scale*diag((d+1)*k1)
W1 <- matrix(mvrnorm(1, mu1, Sigma1), nrow = (d+1), byrow = F) #cada coluna diz respeito a uma variável de entrada
#inicializando w2
mu2 <- rep(0,(k1+1)*k2)
Sigma2 <- scale*diag((k1+1)*k2)
w2 <- mvrnorm(1, mu2, Sigma2)
#Avaliando g
g(y,X,W1,w2,alpha)
Vamos agora aplicar o mini-batch gradient ascent para maximizar nossa função $g$:
lr=0.001
epochs=500
g_hist=c()
batch=250
mom=.1
iter=round(nrow(data)/batch)
for(i in 1:epochs){
#shuffle dataframe
data <- data[sample(nrow(data)),]
X <- as.matrix(data[,-NCOL(data)])
y <- as.matrix(data[,NCOL(data)])
for(j in 1:iter){
#Index
index <- ((j-1)*batch+1):(min(j*batch,nrow(data)))
#Update
w2 <- w2 + lr*grad2(y[index],X[index,],W1,w2,alpha)
W1 <- W1 + lr*grad1(y[index],X[index,],W1,w2,alpha)
}
#Loss
g_hist <- c(g_hist,g(y,X,W1,w2,alpha))
}
Visualizando o histórico de treino:
plot(g_hist)
Avaliando nossa função novamente, mas agora com os pesos treinado e em nossa base de testes:
#Criando base para testes
X0 <- rbind(mvrnorm(n = n0/3, mu=c(1,2), matrix(c(2,-1,-1,1),2,2)),
mvrnorm(n = n0/3, mu=c(-3,-1), matrix(c(2,-1,-1,1),2,2)),
mvrnorm(n = n0/3, mu=c(-2,2), matrix(c(2,1,1,1),2,2)))
X1 <- mvrnorm(n = n1, mu=c(-0.5,0), matrix(c(1.5,0,0,1),2,2))
y0 <- rep(0,n0)
y1 <- rep(1,n1)
X_test <- rbind(X0,X1)
y_test <- c(y0,y1)
#Avaliando g
g(y_test,X_test,W1,w2,alpha)
Observe que a função $g$ nos retorna resultados muitos parecidos tanto para o conjunto de treino quanto para o conjunto de teste, o que nos dá indÃcios que nossa rede 'overfita' pouco ou não 'overfita'. Vendo as probabilidades preditas pela nossa rede quando fixamos os pesos em $\mathbf{w}_{\text{MAP}}$ (uma solução pontual para o problema):
data <- as.data.frame(cbind(X_test, y_test, h(X_test,W1,w2)))
names(data) <- c('x1','x2','y','prob')
qplot(x1, x2, colour = prob, data = data)
Agora que encontramos $\mathbf{w}_{\text{MAP}}$ nos resta o cálculo de $A=-\frac{\partial ^2 }{\partial \mathbf{w}^2}\text{ln}p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)$. Veja que o que fizemos até o momento é o cálculo de $\frac{\partial }{\partial \mathbf{w}}\text{ln}p(\mathbf{w}_{\text{MAP}}|\mathbf{y},\mathbf{x},\alpha)$, que é o gradiente. Como teremos que realizar o cálculo de $A$ somente uma vez e a nossa rede neural não é muito grande, aproximaremos os valores da matriz numericamente e depois a invertaremos para obter a covariância da posteriori aproximada. Para o cálculo numérico utilizaremos as funções grad1 e grad2 que já foram codificadas - variando elementos de $\mathbf{W}_1$ e de $\mathbf{w}_2$ de forma 'infinitesimal' guardaremos as taxas de variação nos outputs de grad1 e grad2 como as nossas segundas derivadas (veja que a matriz $A$ deve ser simétrica dada a natureza do nosso problema).
reshape_W1_1<-function(W1){matrix(W1)} # (d+1) X k1 >>> (d+1)*k1 X 1
reshape_W1_2<-function(W1){matrix(W1,d+1,k1)} # (d+1)*k1 X 1 >>> (d+1) X k1
grad1_2<-function(y,X,W1,w2,alpha){return(reshape_W1_1(grad1(y,X,reshape_W1_2(W1),w2,alpha)))}
#d=2
#Variando W1
#Bloco 1 (up left)
del=10^(-10)
r_W1=reshape_W1_1(W1)
block1<-c()
for(i in 1:(prod(dim(W1)))){
u=rep(0,prod(dim(W1)))
u[i]=del
block1<-cbind(block1, (grad1_2(y,X,r_W1+u,w2,alpha)-grad1_2(y,X,r_W1-u,w2,alpha))/(2*del))
}
#Variando w2
#Bloco 3 (bottom left)
del=10^(-10)
r_W1=reshape_W1_1(W1)
block2<-c()
for(i in 1:length(w2)){
u=rep(0,length(w2))
u[i]=del
block2<-cbind(block2,(grad1_2(y,X,r_W1,w2+u,alpha)-grad1_2(y,X,r_W1,w2-u,alpha))/(2*del))
}
#Variando w2
#Bloco 3 (bottom left)
del=10^(-10)
block3<-c()
for(i in 1:length(w2)){
u=rep(0,length(w2))
u[i]=del
block3<-cbind(block3, (grad2(y,X,W1,w2+u,alpha)-grad2(y,X,W1,w2-u,alpha))/(2*del))
}
Matriz A:
block2<-0*block2
A <- -cbind(rbind(block1,t(block2)),rbind(block2,block3))
dim(A)
Calculando a covariância invertendo A:
cov<-solve(A)
for(i in 1:NROW(cov)){
for(j in 1:i){
m<-mean(cov[i,j],cov[j,i])
cov[i,j]<-m
cov[j,i]<-m
}
}
Checando os autovalores de Cov (todos devem ser não negativos, pois se trata de uma matrix de covariância):
e <- eigen(cov)
e$values
Concatenando nossos pesos para obter o $\mathbf{w}_{\text{MAP}}$:
w_map<-rbind(reshape_W1_1(W1),w2)
length(w_map)
dim(cov)
Gerando uma instância de $\mathbf{w}$ como exemplo:
mvrnorm(1, w_map, cov)
Assim como já fizemos, amostraremos da distribuição dos pesos (agora da posteriori) e vamos obter 1000 predições para cada indivÃduo que dependem de cada um dos valores amostrados para $\mathbf{w}$:
iter <- 1000
simul <- c()
mu1 <- rep(0,(d+1)*k1)
Sigma1 <- alpha*diag((d+1)*k1)
for(i in 1:iter){
#Amostrando w
w <- mvrnorm(1, w_map, cov)
W1 <- reshape_W1_2(w[1:((d+1)*k1)])
w2 <- w[-(1:((d+1)*k1))]
#Rbind
simul <- cbind(simul, h(X_test,W1,w2))
}
Agora vamos ver as distirbuições de predições para dois indivÃduos, começando pelo indivÃduo 3:
hist(simul[3,],main='Histograma da probabilidade de sucesso para o indivÃduo 3')
É possÃvel ver que a rede está muito certa que este indivÃduo é da classe 0. Olhando agora para o indivÃduo 8, podemos ver que a rede está menos certa a respeito de que classe pertence o indivÃduo:
hist(simul[8,], main='Histograma da probabilidade de sucesso para o indivÃduo 8')
O fato de a rede estar mais ou menos certa de que classes os indivÃduos são tem muito a ver com a localização do indivÃduo. Veja que o indivÃduo 8 está em uma zona de divisão entre as classes, enquando o indivÃduo 3 está em uma zona que é predominantemente da classe 0:
X_test[3,]
X_test[8,]
data <- as.data.frame(cbind(X_test, y_test))
names(data) <- c('x1','x2','y')
qplot(x1, x2, colour = y, data = data)
