Commit 1e2e9098 authored by Simon Ciranka's avatar Simon Ciranka
Browse files

added latest probabilistic sims and new figures

parent 57c72c28
Adaptive\_Adolescence\_multi
================
Simy
28/08/2020
# I do everything exactly the same as before
But I make multiple simulations because there is so much stochastizity
in the simluations and we need to get a feeling for the expectation
value.
``` r
legend=cowplot::get_legend(Kid)
cowplot::plot_grid(Kid+theme(legend.position="none"),
Adolescent+theme(legend.position="none",
axis.title.y = element_text(colour="white")),
legend,ncol =3,rel_widths = c(1,1,0.2))
```
![](Adaptive_Adolescence_Multi_files/figure-gfm/unnamed-chunk-2-1.png)<!-- -->
``` r
ggsave(filename = "../X_Figures/Environments.png")
```
## Saving 12 x 6 in image
# load the kalman filter and the UCB rule
Here i define the kalman filter and the choice rule. The Kalman Agent
explores the bandits by representing their mean reward rate and its
uncertainty about that mean. Both, beliefs about the mean reward rate
and uncertainty about that belief are updated every trials after
observing a new outcome. The Agents behavior can be gouverend by 3
parameters:
- mu0, or Optimism: its prior assumption about how rewarding the
envirionment will be on average. This Paramter will also gouvern the
extend of exploration
- var0, or speed of learning: its prior uncertainty about that mean.
the higher prior uncertainty is, the more extremely will the agent
update its beliefs in the beginning after observing a new sample
- vare, or imprecision: Imprecision will impact the speed of updating
as well. The higher the imprecision parameter, the more uncertain
the agent will remain about one options mean, even after observing a
lot of samples.
New samples are generated using an upper (lower) confidence bound rule
where the uncertainty about an option is added to (subtracted from) the
beleif about the mean of this option. This policy regulates exploration
as well where we can vary uncertaintyphile or phobe behavioral ploicies.
This proces is subject to paramter
- beta, or uncertainty aversion
(seeking)
<!-- end list -->
``` r
bayesianMeanTracker <- function(x, y, theta, prevPost=NULL,mu0Par,var0Par){
#Updates the previous posterior based on a single observation
#parameters
mu0 <- mu0Par #prior mean
var0 <- var0Par #prior variance
vare <- theta[1] #error varriance
if (is.null(prevPost)){#if no posterior prior, assume it is the first observation
predictions <- data.frame(mu=rep(mu0,144), sig=rep(var0,144))
}else{#if previous posterior is provided, update
predictions <- prevPost
}
#Which of the 121 options were chosen at time?
allopts<-expand.grid(1:12, 1:12)
chosen <- which(allopts$Var1==x[1] & allopts$Var2==x[2])
#Kalman gain
kGain <- predictions$sig[chosen] / (predictions$sig[chosen] + 3600)# feed the uncertainty in here.
#update mean
predictions$mu[chosen] <- predictions$mu[chosen] + (kGain * (y-predictions$mu[chosen]))
#update variance for observed arm
predictions$sig[chosen] <- predictions$sig[chosen] * (1 - kGain)
#return output
return(predictions)
}
class(bayesianMeanTracker)<- c(class(bayesianMeanTracker), "KalmanFilter")
ucb<-function(out, pars, refactor=F){
if (refactor==TRUE){
gamma <- pars[1]
beta_star<-pars[2]
#calulate all the upper confidence bounds
outtotal<-(gamma*out$mu)#+(beta_star*sqrt(out$sig)) #refactored parameters in combination with softmax tau, where gamma = 1/tau and beta_star = beta/tau
#avoid borderline cases
#outtotal[outtotal<=0]<-0.0001
#outtotal[outtotal>100]<-100
outtotal<-matrix(outtotal, ncol=1, byrow=TRUE)
}else{
beta <- pars[1]
#calulate all the upper confidence bounds
outtotal<-out$mu+(beta*sqrt(out$sig)) #refactored parameters in combination with softmax tau, where gamma = 1/tau and beta_star = beta/tau
#avoid borderline cases
#outtotal[outtotal<=0]<-0.0001
#outtotal[outtotal>99]<-99
outtotal<-matrix(outtotal, ncol=1, byrow=TRUE)
}
#return them
return(outtotal)
}
```
### some values to set.
``` r
# get sd of whole environemt for normalizing model input
set.seed(as.numeric(Sys.time()))
cntrl=list(
#get lambda
lambda=0.8,
#get beta
beta=0,# this scales risk attitude.
#get tau
tau=0.8,
mu0=100,#exploration bonus
var0=40,
#create a parameter vector
parVec = c(0.8, 0.8, 1, .0001) ,
#
ExploreBonus=0,
#kernel is RBF
#k<-rbf
#loop through trials
out=NULL,
AllChoices=NULL,
dummy=NULL,
overallCnt=1,
dat=expand.grid(x1=1:12,x2=1:12)
)
##
##
## Here i store the multiple Sims
ntrialss=100
list_Iter <- vector(mode = "list", length = ntrialss)
```
# Make observations
Here i let the agent learn about the envirionment. There are three
phases of the simulation. Phase one is “childhood”. During childhood
Agents can explore only the lower middle quadrant. Here all decisions
have low risk and there are some low rewards possible. Phase two, that
occurs after some learning experience (400 samples), can be understood
as the onset of adolescence. Here the whole environment becomes
availible but the agent does not know about it so they have to learn.
Then after another 400 samples, the agent transitions into “adulthood”
where the same environment is still present but the learning expierience
now lead to greater exploitation of presumably advantengeous options.
#### setup parallel
``` r
if(loadfromdisk==F){
library(doParallel)
#####setup parralel
cl <- makeCluster(40)
registerDoParallel(cl)
}
```
``` r
exploreEnv<-function(explore_func,choiceRule,env2,env1,cntrl,iter){
#for (rep in 1:ntrialss){
#unpack
lambda=cntrl$lambda
#get beta
beta<-cntrl$beta# this scales risk attitude.
#get tau
tau<-cntrl$tau
mu0<-cntrl$mu0#exploration bonus
var0<-cntrl$var0
#create a parameter vector
parVec <- cntrl$parVec
#
ExploreBonus=cntrl$ExploreBonus
#kernel is RBF
#k<-rbf
#loop through trials
out=cntrl$out
AllChoices=cntrl$AllChoices
dummy=cntrl$dummy
overallCnt=cntrl$overallCnt
dat=cntrl$dat
for (nround in 1:3){
#get parameters for participant on that round
if (nround==1){
# define vectors that are used by the kalman filter
lowestx=4
highestx=9
sampleVec=as.numeric(rownames(dat[dat$x1>=lowestx & dat$x1<=highestx & dat$x2<7,]))# here you define where a child should sample from
ind<-sample(sampleVec,1)
nTrials=400
}else {
ind<-sample(1:144,1)
nTrials=400
}
#random initialization as observation t=0
#y matrix
if (nround==1 & overallCnt==1){
X<-as.matrix(dat[ind,1:2])# generate a new vector of Xs
y<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
}else if(overallCnt==1) {
print("Youre an adolescent now")
X<-as.matrix(dat[ind,1:2])# generate a new vector of Xs
y<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
}
#X-start, i.e. all possible observations
Xstar<-as.matrix(dat[,1:2])
for (trial in 1:nTrials){
#output by GP with particular parameter settings
#don't forget mean centering and standardization.... mean is already 0 :)
if (overallCnt>1){
out<-bayesianMeanTracker(X[overallCnt,1:2],y[overallCnt], theta=lambda, prevPost=out,mu0Par=mu0,var0Par = var0)
}else{
out<-bayesianMeanTracker(X[overallCnt,1:2],y[overallCnt],theta=lambda, prevPost=NULL,mu0Par=mu0,var0Par=var0)
}
#utility vector. transpose if you use greedyMean
#where is everybody?
#here i need a function that calls bayesianMeanTracker. n times and returns the values X for each n. Also, i need some kind of list, where i save the prior for each instance....
#
# browser()
utilityVec<-ucb(out,beta)
utilities <- utilityVec - max(utilityVec)
#softmaximization
p <- exp(utilities/tau)
#probabilities
p <- p/colSums(p)
#numerical overflow
p <- (pmax(p, 0.00001))
p <- (pmin(p, 0.99999))
#index is sampled proprotionally to softmaxed utitily vector
if (nround==1){# subset the probability vector so that it corresponds to the right tiles.
ind<-sample(sampleVec,1,prob=p[dat$x1>=lowestx & dat$x1<=highestx & dat$x2<7,])# sample from a childhood environemnt
#this monster just scales exploration boni
}else {
ind<-sample(1:144, 1, prob=p)# sample from an adolescent environemnt
# print(ind)
}
X<-rbind(X, as.matrix(dat[ind,1:2]))
#bind y-observations
y<-rbind(y, as.matrix(rnorm(n=1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance)))# change this into a sample.
#if(y[overallCnt]<0){
# y[overallCnt]=-1*y[overallCnt]^2# make losses more severe.
#}
dummy<-data.frame(trial=overallCnt, x=as.numeric(X[overallCnt,1]), y=as.numeric(X[overallCnt,2]),
z=as.numeric(y[overallCnt]),round=nround)
AllChoices<-rbind(AllChoices,dummy)
overallCnt=overallCnt+1
}
#dummy data frame
}
#}
#This Here is for Plotting
Plot_dat=expand.grid(x=1:12,y=1:12,trials=0:max(dummy$trial))
Plot_dat$sample=0
Plot_dat$out=0
for (i in 1:length(AllChoices$x)){
AllChoices$y[i]
AllChoices$x[i]
Plot_dat[Plot_dat$x==AllChoices$x[i] & Plot_dat$y==AllChoices$y[i] & Plot_dat$trials==AllChoices$trial[i],]$sample=1
Plot_dat[Plot_dat$trials==AllChoices$trial[i],]$out=AllChoices$z[i]
}
Plot_dat$iter=iter
return(Plot_dat)
}
if(loadfromdisk==F){
Plot_datAll<-foreach(iter=1:ntrialss, .combine='rbind') %dopar%{
exploreEnv(explore_func=bayesianMeanTracker,choiceRule=ucb,env2=EnvirionemntAdol,env1=EnvirionemntKids,cntrl=cntrl,iter=iter)
}
saveRDS(file="../A_GeneratedFiles/Plot_datAll",object = Plot_datAll)
}else{
Plot_datAll<-readRDS(file="../A_GeneratedFiles/Plot_datAll")
}
# do it.
```
# Now lets contrast the relative & the total number of risky decisions.
Here i compare different developmental phases with respect to how many
rewards and losses have been encounterd as a consequence of exploring
this envirionment. An exploring agent likely samples highly
disadventegous options which agents with more expierience may avoid.
Therefore I also show how severe the losses were on average in each
developmental phase. An interesting divergence occurs here: While
encountered losses decrease and rewards increase across development; the
severity of losses is highest for the adolescent agents who still need
to explore the
environment.
``` r
Plot_datAll%>%group_by(x,y,iter)%>%filter(trials<=400)%>%mutate(cumDens=cumsum(sample),
Outcome=case_when(
(out<0)~"Loss",
(out>=0)~"Gain"
)
)%>%group_by(Outcome,iter)%>%summarise(howMany=n()/144,
howMuch=sum(out))->Kids
```
## `summarise()` regrouping output by 'Outcome' (override with `.groups` argument)
``` r
Plot_datAll%>%group_by(x,y,iter)%>%filter(trials>=400 & trials<800)%>%mutate(cumDens=cumsum(sample),
Outcome=case_when(
(out<0)~"Loss",
(out>=0)~"Gain"
)
)%>%group_by(Outcome,iter)%>%summarise(howMany=n()/144,
howMuch=sum(out))->Adolescents
```
## `summarise()` regrouping output by 'Outcome' (override with `.groups` argument)
``` r
Plot_datAll%>%group_by(x,y,iter)%>%filter(trials>800 & trials<2400)%>%mutate(cumDens=cumsum(sample),
Outcome=case_when(
(out<0)~"Loss",
(out>=0)~"Gain"
)
)%>%group_by(Outcome,iter)%>%summarise(howMany=n()/144,
howMuch=sum(out)
)->YoungAdults
```
## `summarise()` regrouping output by 'Outcome' (override with `.groups` argument)
``` r
Adolescents$Stage="2"
Kids$Stage="1"
YoungAdults$Stage="3"
Adolescents$Soc="0"
Kids$Soc="0"
YoungAdults$Soc="0"
Plot_labels=c("Children","Adolescents","Adults")
rbind(Kids,Adolescents,YoungAdults)->AllIndi
```
``` r
AllIndi%>%ungroup()%>%select(howMany)%>%min()->minScale
AllIndi%>%ungroup()%>%select(howMany)%>%max()->maxScale
#first show distribution of risky decisions as kid
ggplot(Kids,aes(x=Outcome,y=howMany,fill=Outcome))+
geom_jitter(aes(group=iter,color=Outcome),alpha=0.3)+
stat_summary(aes(color=Outcome),geom="pointrange",color="black",fun.data = "mean_cl_boot")+
coord_cartesian(ylim=c(0,400))+
scale_y_continuous(name="#")+
ggtitle("outcomes children\nsolo")+
guides(fill=F)+
theme_minimal(8)->One
#then show distribution of risky decisions as adolescent
ggplot(Adolescents,aes(x=Outcome,y=howMany,fill=Outcome))+
geom_jitter(aes(group=iter,color=Outcome),alpha=0.3)+
stat_summary(aes(color=Outcome),geom="pointrange",color="black",fun.data = "mean_cl_boot")+
coord_cartesian(ylim=c(0,400))+
scale_y_continuous(name="")+
ggtitle("outcomes adolescents\nsolo")+
guides(fill=F)+
theme_minimal(8)->Two
#then Adults distribution of risky decisions as adolescent
ggplot(YoungAdults,aes(x=Outcome,y=howMany,fill=Outcome))+
geom_jitter(aes(group=iter,color=Outcome),alpha=0.3)+
stat_summary(aes(color=Outcome),geom="pointrange",color="black",fun.data = "mean_cl_boot")+
scale_y_continuous(name="")+
coord_cartesian(ylim=c(0,400))+
guides(fill=F)+
ggtitle("outcomes adults\nsolo")+
theme_minimal(8)->Three
#Now, look at the outcome
ggplot(AllIndi[AllIndi$Outcome=="Loss",],aes(x=Stage,y=howMuch,fill=Stage))+
geom_jitter(aes(group=iter,color=Stage),alpha=0.3)+
stat_summary(aes(color=Stage),geom="pointrange",color="black",fun.data = "mean_cl_boot")+
scale_y_continuous(name="Cumulative Loss")+
scale_x_discrete("Developmental Stage",breaks=c("1","2","3"),labels=Plot_labels)+
ggtitle("Severety of bad\noutcomes")+
scale_color_brewer(palette = "Set2")+
coord_cartesian(ylim=c(-2000000,0))+
guides(fill=F)+
theme_minimal(8)->Four
cowplot::plot_grid(One,Two,Three,NULL,Four,ncol=5,rel_widths = c(1,1,1,0.3,1.7))->Solo
ggsave(plot = Solo,"../X_Figures/EmergentBehavior_summary.png",width = 11,height = 4)
print(Solo)
```
![](Adaptive_Adolescence_Multi_files/figure-gfm/unnamed-chunk-7-1.png)<!-- -->
# Number of explorative decisions
Not only the outcomes are interesting but also how much exploraiton
happens. So i now declare an explorative decision as switching a bandit.
This is done by observing whether sampling of one xy pair changes from 0
to 1.
``` r
#this is why i love dplyr!
#create count for each new decision.
# do it in parralel again
compute_sampling_strategy<-function(Plot_datAll,subject){
library(dplyr)
Plot_datAll%>%filter(iter==subject)%>%group_by(x,y)%>%mutate(cumDens=cumsum(sample))%>%ungroup()%>% group_by(x,y, trials) %>%
mutate(newC = ifelse(cumDens == 1 & lag(cumDens)==0, 1, 0))->newExp# if i take the derivative of this this might be exaclty what wouter was intersted in.
newExp[is.na(newExp$newC),]$newC=1#
#now i have to "reset" so that each trial there can only be one new decision and because otherwise this will scale up the cumulative sum too much
newExp$newExplore=0
newExp$generalExplore=0
for( i in 1:length(unique(newExp$trials))){
newExp[newExp$trials==i,]$newExplore=xor(newExp[newExp$trials==i,]$newC,newExp[newExp$trials==i-1,]$newC)#
print(i)
# newExp[newExp$trials==i,]$generalExplore=xor(newExp[newExp$trials==i,]$sample,newExp[newExp$trials==i-1,]$sample)
}
newExp$iter=subject
return(newExp)
}
#saveRDS(newExp,"Derivatives_TrialWise.rds")
if (loadfromdisk==F){
#newExp_all<-NULL
#for(subject in 1:ntrialss){
# print(subject)
newExp_all<-foreach(subject = 1:ntrialss,.combine = "rbind",.verbose=T) %dopar%{
compute_sampling_strategy(Plot_datAll,subject)
}
saveRDS(file = "../A_GeneratedFiles/solotraj2.rds",object = newExp_all)
}else {
newExp_all<-readRDS(file = "../A_GeneratedFiles/solotraj2.rds")
}
# to compute the derivative
```
I then compute the cumulative sum of this count and subtract the
cumulative count from 50 trials earlier from it, which gives me a
“slope” of new explorative decisions. But i also contrast this to
the raw cumulative count, in order for us to inspect how these values
correspond.
``` r
newExp_all%>%group_by(iter)%>%
mutate(newCum=cumsum(newExplore))%>%
ungroup()%>%group_by(trials,iter)%>%
summarise(explore=max(newCum))%>%
ungroup()%>%
arrange(iter)%>%
#ungroup()%>%#### looks about right but doublecheck.
mutate(
explore2=explore-lag(explore,50)
)%>%mutate(
Stage=case_when(
(trials<401)~"Children",
(trials>400&trials<800)~"Adolescents",
(trials>800)~"Adults"
),
Which="Solo"
)%>%#filter(trials>36 & (trials <401 | trials >508))%>%
filter(!is.na(Stage) & trials %in% round(seq(50,1200,length.out = 24)))%>%
ggplot(aes(x=trials,y=explore2))+
stat_summary(geom="ribbon", fun.data="mean_sdl",alpha=0.2)+
stat_summary(geom="line",fun.y = "mean",color="black")+
# geom_jitter()+
#geom_rect(aes(xmin=50,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
#geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
#geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
geom_spline(size=1)+
geom_point(aes(shape=Stage),size=1,alpha=0.1)+
stat_summary(aes(shape=Stage),geom="point",fun.y = "mean",size=2.5,color="black",fill="white")+
scale_shape_manual(name="Developmental\nStage",values=c(22,21,24))+
geom_vline(xintercept = 420, linetype="dotted",color="red")+
geom_vline(xintercept = 410, linetype="dotted",color="red")+
geom_vline(xintercept = 800, linetype="dotted",color="red")+
geom_vline(xintercept = 790, linetype="dotted",color="red")+
annotate("text",x=150,y=60,label=c("Children"))+
annotate("text",x=600,y=60,label=c("Adolescents"))+
annotate("text",x=950,y=60,label=c("Adults"))+
#scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
#coord_cartesian(ylim=c(0,))+
scale_y_continuous(name="Exploration decisions")+
scale_x_continuous(name="trials")+
guides(linetype=F)+
theme_minimal(14)
```
## `summarise()` regrouping output by 'trials' (override with `.groups` argument)
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.
![](Adaptive_Adolescence_Multi_files/figure-gfm/unnamed-chunk-9-1.png)<!-- -->
``` r
#->ExploreSoloPlot
#ggsave(plot=ExploreSoloPlot,filename="../X_Figures/TimecourseExplore.png",width=6.3,height=5)
```
``` r
newExp_all%>%group_by(iter)%>%
mutate(newCum=cumsum(newExplore))%>%
ungroup()%>%group_by(trials,iter)%>%
summarise(explore=max(newCum))%>%
ungroup()%>%
arrange(iter)%>%
#ungroup()%>%#### looks about right but doublecheck.
mutate(
explore2=explore-lag(explore,50)
)%>%mutate(
Stage=case_when(
(trials<401)~"1",
(trials>400&trials<800)~"2",
(trials>800)~"3"
)
)%>%filter(!is.na(Stage) & trials %in% round(seq(50,1200,length.out = 24)))%>%group_by(Stage,iter)%>%
summarize(AUC=trapz(trials,explore2))->AUCSoloTbl
```
## `summarise()` regrouping output by 'trials' (override with `.groups` argument)
## `summarise()` regrouping output by 'Stage' (override with `.groups` argument)
``` r
ylimauc<-AUCSoloTbl%>%ungroup()%>%select(AUC)%>%max()
ggplot(AUCSoloTbl,aes(x=Stage,y=AUC))+
#geom_jitter(aes(group=iter,color=Stage),alpha=0.1)+
stat_summary(aes(x=Stage,y=AUC,fill=Stage),geom="bar",color="black",fun.y = "mean",stat="identity")+
stat_summary(aes(color=Stage),geom="pointrange",color="black",fun.data = "mean_cl_boot")+
scale_x_discrete(name="Developmental Stage",breaks=c(1,2,3),labels=c("Children","Adolescents","Adults"))+
ggtitle("Exploration (AUC) - Solo")+
scale_fill_brewer(palette = "Set2",breaks=c(1,2,3),labels=c("Children","Adolescents","Adults"))+
scale_y_continuous(name=c("AUC Explore"),limits=c(0,ylimauc))+
theme_minimal(14)->AUCSolo
```
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: Ignoring unknown parameters: stat
# Social
One particular feature about adolescents is that they are sensitive to
social information. In order to invesitage how such a sensitivity
impacts 1) exploration behavior and 2) encountered outcomes. For this i
let *n* agents perform the task in parallel and introduce a bonus to the
choice utilites that depends on how many others have visited this option
in the previous trial.
``` r
# get sd of whole environemt for normalizing model input
set.seed(as.numeric(Sys.time()))
loadfromdisk=T
###
###
###
###
### PARAMETERS
###
###
###
#get lambda
#get beta
cntrl_social<-list(
beta=