Commit 002ea0ed authored by Simon Ciranka's avatar Simon Ciranka

i implemented MLE and Wouters version of the Bayesian Updating

parent 4ae8bef6
......@@ -10,8 +10,6 @@ knitr::opts_chunk$set(echo = TRUE)
library(sigmoid)
library(tidyverse)
# this is the range of possible outcome probabilities.
posResp = seq(0,1,0.01);
posResp = round(posResp,2);
# this is Wouters Theme.
my_theme <- theme_bw( )+theme(axis.text = element_text(size = 18),
axis.title = element_text(size = 18),
......@@ -25,9 +23,21 @@ my_theme <- theme_bw( )+theme(axis.text = element_text(size = 18),
```
# Implementing Binomal Updating in Wouters way.
These code chunks and results are implementing the models in wouters mail form the 11.9.2018. Other than before, we now caluculate the log likelihood of our participants estimates by assuming that the range of possible estimates is normally distributed. The probability density under the specific probability estimate of the participant is now the likelihood. For this i need a function that returns discrete values for the probability denisty in question.
I implemented the function you sent to me in your mail. At some point when you talk about the arguments that you pass to the cumulative densitiy function which you use to get a log likelihood; you say sigma is sigma.
In the example you sent to me; you do not use sigma as the variance estimate of the beta updating before. You pass sigma with the function one higher level that you then minimize. Is sigma as you use it a free parameter of the model as well? So you dont have one but two free parameters?
discretenormalMarble005(posResp,model_weight, sigma); % THIS IS THE CUSTOM FUNCTION WE BUILT FOR THIS STEP SEE BELOW (posResp is list made above. I guess of all possible resonses a subject can make).. sigma is sigma..
# Make the cumulative Densitiy
```{r Discrete Function}
posResp = seq(0,1,0.01);
posResp = round(posResp,2);
#sigma=0.1
#here i define a function that gives us the probability densitiy of all possible values given the mean and the standart deviation
#of the posterior beta distribution that results from binomal updating.
discretenormalMarble005<-function (possibleResp, mean, stdev){
......@@ -61,6 +71,7 @@ To allow for individual and developmental differences in the updating process I
```{r linear,message=FALSE, warning=FALSE}
linearWeights<- function(v){
sigma=1;
lr<-v[1]
for (i in 1:nrow(subjectLevel)){
#reset shape paramters
......@@ -86,7 +97,7 @@ linearWeights<- function(v){
#calculate the sufficient statistics of the beta distribution.
mu<-aBeta/(aBeta+bBeta)# calculate only the last mu.
sigma<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
var<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
subEstimate=subjectLevel$p.response[i]/100
......@@ -96,8 +107,8 @@ linearWeights<- function(v){
# get pr and m. I need to dig into what this means exactly.
# i use the estiamted mean and sigma to be able to get descrete log likelihoods.
pr = discretenormalMarble005(posResp,mu, sigma)[[1]];
m = discretenormalMarble005(posResp,mu, sigma)[[2]];
pr = discretenormalMarble005(posResp,mu, 0.1)[[1]];
m = discretenormalMarble005(posResp,mu, 0.1)[[2]];
#check where the possible resposnes equal the exact response of the subject. now we have the response liklihood.
probEstimate = pr[which(posResp==subEstimate)]; # SEE HERE YOU FIND THE PROB EST THAT FITS THE GIVEN SUBJECTIVE VALUE
#zero doesnt work.
......@@ -127,6 +138,7 @@ $$ Posterior_{p(blue)} \sim Beta\bigg(\alpha+ k ^\delta , \beta + (N-k)^\delta\
This assumes a more holistic processing but i think it may actually be a good heuristic given that the stimuli are presented such a short amount of time.
```{r expo,message=FALSE, warning=FALSE}
exponentialWeights<- function(v){
sigma=1;
lr<-v[1]
for (i in 1:nrow(subjectLevel)){
#reset shape paramters
......@@ -152,8 +164,7 @@ exponentialWeights<- function(v){
#calculate the sufficient statistics of the beta distribution.
mu<-aBeta/(aBeta+bBeta)# calculate only the last mu.
sigma<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
var<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
subEstimate=subjectLevel$p.response[i]/100
# correct extreme cases.
......@@ -162,8 +173,8 @@ exponentialWeights<- function(v){
# get pr and m. I need to dig into what this means exactly.
# i use the estiamted mean and sigma to be able to get descrete log likelihoods.
pr = discretenormalMarble005(posResp,mu, sigma)[[1]];
m = discretenormalMarble005(posResp,mu, sigma)[[2]];
pr = discretenormalMarble005(posResp,mu, 0.1)[[1]];
m = discretenormalMarble005(posResp,mu, 0.1)[[2]];
#check where the possible resposnes equal the exact response of the subject. now we have the response liklihood.
probEstimate = pr[which(posResp==subEstimate)]; # SEE HERE YOU FIND THE PROB EST THAT FITS THE GIVEN SUBJECTIVE VALUE
#zero doesnt work.
......@@ -272,7 +283,7 @@ data %>% gather( key = ModelLik, value = GSquared, LLLinear, LLExp) %>%
xlab("Models")+
scale_color_discrete(
name="Models",
breaks = c("LLDisc", "LLSimple"),
breaks = c("LLLinear", "LLExp"),
labels = c("Exponential Weight Beta Update", " Weight Beta Update"))+
my_theme
```
......@@ -290,7 +301,7 @@ learningPlot<-ggplot(,aes(x=1,y=unique(data$learningRateLinear)))+geom_violin(co
annotate("text", x=0.65, y=1, vjust = -1, label = "Ideal Observer")+
my_theme
learningDiscounting<-ggplot(,aes(x=1,y=unique(data$learningRateExp)))+geom_violin(color="blue",fill="red",alpha=0.1)+geom_jitter(width = 0.2)+geom_boxplot(width=0.1, alpha=0.5)+
learningDiscounting<-ggplot(,aes(x=1,y=unique(data$learningRateExp)))+geom_violin(color="blue",fill="blue",alpha=0.1)+geom_jitter(width = 0.2)+geom_boxplot(width=0.1, alpha=0.5)+
ggtitle("LearningRate Distribution Marbles")+
ylab("Learning Rate Estimate: Sequential")+
xlab("Data Density")+geom_hline(yintercept = 1) +
......@@ -302,14 +313,17 @@ learningDiscounting<-ggplot(,aes(x=1,y=unique(data$learningRateExp)))+geom_violi
### Marble Estimate Distribution
This plot shows how the learning rates are distributed in all subjects. We can see that most of the subjects seem to overweight new pieces of information
relative to an Ideal observer.
The pattern which we had before remains. If i use the linear weighting, it looks like new information is overweighted.
```{r ShowPlot, results="markup",echo=FALSE}
print(learningPlot)
```
### Marble Estimate Distribution
What i get now if i look at the parameter estimate of the exponential model is an **underweghting** of all information considered.
```{r Show second Plot, results="markup",echo=FALSE}
print(learningDiscounting)
```
The HGF like model has no Benchmark and the "learning rate can be
\ No newline at end of file
\ No newline at end of file
This diff is collapsed.
......@@ -9,12 +9,31 @@ under the specific probability estimate of the participant is now the
likelihood. For this i need a function that returns discrete values for
the probability denisty in question.
I implemented the function you sent to me in your mail. At some point
when you talk about the arguments that you pass to the cumulative
densitiy function which you use to get a log likelihood; you say sigma
is sigma. In the example you sent to me; you do not use sigma as the
variance estimate of the beta updating before. You pass sigma with the
function one higher level that you then minimize. **Is sigma as you use it
a free parameter of the model as well?** So you dont have one but two free
parameters?
discretenormalMarble005(posResp,model_weight, sigma); % THIS IS THE CUSTOM FUNCTION WE BUILT FOR THIS STEP SEE BELOW (posResp is list made above. I guess of all possible resonses a subject can make).. sigma is sigma..
Make the cumulative Densitiy
============================
```r
posResp = seq(0,1,0.01);
posResp = round(posResp,2);
#sigma=0.1
#here i define a function that gives us the probability densitiy of all possible values given the mean and the standart deviation
#of the posterior beta distribution that results from binomal updating.
discretenormalMarble005<-function (possibleResp, mean, stdev){
lower = min(possibleResp);
upper = max(possibleResp);
top = seq((lower+0.005),(upper+0.005),0.01);
bottom = seq((lower-0.005),(upper-0.005),0.01);
......@@ -24,6 +43,7 @@ the probability denisty in question.
m = p*possibleResp;
return(list(p,m))
}
```
Sequential Updating
-------------------
......@@ -59,7 +79,10 @@ than 1 and overweighting larger than 1.
The mean of this posterior is then compared to the actual probability
estimate of the participants to create logliks.
```r
linearWeights<- function(v){
sigma=1;
lr<-v[1]
for (i in 1:nrow(subjectLevel)){
#reset shape paramters
......@@ -85,7 +108,7 @@ estimate of the participants to create logliks.
#calculate the sufficient statistics of the beta distribution.
mu<-aBeta/(aBeta+bBeta)# calculate only the last mu.
sigma<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
var<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
subEstimate=subjectLevel$p.response[i]/100
......@@ -95,8 +118,8 @@ estimate of the participants to create logliks.
# get pr and m. I need to dig into what this means exactly.
# i use the estiamted mean and sigma to be able to get descrete log likelihoods.
pr = discretenormalMarble005(posResp,mu, sigma)[[1]];
m = discretenormalMarble005(posResp,mu, sigma)[[2]];
pr = discretenormalMarble005(posResp,mu, 0.1)[[1]];
m = discretenormalMarble005(posResp,mu, 0.1)[[2]];
#check where the possible resposnes equal the exact response of the subject. now we have the response liklihood.
probEstimate = pr[which(posResp==subEstimate)]; # SEE HERE YOU FIND THE PROB EST THAT FITS THE GIVEN SUBJECTIVE VALUE
#zero doesnt work.
......@@ -115,6 +138,7 @@ estimate of the participants to create logliks.
G2 <- -2*G2
G2
}
```
Exponential Discount Factor.
----------------------------
......@@ -127,7 +151,10 @@ it to the power of some discount factor *δ*
a good heuristic given that the stimuli are presented such a short
amount of time.
```r
exponentialWeights<- function(v){
sigma=1;
lr<-v[1]
for (i in 1:nrow(subjectLevel)){
#reset shape paramters
......@@ -153,8 +180,7 @@ amount of time.
#calculate the sufficient statistics of the beta distribution.
mu<-aBeta/(aBeta+bBeta)# calculate only the last mu.
sigma<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
var<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
subEstimate=subjectLevel$p.response[i]/100
# correct extreme cases.
......@@ -163,8 +189,8 @@ amount of time.
# get pr and m. I need to dig into what this means exactly.
# i use the estiamted mean and sigma to be able to get descrete log likelihoods.
pr = discretenormalMarble005(posResp,mu, sigma)[[1]];
m = discretenormalMarble005(posResp,mu, sigma)[[2]];
pr = discretenormalMarble005(posResp,mu, 0.1)[[1]];
m = discretenormalMarble005(posResp,mu, 0.1)[[2]];
#check where the possible resposnes equal the exact response of the subject. now we have the response liklihood.
probEstimate = pr[which(posResp==subEstimate)]; # SEE HERE YOU FIND THE PROB EST THAT FITS THE GIVEN SUBJECTIVE VALUE
#zero doesnt work.
......@@ -183,6 +209,7 @@ amount of time.
G2 <- -2*G2
G2
}
```
Data Loading
------------
......@@ -191,16 +218,20 @@ In this Chunk of Code i load the Data which i made with first loading
the rawdata in matblab, and squeezing the struct into two dimensional
data and run the script [01\_makeDataFrame.R](01_makeDataFrame.R)
```r
load("RobertsMarbleDf.RData")
data$sub.id<-as.numeric(data$sub.id)
Subs<-unique(data$sub.id)
data$sequence.marbles.color1<-as.character(data$sequence.marbles.color1) #blue
data$sequence.marbles.color2<-as.character(data$sequence.marbles.color2) #red
sub.list<-list()
```
Here i Fit the Simple LearningRate Model.
-----------------------------------------
```r
for (i in 1:length(Subs)){
subjectLevel<-data[data$sub.id==Subs[i],]
output<-optim(c(1), fn = linearWeights, method = c("Brent"),upper = 5,lower = 0)
......@@ -226,10 +257,12 @@ Here i Fit the Simple LearningRate Model.
data$learningRateLinear[i] = toMerge[toMerge$PPN == data$sub.id[i], ]$lr
data$LLLinear[i] = toMerge[toMerge$PPN == data$sub.id[i], ]$LL_win
}
```
Here i Fit the Discount LearningRate Model.
-------------------------------------------
``` r
for (i in 1:length(Subs)){
subjectLevel<-data[data$sub.id==Subs[i],]
output<-optim(c(1), fn = exponentialWeights, method = c("Brent"),upper = 10,lower = 0)
......@@ -255,6 +288,7 @@ Here i Fit the Discount LearningRate Model.
data$learningRateExp[i] = toMerge[toMerge$PPN == data$sub.id[i], ]$delta
data$LLExp[i] = toMerge[toMerge$PPN == data$sub.id[i], ]$LL_win
}
```
Model Comparison
----------------
......@@ -264,6 +298,7 @@ Model”, the “sequential Updating” and the “exponential discounting”
model. The discounting model and but the HGF Like Model are quite close.
Seqential Updating is bad.
```r
data %>% gather( key = ModelLik, value = GSquared, LLLinear, LLExp) %>%
distinct(GSquared,ModelLik) %>%
ggplot(aes(x=as.factor(ModelLik),y=GSquared,color=as.factor(ModelLik)))+
......@@ -272,9 +307,10 @@ Seqential Updating is bad.
xlab("Models")+
scale_color_discrete(
name="Models",
breaks = c("LLDisc", "LLSimple"),
breaks = c("LLLinear", "LLExp"),
labels = c("Exponential Weight Beta Update", " Weight Beta Update"))+
my_theme
```
![](LiklelihoodUpdate_files/figure-markdown_strict/unnamed-chunk-1-1.png)
......@@ -285,12 +321,15 @@ So now lets look at the learning rates.
This plot shows how the learning rates are distributed in all subjects.
We can see that most of the subjects seem to overweight new pieces of
information relative to an Ideal observer.
information relative to an Ideal observer. The pattern which we had
before remains. If i use the linear weighting, it looks like new
information is overweighted.
![](LiklelihoodUpdate_files/figure-markdown_strict/ShowPlot-1.png)
\#\#\# Marble Estimate Distribution What i get now if i look at the
parameter estimate of the exponential model is an **underweghting** of
all information considered.
![](LiklelihoodUpdate_files/figure-markdown_strict/Show%20second%20Plot-1.png)
The HGF like model has no Benchmark and the “learning rate can be
### Marble Estimate Distribution
What i get now if i look at the parameter estimate of the exponential
model is an **underweghting** of all information considered.
![](LiklelihoodUpdate_files/figure-markdown_strict/Show%20second%20Plot-1.png)
......@@ -9,14 +9,31 @@ under the specific probability estimate of the participant is now the
likelihood. For this i need a function that returns discrete values for
the probability denisty in question.
```r
I implemented the function you sent to me in your mail. At some point
when you talk about the arguments that you pass to the cumulative
densitiy function which you use to get a log likelihood; you say sigma
is sigma. In the example you sent to me; you do not use sigma as the
variance estimate of the beta updating before. You pass sigma with the
function one higher level that you then minimize. **Is sigma as you use it
a free parameter of the model as well?** So you dont have one but two free
parameters?
discretenormalMarble005(posResp,model_weight, sigma); % THIS IS THE CUSTOM FUNCTION WE BUILT FOR THIS STEP SEE BELOW (posResp is list made above. I guess of all possible resonses a subject can make).. sigma is sigma..
Make the cumulative Densitiy
============================
```r
posResp = seq(0,1,0.01);
posResp = round(posResp,2);
#sigma=0.1
#here i define a function that gives us the probability densitiy of all possible values given the mean and the standart deviation
#of the posterior beta distribution that results from binomal updating.
discretenormalMarble005<-function (possibleResp, mean, stdev){
lower = min(possibleResp);
upper = max(possibleResp);
top = seq((lower+0.005),(upper+0.005),0.01);
bottom = seq((lower-0.005),(upper-0.005),0.01);
......@@ -63,7 +80,9 @@ than 1 and overweighting larger than 1.
estimate of the participants to create logliks.
```r
linearWeights<- function(v){
sigma=1;
lr<-v[1]
for (i in 1:nrow(subjectLevel)){
#reset shape paramters
......@@ -89,7 +108,7 @@ estimate of the participants to create logliks.
#calculate the sufficient statistics of the beta distribution.
mu<-aBeta/(aBeta+bBeta)# calculate only the last mu.
sigma<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
var<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
subEstimate=subjectLevel$p.response[i]/100
......@@ -99,8 +118,8 @@ estimate of the participants to create logliks.
# get pr and m. I need to dig into what this means exactly.
# i use the estiamted mean and sigma to be able to get descrete log likelihoods.
pr = discretenormalMarble005(posResp,mu, sigma)[[1]];
m = discretenormalMarble005(posResp,mu, sigma)[[2]];
pr = discretenormalMarble005(posResp,mu, 0.1)[[1]];
m = discretenormalMarble005(posResp,mu, 0.1)[[2]];
#check where the possible resposnes equal the exact response of the subject. now we have the response liklihood.
probEstimate = pr[which(posResp==subEstimate)]; # SEE HERE YOU FIND THE PROB EST THAT FITS THE GIVEN SUBJECTIVE VALUE
#zero doesnt work.
......@@ -133,7 +152,9 @@ a good heuristic given that the stimuli are presented such a short
amount of time.
```r
exponentialWeights<- function(v){
sigma=1;
lr<-v[1]
for (i in 1:nrow(subjectLevel)){
#reset shape paramters
......@@ -159,8 +180,7 @@ amount of time.
#calculate the sufficient statistics of the beta distribution.
mu<-aBeta/(aBeta+bBeta)# calculate only the last mu.
sigma<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
var<-(aBeta*bBeta)/((aBeta+bBeta)^2*(aBeta+bBeta+1))#calculate the sigma.
subEstimate=subjectLevel$p.response[i]/100
# correct extreme cases.
......@@ -169,8 +189,8 @@ amount of time.
# get pr and m. I need to dig into what this means exactly.
# i use the estiamted mean and sigma to be able to get descrete log likelihoods.
pr = discretenormalMarble005(posResp,mu, sigma)[[1]];
m = discretenormalMarble005(posResp,mu, sigma)[[2]];
pr = discretenormalMarble005(posResp,mu, 0.1)[[1]];
m = discretenormalMarble005(posResp,mu, 0.1)[[2]];
#check where the possible resposnes equal the exact response of the subject. now we have the response liklihood.
probEstimate = pr[which(posResp==subEstimate)]; # SEE HERE YOU FIND THE PROB EST THAT FITS THE GIVEN SUBJECTIVE VALUE
#zero doesnt work.
......@@ -207,9 +227,10 @@ data and run the script [01\_makeDataFrame.R](01_makeDataFrame.R)
sub.list<-list()
```
# Here i Fit Both Models
Here i Fit the Simple LearningRate Model.
-----------------------------------------
``` r
```r
for (i in 1:length(Subs)){
subjectLevel<-data[data$sub.id==Subs[i],]
......@@ -236,10 +257,12 @@ data and run the script [01\_makeDataFrame.R](01_makeDataFrame.R)
data$learningRateLinear[i] = toMerge[toMerge$PPN == data$sub.id[i], ]$lr
data$LLLinear[i] = toMerge[toMerge$PPN == data$sub.id[i], ]$LL_win
}
########## this fits the exponential model
```
Here i Fit the Discount LearningRate Model.
-------------------------------------------
``` r
for (i in 1:length(Subs)){
subjectLevel<-data[data$sub.id==Subs[i],]
output<-optim(c(1), fn = exponentialWeights, method = c("Brent"),upper = 10,lower = 0)
......@@ -267,10 +290,13 @@ data and run the script [01\_makeDataFrame.R](01_makeDataFrame.R)
}
```
# Model Comparison
Model Comparison
----------------
Here i judge via G^2 which model is the best. Lower values is better fit.
Here i judge via G^2 which model is the best. I compare the “HGF Like
Model”, the “sequential Updating” and the “exponential discounting”
model. The discounting model and but the HGF Like Model are quite close.
Seqential Updating is bad.
```r
data %>% gather( key = ModelLik, value = GSquared, LLLinear, LLExp) %>%
......@@ -281,7 +307,7 @@ Here i judge via G^2 which model is the best. Lower values is better fit.
xlab("Models")+
scale_color_discrete(
name="Models",
breaks = c("LLDisc", "LLSimple"),
breaks = c("LLLinear", "LLExp"),
labels = c("Exponential Weight Beta Update", " Weight Beta Update"))+
my_theme
```
......@@ -295,11 +321,15 @@ So now lets look at the learning rates.
This plot shows how the learning rates are distributed in all subjects.
We can see that most of the subjects seem to overweight new pieces of
information relative to an Ideal observer.
information relative to an Ideal observer. The pattern which we had
before remains. If i use the linear weighting, it looks like new
information is overweighted.
## This corresponds to the linear weighting model
![](LiklelihoodUpdate_files/figure-markdown_strict/ShowPlot-1.png)
## This corresponds to the exponential weighting model.
![](LiklelihoodUpdate_files/figure-markdown_strict/Show%20second%20Plot-1.png)
### Marble Estimate Distribution
What i get now if i look at the parameter estimate of the exponential
model is an **underweghting** of all information considered.
![](LiklelihoodUpdate_files/figure-markdown_strict/Show%20second%20Plot-1.png)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment