Exploração Uniforme

Nas primeiras N rodadas é escolhido um braço aleatoriamente e é coletado sua recompensa com base na distribuição de probabilidade desse braço. A partir das demais rodadas é escolhido o braço com maior probabilidade de recompensa.

#Definição do Epsilon
policy <- EpsilonFirstPolicy$new(epsilon = 0.25, N = 100)


#Foi definido que a simulação será realizada com 3 "braços". Ambos com distribuição Bernoulli com Média/Esperança igual a respectivamente: 0.4, 0.5, 0.3 
bandit             <- BasicBernoulliBandit$new(weights = c(0.5, 0.2, 0.1))

#Definição do Método 
agent              <- Aef_agent <- Agent$new(policy, bandit, "Exploration")

# Rodar a simulação
simulator          <- Simulator$new(agents      = agent,
                                    horizon     = 100, # Quantas vezes será realizado o experimento
                                    set_seed    = 8, # Semente para sempre vir os mesmos resultados
                                    save_theta  = TRUE,
                                    simulations = 100, # Quantas simulações realizadas
                                    do_parallel = TRUE
)

history <- simulator$run()


# Acessando Informações da Simulação:
summary(history)

## 
## Agents:
## 
##   Exploration
## 
## Cumulative regret:
## 
##        agent   t sims cum_regret cum_regret_var cum_regret_sd
##  Exploration 100  100        8.7       64.13131      8.008203
## 
## 
## Cumulative reward:
## 
##        agent   t sims cum_reward cum_reward_var cum_reward_sd
##  Exploration 100  100      41.16       85.10545      9.225262
## 
## 
## Cumulative reward rate:
## 
##        agent   t sims cur_reward cur_reward_var cur_reward_sd
##  Exploration 100  100     0.4116      0.8510545    0.09225262

# Gráfico de Escolha dos Braços
plot(history, type = "arms")

#Gráfico de avaliação da probabilidade média da recompensa
plot(history, type = "average", regret = FALSE)

# Grafico de avaliação da probabilidade acumulativa
plot(history, type = "cumulative", regret = FALSE, rate = TRUE)

Epsilon Greedy

Defina-se um Epsilon, dado que 0<\(\epsilon\)<1. Após isso, sorteia-se um número aleatório “S” entre 0 e 1. Quando:

• S>=\(\epsilon\) : Escolha o braço com a melhor recompensa esperada atual.

• S<\(\epsilon\): Escolher outro braço aleatoriamente.

Dessa forma, temos que a probabilidade de se escolher o maior braço equivale a 1 - \(\epsilon\) e, consequentemente, a probabilidade de se escolher um braço aleatóriamente é dada como \(\epsilon\).

#Definição do Epsilon
policy             <- EpsilonGreedyPolicy$new(epsilon = 0.1) 

#Foi definido que a simulação será realizada com 3 "braços". Ambos com distribuição Bernoulli com Média/Esperança igual a respectivamente: 0.4, 0.5, 0.3 
bandit             <- BasicBernoulliBandit$new(weights = c(0.5, 0.2, 0.1))

#Definição do Método - E-Greedy
agent              <- Agent$new(policy,bandit, "EG")

# Rodar a simulação
simulator          <- Simulator$new(agents      = agent,
                                    horizon     = 100, # Quantas vezes será realizado o experimento
                                    set_seed    = 8, # Semente para sempre vir os mesmos resultados
                                    save_theta  = TRUE,
                                    simulations = 100 # Quantas simulações realizadas
                                    )

history <- simulator$run()


# Acessando Informações da Simulação:
summary(history)

## 
## Agents:
## 
##   EG
## 
## Cumulative regret:
## 
##  agent   t sims cum_regret cum_regret_var cum_regret_sd
##     EG 100  100        8.6       80.94949      8.997194
## 
## 
## Cumulative reward:
## 
##  agent   t sims cum_reward cum_reward_var cum_reward_sd
##     EG 100  100      41.31       113.6302      10.65975
## 
## 
## Cumulative reward rate:
## 
##  agent   t sims cur_reward cur_reward_var cur_reward_sd
##     EG 100  100     0.4131       1.136302     0.1065975

# Gráfico de Escolha dos Braços
plot(history, type = "arms")

#Gráfico de avaliação da probabilidade média da recompensa
plot(history, type = "average", regret = FALSE)

# Grafico de avaliação da probabilidade acumulativa
plot(history, type = "cumulative", regret = FALSE, rate = TRUE)

UCDB1

Princípio do otimismo diante da incerteza: braços que foram pouco explorados são escolhidos com maior frequência para que as incertezas com relação as suas recompensas sejam reduzidas.

Esse algoritmo é considerado o mais simples de um grupo de outros algoritmos cujos os objetivos são os mesmos, utilizando o princípio do otimismo diante da incerteza.

Considere em um determinado intervalo de tempo que a função de recompensa de cada braço pode ser percebida como uma estimativa pontual com base na taxa média de recompensa observada. Extraindo a intuição dos intervalos de confiança, para cada estimativa pontual, também podemos incorporar alguma forma de limite de incerteza em torno da estimativa pontual. Nesse sentido, temos o limite inferior e o limite superior para cada braço.

O algoritmo UCB é apropriadamente nomeado porque estamos preocupados apenas com o limite superior, já que estamos tentando encontrar o braço com a maior taxa de recompensa. No contexto da UCB, em cada rodada de n testes, a recompensa UCB de todas as armas é representada pela seguinte fórmula:

\(\mu_{i} + \sqrt{\frac{2ln\left ( n \right )}{n_{i}}}\)

Para cada passo de tempo, escolhemos o braço que possui o maior valor de UCB (média estimada mais seu limite de confiança) e tentamos observar a recompensa apenas para esse braço. A cada iteração, atualizamos a média de retorno de cada braço por meio da equação:

A média de cada braço \(\mu_{i}\) aumenta com recompensa positiva (1) e diminui sem recompensa (0). O limite de confiança é atualizado por \(n\) (número de etapas de tempo gerais) e \(n_{i}\) (número de vezes que escolhemos esse braço). No numerador, observamos a complexidade logarítmica do tempo onde a fronteira aumenta \(ln(n)\) com a passagem do tempo. Com o denominador, podemos ver que o limite diminui em complexidade de tempo linear com \(n_{i}\). O algoritmo, portanto, explora os braços que podem não ter um bom desempenho sobre os quais temos alta incerteza (portanto, um alto UCB), enquanto explora os braços que têm retornos médios superiores, mas com alta certeza.

#Definição da Regra/Método - UCB1
policy             <- UCB1Policy$new() 

#Foi definido que a simulação será realizada com 3 "braços". Ambos com distribuição Bernoulli com Média/Esperança igual a respectivamente: 0.4, 0.5, 0.3 
bandit             <- BasicBernoulliBandit$new(weights = c(0.5, 0.2, 0.1))

#Definição do Agente
agent              <- Agent$new(policy,bandit, "UCB1")

# Rodar a simulação
simulator          <- Simulator$new(agents      = agent,
                                    horizon     = 100, # Quantas vezes será realizado o experimento
                                    set_seed    = 8, # Semente para sempre vir os mesmos resultados
                                    save_theta  = TRUE,
                                    simulations = 100 # Quantas simulações realizadas
)

history <- simulator$run()


# Acessando Informações da Simulação:
summary(history)

## 
## Agents:
## 
##   UCB1
## 
## Cumulative regret:
## 
##  agent   t sims cum_regret cum_regret_var cum_regret_sd
##   UCB1 100  100      13.12       13.62182      3.690775
## 
## 
## Cumulative reward:
## 
##  agent   t sims cum_reward cum_reward_var cum_reward_sd
##   UCB1 100  100      36.39       32.48273      5.699362
## 
## 
## Cumulative reward rate:
## 
##  agent   t sims cur_reward cur_reward_var cur_reward_sd
##   UCB1 100  100     0.3639      0.3248273    0.05699362

# Gráfico de Escolha dos Braços
plot(history, type = "arms")

#Gráfico de avaliação da probabilidade média da recompensa
plot(history, type = "average", regret = FALSE)

# Grafico de avaliação da probabilidade acumulativa
plot(history, type = "cumulative", regret = FALSE, rate = TRUE)

Thompson Sampling

O Thompson Sampling (TS) parte de uma ideia de modelar uma distribuição de probabilidade a partir das recompensas e usar esse modelo para escolher o braço com a maior esperança de recompensa.

Uma distribuição muito utilziada é a distribuição Beta. Logo, é admitido que o \(\theta\) de cada braço possui uma distribuição beta com alfa e beta iguais a 1. Dessa forma definimos a distribuição a priori dada como:

\(\theta_{i} \sim (1,1)\)

Através dessa definição e da distribuição de cada braço já conhecida, podemos estimar a distribuição posteriori.

Algoritmo:

• Inicializar o alpha e beta igual a 1, o que inicia uma distribuição com probabilidade de recompensa de 50%.
• A cada rodada, amostramos uma recompensa através da distribuição que estamos criando de cada braço e selecionamos a que tem a maior recompensa esperada.
• Após receber a verdadeira recompensa atualizamos os valores de alpha e beta para que a próxima estimativa esteja mais correta.

#Definição do Método - TS
policy             <- ThompsonSamplingPolicy$new()

#Foi definido que a simulação será realizada com 3 "braços". Ambos com distribuição Bernoulli com Média/Esperança igual a respectivamente: 0.4, 0.5, 0.3 
bandit             <- BasicBernoulliBandit$new(weights = c(0.5, 0.2, 0.1)) #Distribuição Real

#Definição do Agente
agent              <- Agent$new(policy,bandit, "TS")

# Rodar a simulação
simulator          <- Simulator$new(agents      = agent,
                                    horizon     = 100, # Quantas vezes será realizado o experimento
                                    set_seed    = 8, # Semente para sempre vir os mesmos resultados
                                    save_theta  = TRUE,
                                    simulations = 100 # Quantas simulações realizadas
)

history <- simulator$run()


# Acessando Informações da Simulação:
summary(history)

## 
## Agents:
## 
##   TS
## 
## Cumulative regret:
## 
##  agent   t sims cum_regret cum_regret_var cum_regret_sd
##     TS 100  100       7.17       9.637475      3.104428
## 
## 
## Cumulative reward:
## 
##  agent   t sims cum_reward cum_reward_var cum_reward_sd
##     TS 100  100      42.21       44.14737      6.644349
## 
## 
## Cumulative reward rate:
## 
##  agent   t sims cur_reward cur_reward_var cur_reward_sd
##     TS 100  100     0.4221      0.4414737    0.06644349

# Gráfico de Escolha dos Braços
plot(history, type = "arms")

#Gráfico de avaliação da probabilidade média da recompensa
plot(history, type = "average", regret = FALSE)

# Grafico de avaliação da probabilidade acumulativa
plot(history, type = "cumulative", regret = FALSE, rate = TRUE)

Comparando resultados

bandit <- BasicBernoulliBandit$new(weights = c(0.5, 0.2, 0.1)) #Distribuição Real


EXpolicy <- EpsilonFirstPolicy$new(epsilon = 0.25, N = 100)
EXagent  <- Agent$new(EXpolicy,bandit, "Exploration")


EGpolicy <- EpsilonGreedyPolicy$new(epsilon = 0.1) 
EGagent  <- Agent$new(EGpolicy,bandit, "E-Greedy")

UCBpolicy <- UCB1Policy$new()
UCBagent  <- Agent$new(UCBpolicy,bandit, "UCB1")

TSpolicy <- ThompsonSamplingPolicy$new()
TSagent  <- Agent$new(TSpolicy,bandit, "TS")

agents <- list(EXagent, EGagent,UCBagent,TSagent)

simulator <- Simulator$new(agents,
                           horizon=100,
                           simulations = 10,
                           set_seed    = 8,
                           do_parallel = TRUE)

history_all <- simulator$run() 

#Gráfico de avaliação da probabilidade média da recompensa
plot(history_all, type = "average", regret = FALSE)

# Grafico de avaliação da probabilidade acumulativa
plot(history_all, type = "cumulative", regret = FALSE, rate = TRUE)

Através dos resultados apresentados temos que o método de \(\epsilon\)-Greedy foi o que obteve melhores resultados. Mas para uma estimativa melhor, devem ser realizadas mais simulações.

Exemplos Bases Reais - Movielens

Iremos aplicar a metodologia na base de dados Movielens com o objetivo de recomendar filmes para os usuários presentes na base. Essa base de dados contém 10 milhões de avaliações realizadas. Para acessar a base, iremos realizar download e logo em seguida iremos estruturá-la para a aplicação dos algoritmos.

Primeiro temos uma base de dados contendo informações dos filmes presentes no sistema. Essa base contém informações como: identificador do filme, nome do filme e gênero.

movies_dat      <- "http://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-10M/movies.dat"


movies_dat      <- readLines(movies_dat)
movies_dat      <- gsub( "::", "~", movies_dat )
movies_dat      <- paste0(movies_dat, collapse = "\n")
movies_dat      <- fread(movies_dat, sep = "~", quote="")
setnames(movies_dat, c("V1", "V2", "V3"), c("MovieID", "Name", "Type"))
movies_dat      <- splitstackshape::cSplit_e(movies_dat, "Type", sep = "|", type = "character",
                                             fill = 0, drop = TRUE)
movies_dat[[3]] <- NULL

kable(head(movies_dat,10))

Depois temos outra base contendo as avaliações realizadas pelos usuários cadastrados. Essa base contém informações como: identificador do usuário, identificador do filme avaliado pelo usuário, a avaliação e o timestamp da avaliação (dia e hora em que a avaliação foi realizada).

ratings_dat     <- "http://d1ie9wlkzugsxr.cloudfront.net/data_movielens/ml-10M/ratings.dat"
ratings_dat     <- RCurl::getURL(ratings_dat)
ratings_dat     <- readLines(textConnection(ratings_dat))
ratings_dat     <- gsub( "::", "~", ratings_dat )
ratings_dat     <- paste0(ratings_dat, collapse = "\n")
ratings_dat     <- fread(ratings_dat, sep = "~", quote="")
setnames(ratings_dat, c("V1", "V2", "V3", "V4"), c("UserID", "MovieID", "Rating", "Timestamp"))

all_movies      <- ratings_dat[movies_dat, on=c(MovieID = "MovieID")]

all_movies      <- na.omit(all_movies,cols=c("MovieID", "UserID"))

rm(movies_dat,ratings_dat)

kable(head(all_movies,10))

Para realizar esse experimento, iremos realizar a recomendação dos 50 filmes mais consumidos pelos usuários. Para isso, iremos construir uma base contendo todas as avaliações dos 50 filmes mais consumidos. Essa base também irá conter as informações dos filmes, como por exemplo, o gênero. Esses 50 filmes serão nossos “braços”.

# Data wrangling --------------------------------------------------------------------------------

all_movies[, UserID   := as.numeric(as.factor(UserID))]

count_movies    <- all_movies[,.(MovieCount = .N), by = MovieID]
top_50          <- as.vector(count_movies[order(-MovieCount)][1:50]$MovieID)
not_50          <- as.vector(count_movies[order(-MovieCount)][51:nrow(count_movies)]$MovieID)

top_50_movies   <- all_movies[MovieID %in% top_50]

kable(head(top_50_movies,10))

Para utilizar as informações de gênero dos filmes na recomendação, separamos os dados numa determinada base. Ou seja, para cada MovieID temos suas informações.

# Create feature lookup tables - to speed up, MovieID and UserID are
# ordered and lined up with the (dt/matrix) default index.

# Arm features

# MovieID of top 50 ordered from 1 to N as arms (Arms always start from 1)
top_50_movies[, MovieID   := as.numeric(as.factor(MovieID))]
arm_features    <- top_50_movies[,head(.SD, 1),by = MovieID][,c(1,6:24)]
setorder(arm_features,MovieID)

kable(head(arm_features,10))

Porém podemos extrair algumas informações sobre cada usuário único no sistema e utilizar isso a favor da recomendação. Utilizando as informações das avaliações dos filmes não selecionados (ou seja, sem ser os top 50), é possível avaliar qual gênero é mais consumido por cada usuário e utilziar essa informação para auxiliar na recomendação.

# User features

# Count of categories for non-top-50 movies normalized per user
user_features   <- all_movies[MovieID %in% not_50]
user_features[, c("MovieID", "Rating", "Timestamp", "Name"):=NULL]
user_features   <- user_features[, lapply(.SD, sum, na.rm=TRUE), by=UserID ]
user_features[, total  := rowSums(.SD, na.rm = TRUE), .SDcols = 2:20]
user_features[, 2:20 := lapply(.SD, function(x) x/total), .SDcols = 2:20]
user_features$total <- NULL

Com o objetivo de manter coerência com os dados, são acrescentados os usuários que não estavam sendo contabilizados acima. Ou seja, que dos filmes avaliados, nenhum foi visto por eles, o que estava ocorrendo com 4 usuários.

# Add users that were not in the set of non-top-50 movies 
all_users <- as.data.table(unique(all_movies$UserID))
user_features <- user_features[all_users, on=c(UserID = "V1")]
user_features[is.na(user_features)] <- 0

setorder(user_features,UserID)
rm(all_movies, not_50, top_50, count_movies)

kable(head(user_features,10))

Agora, visto que calculamos as informações e definimos nosso grupo, iremos montar e aplicar algoritmos de MAB.

# Contextual format

top_50_movies[, t := .I]
top_50_movies[, sim := 1]
top_50_movies[, agent := "Offline"]
top_50_movies[, choice := MovieID]
top_50_movies[, reward := ifelse(Rating <= 4, 0, 1)]

setorder(top_50_movies,Timestamp,Name)

kable(head(top_50_movies,10))

# Run simulation --------------------------------------------------------------------------------

simulations <- 10
horizon     <- 1000

bandit      <- OfflineLookupReplayEvaluatorBandit$new(top_50_movies,
                                                      k             = 50,
                                                      unique_col    = "UserID",
                                                      shared_lookup = arm_features,
                                                      unique_lookup = user_features)
agents      <-
  list(Agent$new(EpsilonFirstPolicy$new(epsilon = 0.25), bandit, "Uniform Exploration"),
       Agent$new(EpsilonGreedyPolicy$new(epsilon = 0.1) , bandit, "E-Greedy"),
       Agent$new(ThompsonSamplingPolicy$new(), bandit, "Thompson"),
       Agent$new(UCB1Policy$new(), bandit, "UCB1"),
       Agent$new(RandomPolicy$new(), bandit, "Random"),
       Agent$new(LinUCBHybridOptimizedPolicy$new(0.9), bandit, "LinUCB Hyb 0.9"),
       Agent$new(LinUCBDisjointOptimizedPolicy$new(2.1), bandit, "LinUCB Dis 2.1"))

simulation  <-
  Simulator$new(
    agents           = agents,
    simulations      = simulations,
    horizon          = horizon,
    set_seed    = 8
    
  )

results  <- simulation$run()

plot(results, type = "cumulative", regret = FALSE,
     rate = TRUE, legend_position = "topleft")

LinUCB Disjunt

Um algoritmo utilizado no exemplo anterior foi o de LinUCB. Iremos apresentar um pouco da sua metodologia.

De forma geral o LinUCB é um algoritmo de Bandits determinista e contextual que utiliza a regressão de Ridge para estimar os payoffs dos braços com base em contextos. Como demonstrado em [1], a utilização de uma solução de bandido contextual multi-armas em que a recompensa esperada \(r_{t,a}\) de um braço \(a\) é linear no seu contexto covariáveis \(x_{t,a}\) no momento t. Para essa situação, podemos chamar de Disjoint Linear UCB (LinUCB).

Na fórmula matricial, podemos representar o modelo de recompensa da seguinte forma:

\(E\left [ r_{t,a}|x_{t,a} \right ] = x_{t,a}^{T}\theta_{a}^{*}\)

Simplificando, pense no uso de um modelo de regressão linear em que a variável de resposta para cada braço é um retorno de recompensa com base nas covariáveis de entrada. Os coeficientes \(\theta_{a}^{*}\) são desconhecidos, mas com o algoritmo de aprendizado online, realizamos Regressão Ridge para obtê-lo em cada passo de tempo. De forma geral, a Regressão Ridge nos permite visualizar uma distribuição a posteriore dos coeficientes \(\theta_{a}^{*}\) e assim calcular \(p(\theta_{a}^{*})\) como uma estimativa pontual Bayesiana com uma média gaussiana \(\hat{\theta_{a}}\) e covariância \(A^{-1}_{a}\).

Dessa vez, dado que temos o valor esperado e um desvio padrão dado como \(\sqrt{x_{t,a}^{T} A^{-1}_{a}x_{t,a}}\)

No final do algoritmo, desejamos encontrar o braço com maior UCB no tempo t:

Ou seja, escolhe-se o braço com base na soma de:

• A estimativa média de recompensa de cada braço (caixa vermelha) e
• O limite de confiança superior correspondente (caixa azul), onde \(\alpha\) é um hiperparâmetro.

Quanto mais alto \(\alpha\), mais amplos se tornam os limites de confiança. Assim, resulta em uma maior ênfase colocada na exploração em vez da explotação.

LinUCB Hybrid

Outro algoritmo utilizado no exemplo dos filmes é o LinUCB Hybrid. Ele é considerado uma outra variante de um CB(Contextual Bandits) linear. Em certas aplicações, pode-se considerar que os braços não são mutuamente exclusivos em suas propriedades. Por exemplo, em um sistema de recomendação de notícias onde os artigos são braços recomendados aos usuários (com seu próprio contexto) para obter engajamento,temos os fatos de que alguns artigos podem ser semelhantes entre si e, portanto, possuem recursos compartilhados. Isso nos leva ao conceito de algoritmo LinUCB híbrido com modelos lineares híbridos (LinUCB Hybrid). O modelo híbrido afirma que a recompensa de cada braço é uma função linear de componentes compartilhados e não compartilhados , conforme mostrado por:

\(E\left [ r_{t,a}|x_{t,a} \right ] = z_{t,a}^{T}\beta_{*} + x_{t,a}^{T}\theta_{a}^{*}\)

onde,

• \(x_{t,a}\) representa as informações contextuais dos novos dados de observação (conforme apresentado no algoritmo LinUCB Disjoint)

• \(z_{t,a}\) representa a combinação de recursos do braço e atributos contextuais.

O modelo híbrido, como apresentado, tem um termo de recurso adicional \(z_{t}\). Dado que temos a noção de que cada braço tem seus próprios atributos, podemos criar um perfil para todos os braços sob as dimensões dos atributos. Com base nesses atributos, é possível identificar a semelhança ou dissimilaridade entre os braços. Usando as dimensões das informações do contexto de observação \(x_{t}\) e os atributos específicos para cada braço, podemos criar \(z_{t,a}\) um produto externo dessas duas informações.

Na equação do modelo linear, assim como \(x_{t,a}\) tem seus coeficientes correspondentes \(\theta^{*}_{a}\), \(z_{t,a}\) também tem seus coeficientes correspondentes \(\beta^{*}\). A principal diferença é que, embora \(\theta^{*}_{a}\) seja específico do braço , \(\beta^{*}\) é para todos os braços e está embutido na política Híbrida como um conjunto global de pesos. Essa intuição é capturada nas etapas do algoritmo, onde as informações de recompensa são usadas para atualizar para um braço escolhido. Assim, o braço escolhido atualiza seus coeficientes \(\theta^{*}_{a}\) com \(x_{t,a}\), mas, ao mesmo tempo, a política \(\beta^{*}\) também é atualizada.

Simplificando, no algoritmo Disjoint, atualizamos apenas as informações sobre as armas escolhidas em termos de recompensa como uma função linear das informações contextuais. No algoritmo Híbrido, atualizamos as informações não apenas para os braços escolhidos, mas também atualizamos como a recompensa muda com recursos do braço em comum.