---
title: "Computational Cognitive Science 2022-2023"
output:
  pdf_document: default
  word_document: default
  html_document: default
---
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(DirichletReg) # We'll be using Dirichlet distributions. gtools includes rdirichlet and ddirichlet but not with log densities, which are useful.
```

# Tutorial 4: Bayesian Model Comparison

## Nisbett and Wilson: Are people biased?

We previously used count data from a Nisbett and Wilson study to look at different approaches to parameter estimation. We will now use the same data to look at model comparison.

Options A-D are identical consumer products arranged from left to right in a display.
Observations are the number of participants who chose each item in the array.

---------------------------------------------
Choice             A       B     C      D    
----------------- ------- ----- ------ ------
 Observations     6       9     16     21  
---------------------------------------------

```{r}
choice_data <- c(6,9,16,21)
```

Let's compare two models:

* M1: Uniform categorical distribution. This captures the hypothesis that people are unbiased in their selections.
* M2: Dirichlet with $\alpha=2$ priors. This captures the hypothesis that there is bias in participants' choices. We aren't using Jeffreys priors because they attach substantial weight to extreme biases, which are unlikely.

We will now try to see which model is supported by the data. We'll ignore the fact that the answer is fairly obvious, and a null hypothesis significance test would probably be adequate here.

**Question 1**: Are these descriptive or process models?

**Question 2**: We have chosen these models for simplicity. If our goal is to detect left/right bias, can you think of changes that make one or both models more reasonable, without adding much complexity?

**Question 3**: What is the marginal log likelihood for M1?

**Question 4**: Finish implementing the importance sampler below.

```{r}

impSamp <- function(logTargD, # log target density
                    propR,    # a function to sample proposals
                    logPropD, # the log density for our proposal distribution
                    # Somewhat arbitarily, we are computing
                    # the expectation of the 4th param
                    ef=function(x) {x[,4]}) {
  nSamps <- 100000 # The more the better
  proposals <- propR(nSamps)
  pDens <- logPropD(proposals)
  unnP <- logTargD(proposals)
  w <- exp(unnP-pDens)
  print(paste("Expected value of target function:",
    sprintf("%2.3f",sum(w*ef(proposals))/sum(w))))
  print(paste("Log of average importance weight:",
    sprintf("%2.6f",log(sum(w)/nSamps))))
}
alpha <- 2

# Define logTargD, propR, logPropD. You can use ef to sanity-check --
#   it returns the posterior expectation 4th alpha entry, which you can
#   compute exactly: alpha_k/(sum_i alpha_i), combining real counts and
#   pseudocounts.
# Then call impSamp on these arguments.

# logTargD <- ??
# propR <- ??
# logPropD <- ??
# 
# impSamp(logTargD,propR,logPropD,function(x) {x[,4]})
```

**Question 5**: What is the Bayes factor favoring the more complex model? How should we interpret this?

**Question 6**: Which model is favored by the BIC? Does the BIC make sense to use, here?