lab3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating Language, Lab 3, Regularisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This week we'll be working with a simple Bayesian model of frequency learning.  This simulation allows you to explore the effects of the prior and the data on frequency learning, as discussed in the lecture as a model of the Hudson-Kam & Newport (2005) experiment; next week we'll use the same model to move beyond studying individuals and start looking at cultural evolution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The code\n",
    "\n",
    "The basic framework is the same as the word learning code: after some preliminaries we lay out a hypothesis space, specify the likelihood and the prior, then we have everything we need to calculate the posterior and do Bayesian inference. All the details are different from the word learning lab, because we are modelling a different aspect of language learning, but the code will hopefully already look somewhat familiar.\n",
    "\n",
    "### Libraries etc\n",
    "First, we'll load the `random` library (for generating random numbers) and the `prod` function (for multiplying a list of numbers), plus one more library for doing stuff with beta distributions, which we are using for our prior. We also have to load the plotting library and set up inline plots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "from numpy import prod\n",
    "from scipy.stats import beta\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import set_matplotlib_formats\n",
    "set_matplotlib_formats('svg', 'pdf')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some useful functions for dealing with probabilities\n",
    "\n",
    "The code starts with two functions we need for doing stuff with probabilities. You saw `normalize_probs` last week, it will take a list of numbers and normalise them for us (i.e. scaling them so they sum to 1). `roulette_wheel` takes a list of probabilities and selects a random index from that list, with probability of any particular index being selected being given by its probability (i.e. if index 0 has twice the probability as index 1, it's twice as likely to be selected). These functions are used elsewhere in the code, but it is not particularly important that you understand exactly how they work. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize_probs(probs):\n",
    "    total = sum(probs) #calculates the summed probabilities\n",
    "    normedprobs = []\n",
    "    for p in probs:\n",
    "        normedprobs.append(p / total) #normalise - divide by summed probs\n",
    "    return normedprobs\n",
    "\n",
    "def roulette_wheel(normedprobs):\n",
    "    r=random.random() #generate a random number between 0 and 1\n",
    "    accumulator = normedprobs[0]\n",
    "    for i in range(len(normedprobs)):\n",
    "        if r < accumulator:\n",
    "            return i\n",
    "        accumulator = accumulator + normedprobs[i + 1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The hypothesis space\n",
    "The main part of the code starts by laying out our hypothesis space, the candidate hypotheses that our learner is going to consider. As discussed in class, we are going to turn the problem of inferring a potentially continuous value (the probability with which your teacher uses word 1) into the problem of inferring one of a limited set of possible values (either your teacher is using the word with probability 0.005, or 0.015, or 0.025, etc). In the code we will refer to a certain probability of using word 1 as `pW1`.  We will call the set of possible values for `pW1` the grid - you can set the granularity of the grid as high as you like, but 100 works OK without being too slow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_granularity = 100 #How many values of pW1 do you want to consider? Setting this to 100.\n",
    "grid_increment = 1 / grid_granularity\n",
    " \n",
    "# sets up the grid of possible values of pW1 to consider\n",
    "possible_pW1 = []\n",
    "for i in range(grid_granularity):\n",
    "    possible_pW1.append(grid_increment / 2 + (grid_increment * i))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have a look at `possible_pW1`. Does it look like you expected?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One thing you might notice (you might already have noticed it last week), and be a bit surprised about, is that there are rounding errors! Instead of the probabilities in the grid being exactly 0.05 apart, sometimes the values are over or under by a tiny amount. Working with probabilities on a computer can be a problem, because the computer cannot exactly represent real numbers (i.e. numbers we would write in decimal notation, e.g. numbers like 0.1, 3.147). Your computer has a very large memory where it can store and manipulate numbers, but the problem is that this memory is necessarily finite (it has to fit in your computer) and there are infinitely many real numbers. Think of recurring decimal you get by dividing 1 by 3, 0.3333..., where the threes go on forever - it would take an infinite amount of space to exactly represent this number in your computer, and distinguish it from a very similar number, e.g. 0.33333... where the threes go on for a few thousand repetitions only. More relevantly for the rounding errors above, imagine how much memory it would take to distinguish 0.5 from 0.50000000000...00001, where there could be arbitrarily many decimal places. Spoiler: it would take an infinite memory to do this perfectly. So there’s no way your computer can exactly represent every possible real number. What it does instead is store numbers as accurately as it can, which involves introducing small rounding errors, rounding a number it can't represent to a number it can (and sometimes the results there are slightly un-intuitive to us, like rounding 0.035 to 0.034999999999999996). In fact your computer does its best to conceal these errors from you, and often displays numbers in a format that hides exactly what numbers it is actually working with. But you can see those rounding errors here. Next week we are going to be forced to introduce a technique to deal with these rounding errors otherwise they start causing glitches in the code, but for this week we'll try to live with them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next up come the various functions we need for Bayesian inference. I will step through these gradually."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The prior\n",
    "For this model our prior is more complicated than last week - we are modelling the prior using a *beta distribution*, which is a family of probability distributions that can capture a uniform prior (representing an unbiased learner), a prior favouring regularity, or a prior favouring variability. \n",
    "\n",
    "The `calculate_prior` function calculates the prior probability of each of our possible values of `pW1`. The beta distribution, which is what we are using for our prior, is a standard probability distribution, so we can just use a function from a library (`beta.pdf`) to get the probability density for each value of `pW1`, then normalise those to convert them to probabilities. The `alpha` parameter determines the shape of our prior, and therefore the bias of our learners when it comes to inferring `pW1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_prior(alpha):\n",
    "    prior = []\n",
    "    for pW1 in possible_pW1:\n",
    "        prior.append(beta.pdf(pW1, alpha, alpha)) #look up the value using beta.pdf and add to our growing prior\n",
    "    return normalize_probs(prior) #normalize the final list so they are all probabilities\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot some different prior probability distributions. To get a line graph, try typing e.g. \n",
    "```python\n",
    "plt.plot(possible_pW1, calculate_prior(0.1))\n",
    "``` \n",
    "to see the prior probability distribution over various values of pW1 for the alpha=0.1 prior. Or, if you prefer a barplot (one bar per value of pW1) do \n",
    "```python\n",
    "plt.bar(possible_pW1,calculate_prior(0.1),align='center',width=1./grid_granularity)\n",
    "``` \n",
    "(`align='center'` and `width=1./grid_granularity` makes sure your x-axis looks right and the bars line up where they should). \n",
    "\n",
    "Play around with the code, trying different values for the alpha parameter, to answer these three questions: \n",
    "- What values of alpha lead to a prior bias for regularity (i.e. higher prior probability for pW1 closer to 0 or 1)? \n",
    "- What values of alpha lead to a prior bias for variability (i.e. higher prior probability for pW1 closer to 0.5)? \n",
    "- What values of alpha lead to a completely unbiased learner (i.e. all values of pW1 are a priori equally likely)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood and production\n",
    "\n",
    "In order to do Bayesian inference, we need a likelihood function that tells us how probable some data is given a certain hypothesis (a value of `pW1`) - that's our `likelihood` function. Next week we are also going to need a way of modelling production - taking an individual, with a value of `pW1` in their head, and having them produce data that someone else can learn from. We'll specify that now too, it's called `produce`. \n",
    "\n",
    "We are going to model data - a sequence of utterances - as a simple list of 0s and 1s: the 0s correspond to occurrences of word 0, the 1s correspond to occurrences of word 1. For example, this is how we will represent a data set consisting of one occurence of word 0 and one occurence of word 1:\n",
    "```python\n",
    "small_data = [0,1]\n",
    "```\n",
    "\n",
    "How would you represent a dataset consisting of two occurences of word 0 and two occurences of word 1?\n",
    "\n",
    "Both the `likelihood` function and the `produce` function take as an argument the probability of word 1 being produced, `pW1`, and use that to calculate the probability of word 0 being produced (which is `1 - pW1`: in this model there are only two words, every time you produce an utterance you produce one or the other). The `likelihood` function calculates the likelihood of `data`, a list of utterances, given a particular value of `pW1`. The `produce` function generates some data for a speaker with a specific value of `pW1` in their head - you tell it how many productions you want (`n_productions`) and it spits out a list of data for you. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(data, pW1):\n",
    "    pW0 = 1 - pW1 #probability of w0 is 1-prob of w1\n",
    "    probs = [pW0, pW1]\n",
    "    likelihoods = []\n",
    "    for d in data:\n",
    "        likelihood_this_item = probs[d] #d will be either 0 or 1, so we can use as an index\n",
    "        likelihoods.append(likelihood_this_item)\n",
    "    return prod(likelihoods) #multiply the probabilities of the individual data items\n",
    "    \n",
    "def produce(pW1, n_productions):\n",
    "    pW0 = 1- pW1\n",
    "    probs = [pW0, pW1]\n",
    "    data = []\n",
    "    for p in range(n_productions):\n",
    "        data.append(roulette_wheel(probs))\n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Test out the `produce` function - decide on a probability for w1 and then specify how many utterances you would like to produce. What kind of data will be produced if the probability of w1 is low? Or if it is high? Hint: `produce(0.1,10)` will produce 10 utterances with the probability of producing word 1 on each utterance being 0.1. \n",
    "- Next, check out the likelihood function - how does the likelihood of a set of data depend on the data and the probability of word 1? Hint: `likelihood([0,0,1,1],0.5)` will tell you the likelihood of producing word 0 twice then word 1 twice when the probability of producing word 1 each time is 0.5. Try plugging in different sequences of data and different values for `pW1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have all the bits we need to calculate the posterior probability distribution - that's what the `posterior` function does, and it works in exactly the same way as the `posterior` function from lab 2. We are also going to make a function, with we will call `learn`, which picks a hypothesis (a value of pW1) based on its posterior probability - again, we'll be needing this next week, but you can play with it now. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def posterior(data, prior):\n",
    "    posterior_probs = [] #this list will hold the posterior for each possible value of pW1\n",
    "    for i in range(len(possible_pW1)): #work through the list of pW1 values, by index\n",
    "        pW1 = possible_pW1[i] #look up that value of pW1\n",
    "        p_h = prior[i] #look up the prior probability of this pW1\n",
    "        p_d = likelihood(data, pW1) #calculate the likelihood of data given this pW1\n",
    "        p_h_given_d = p_h * p_d #multiply prior x likelihood\n",
    "        posterior_probs.append(p_h_given_d) \n",
    "    return normalize_probs(posterior_probs) #normalise\n",
    "     \n",
    "def learn(data,prior):\n",
    "    posterior_probs = posterior(data, prior) #calculate the posterior\n",
    "    selected_index = roulette_wheel(posterior_probs) #select a random index from the posterior\n",
    "    return possible_pW1[selected_index] #look up the corresponding value of pW1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `posterior` function implements Bayes Rule for us - the posterior probability of each value of `pW1` given the data is the likelihood of the data given that value of `pW1` multiplied by the prior probability of that value of `pW1`. \n",
    "\n",
    "Test out the `posterior` function. To do this you will need to build a prior, and make some data for your learner to learn from; you might also want to plot the posterior rather than staring at a big list of numbers. For example:\n",
    "```python\n",
    "my_prior = calculate_prior(1)\n",
    "my_data = [0] * 2 + [1] * 2\n",
    "print('data = ',my_data)\n",
    "my_posterior = posterior(my_data,my_prior)\n",
    "plt.plot(possible_pW1, my_posterior)\n",
    "```\n",
    "\n",
    "This example shows how to test a uniform prior (alpha = 1) and data consisting of two 1s and two 0s (note that there is a cute little trick there for creating lists of duplicates and sticking two lists together). Then it calculates the posterior probability for all possible values of `pW1` and plots it (the x-axis in the graph shows value of `pW1`, the y-axis shows the posterior probability of that value of `pW1` given `my_data`). \n",
    "\n",
    "Start with a uniform prior and see how the data affects the learner’s hypothesis. \n",
    "\n",
    "- What does adding more data do? \n",
    "- What does making the data highly skewed in favour of one word do? \n",
    "- Try different priors - what does a strong prior in favour of regularity do? \n",
    "- What does a strong prior in favour of variability do? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also play about with the `learn` function, which picks a single value of `pW1` based on the posterior. e.g. try `learn(my_data,my_prior)` several times. Do you always get the same answer? Why not?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions\n",
    "\n",
    "The priority for this worksheet is to work through the in-text questions above: experimenting with different shapes for the prior, checking that the `likelihood` and `production` functions makes sense to you, checking you understand how the posterior depends on the prior and the data and that you understand why`learn` behaves as it does. Once you are happy with that, try these questions:\n",
    "\n",
    "1. Can you produce a result like the Hudson Kam & Newport (2005) results *for adults*, i.e. that adult learners fairly accurately track the frequency of a linguistic variant in their input? What kinds of priors and what kinds of data does this work for? *Be aware: when using large amounts of data, you might run into the kinds of rounding error I mentioned earlier!* \n",
    "2. Can you produce a result like the Hudson Kam & Newport (2005) results *for children*, i.e. that children tend to regularise, sometimes producing only one variant even when their data contains variation? Again, what kinds of priors and what kinds of data does this work for?\n",
    "3. Do we have any kind *innateness* in our model? Are there languages our model learner could never learn, no matter what kind of data we gave them? Are there languages they are predisposed to learn, or languages that they are predisposed to not learn?\n",
    "4. [Harder, mainly a coding problem] It turns out that one reason this type of model (known as the beta-binomial model) is so popular is that there is a handy shortcut for working out the posterior: if your prior is beta($\\alpha, \\alpha$) and you see $a$ occurences of word 0 and $b$ occurences of word 1 then your posterior distribution will also be a beta distribution, beta($\\alpha + b$, $\\alpha + a$). We won't be using this shortcut because we'll be moving on from beta-binomial models shortly and our grid technique works for all the examples we are going to be dealing with, but for fun: can you use the `beta.pdf` function to write a new version of the `posterior` function? And can you use the function`beta.rvs(x,y)` (which samples a random number from a beta distribution with parameters `x` and `y`) to replace our `learn` function above?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "Hudson Kam, C. L., & Newport, E. L. (2005). Regularizing unpredictable variation: The roles of adult and child learners in language formation and change. *Language Learning and Development, 1,* 151–195."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}