<jupyter_start><jupyter_text>Simulating Language, Lab 4, Iterated Bayesian Learning This simulation features a replication of the Reali & Griffiths iterated learning model of the evolution of frequency distributions, and is built around the Bayesian model of frequency learning/regularisation from Lab 3.  This simulation allows you to explore the effects of learning bias on cultural evolution. But before we get onto the model itself, we need to talk about log probabilities.<jupyter_code>from math import log, exp<jupyter_output><empty_output><jupyter_text>Introduction to log probabilitiesIn the lectures I introduced Bayes’ Rule as a relationship between probabilities: the posterior is proportional to the product of the likelihood and the prior, and all three of these quantities are probabilities. Doing Bayesian models of learning therefore involves manipulating probabilities, numbers between 0 and 1. And some of these probabilities can be very small indeed, because they involve multiplying small numbers lots of times (consider, for instance, how small the probability is of getting 100 heads if you flip a fair coin 100 times: it’s 0.5 x 0.5 x 0.5 ... 100 times, or $0.5^{100}$ if you prefer. That’s a very small number.)As I said in the Lab 3 worksheet (go back and check if you need a refresher), working with small numbers 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). There’s no way your computer can exactly represent every possible real number, so what it does instead is store numbers as accurately as it can, which involves introducing small rounding errors.   Why do you need to care about this? Well, if you are dealing with very very small numbers (as you might do if you were doing a Bayesian model which involves learning from lots of data) then the rounding errors become a real factor - for big numbers the rounding errors are so small (relative to the big numers) we don’t really care, but for very small numbers, the rounding errors might be relatively big. Worse, sometimes the computer will round a very very small number to 0, which can generate unexpected and hard-to-predict errors in your code (e.g. if you try to divide something by a very very small number which gets rounded to 0).  The solution to this is to have the computer work not with probabilities, but with *log probabilities*: we take our probabilities, take the log of those numbers, then carry on as before. We can convert log probabilities back to probabilities with the `exp` function, the inverse of the `log` function.<jupyter_code>print(log(1))
print(log(0.1))
print(log(0.000001))
print(exp(log(0.5)))
print(exp(log(0.1)))<jupyter_output><empty_output><jupyter_text>As you can see from the code above, taking the log of a very small number turns it into a large negative number - these are still real numbers, so the computer still can’t represent them exactly, but in the log domain the rounding errors will be proportionately smaller for very small numbers and the rounding-to-0 problem won’t crop up. Then, if we want to see an actual probability, rather than a log probability, we can reverse this process, using the exp function, to get back raw probabilities. Jumping back and forth from logs can introduce rounding errors of its own, but it’s necessary to avoid the catastrophic rounding errors you can get if you just work with raw probabilities. Some basic arithmetic operations work a bit differently with logs. If you want to multiply two probabilities, you *add* their logarithms; if you want to divide one probability by another, you *subtract* the logarithm of one from another. And there is no direct equivalent of adding and subtracting in the log domain, which involves a little bit of fancy footwork in the code that you don’t have to worry about too much. The important thing is 1) to understand that the code is going to manipulate log probabilities and 2) this changes nothing conceptually, it’s just a matter of implementation.<jupyter_code>print(0.5 * 0.5)
print(exp(log(0.5) + log(0.5)))
print(0.5 / 0.5)
print(exp(log(0.5) - log(0.5)))<jupyter_output><empty_output><jupyter_text>On to the codeFirst, loading in the functions for random numbers and the beta distribution, plus some more that are specifically for doing stuff with logs. For instance, `logsumexp` (which we get from the `scipy.special` library allows us to do addition in the log domain (remember, just using the normal addition operator + with log probabilities is the equivalent of multiplying the raw probabilities).<jupyter_code>import random
from scipy.stats import beta
from scipy.special import logsumexp
from math import log, log1p, exp

%matplotlib inline
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg', 'pdf')<jupyter_output><empty_output><jupyter_text>The code starts with various bits and pieces which we need for working with logs and probability distributions. We define a function called `log_subtract` which allows us to do the equivalent of subtraction in the log domain (because if we just use normal subtraction, -, that’s equivalent to division). Then there are a couple of functions which we need for doing probabilistic sampling the log domain - `normalize_logprobs` will take a list of logs and normalise them for us (the equivalent of taking a list of pseudo-probabilities and rescaling them so they sum to 1, but in the log domain) and `log_roulette_wheel` takes a list of log probabilities and selects a random index from that list, with probability of any particular index being selected being given by its log probability. These are the log equivalents of `normalize_probs` and `roulette_wheel` from the Lab 3 code - they do exactly the same things, but with log probabilities rather than normal probabilities. As I said for the equivalent functions last week, *it is not important that you understand exactly how they work*.<jupyter_code>def log_subtract(x,y):
    return x + log1p(-exp(y - x))

def normalize_logprobs(logprobs):
    logtotal = logsumexp(logprobs) #calculates the summed log probabilities
    normedlogs = []
    for logp in logprobs:
        normedlogs.append(logp - logtotal) #normalise - subtracting in the log domain
                                           #is equivalent to dividing in the normal domain
    return normedlogs

def log_roulette_wheel(normedlogs):
    r=log(random.random()) #generate a random number between 0 and 1, then convert to log
    accumulator = normedlogs[0] 
    for i in range(len(normedlogs)):
        if r < accumulator:
            return i
        accumulator = logsumexp([accumulator, normedlogs[i + 1]])<jupyter_output><empty_output><jupyter_text>As in Lab 3, the main part of the code starts by laying out our hypothesis space, our grid of possible values of $\theta$ or `pW1`. We will need two grids here - one of normal probabilities (because they are easier to look at) and one of log probabilities (since we are going to work with log probabilities when we do our calculations).<jupyter_code>grid_granularity = 100
grid_increment = 1 / grid_granularity
 
# sets up the grid of possible probabilities to consider
possible_pW1 = []
for i in range(grid_granularity):
    possible_pW1.append(grid_increment / 2 + (grid_increment * i))

# sets up the grid of log probabilities
possible_logpW1 = []
for pW1 in possible_pW1:
    possible_logpW1.append(log(pW1))<jupyter_output><empty_output><jupyter_text>Have a look at the two grids (`possible_pW1` and `possible_logpW1`). Do they look like you expected? Next up come the various functions we need for Bayesian inference. Again, these are all the same as last week, but just set up to work with log probabilities rather than normal probabilities.   The prior<jupyter_code>def calculate_logprior(alpha):
    logprior = []
    for pW1 in possible_pW1:
        logprior.append(beta.logpdf(pW1, alpha, alpha)) 
    return normalize_logprobs(logprior)<jupyter_output><empty_output><jupyter_text>This is exactly the same as the `calculate_prior` function from last week: you pass in an alpha parameter which determines the shape of the prior, it returns a list of values, one for each hypothesis in the hypothesis space, but the returned values are log probabilities rather than normal probabilities.   - Plot some different prior probability distributions - for example, try typing `plt.plot(possible_pW1, calculate_logprior(0.1))` to see the prior log-probability distribution over various values of `pW1` for the `alpha=0.1` prior. It should be the same shape as you saw last week, but the values up the y axis will be different - they are now log probabilities rather than probabilities, but all you need to care about is the shape.  Likelihood and productionAs per the Lab 3 model, we are going to model data - sets 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. `loglikelihood` takes a log probability of word 1 being produced, and use that to calculate the probability of word 0 (which is 1 minus the probability of word 1). `logproduce` is the equivalent of `produce` last week, it generates some data based on a log-probability of producing word 1.<jupyter_code>def loglikelihood(data, logpW1):
    logpW0 = log_subtract(log(1), logpW1) #probability of w0 is 1-prob of w1
    logprobs = [logpW0, logpW1]
    loglikelihoods = []
    for d in data:
        loglikelihood_this_item = logprobs[d] #d will be either 0 or 1, 
                                              #so can use as index
        loglikelihoods.append(loglikelihood_this_item)
    return sum(loglikelihoods) #summing log probabilities = 
                               #multiply non-log probabilities
    
def logproduce(logpW1, n_productions):
    logpW0 = log_subtract(log(1), logpW1)
    logprobs = [logpW0, logpW1]
    data = []
    for p in range(n_productions):
        data.append(log_roulette_wheel(logprobs))
    return data<jupyter_output><empty_output><jupyter_text>- Test out the `logproduce` function - remember, you need to feed it a log probability, so decide on a probability for w1 and then convert it to log using the log function. Check that it works as expected based on what you did last week. - Next, check out the `loglikelihood` function - how does the likelihood of a set of data depend on the data and the probability of word 1? It should work in exactly the same way as last week. Remember that the likelihood function returns a log probability, so you can convert this to a probability using the `exp` function if you want.   Learning<jupyter_code>def logposterior(data, logprior):
    posterior_logprobs = []
    for i in range(len(possible_logpW1)):
        logpW1 = possible_logpW1[i] 
        logp_h = logprior[i] #prior probability of this pW1
        logp_d = loglikelihood(data, logpW1) #likelihood of data given this pW1
        posterior_logprobs.append(logp_h + logp_d) #adding logs = 
                                                   #multiplying non-logs
    return normalize_logprobs(posterior_logprobs) 
     
def loglearn(data,prior):
    posterior_logprobs = logposterior(data, prior)
    selected_index = log_roulette_wheel(posterior_logprobs)
    return possible_logpW1[selected_index]<jupyter_output><empty_output><jupyter_text>Now we have all the bits we need to calculate the posterior probability distribution, and therefore to do learning (by picking a hypothesis, a value of pW1, based on its posterior probability). You can verify that the posterior has the right sort of shape by doing something like this (which I adapted from the Lab 3 answer sheet): ```pythonmy_prior = calculate_logprior(1)my_data1 = [0] * 2 + [1] * 2my_data2 = [0] * 4 + [1] * 4my_data3 = [0] * 8 + [1] * 8print('data1 = ',my_data1)print('data2 = ',my_data2)print('data3 = ',my_data3)my_logposterior1 = logposterior(my_data1,my_prior)my_logposterior2 = logposterior(my_data2,my_prior)my_logposterior3 = logposterior(my_data3,my_prior)plt.plot(possible_pW1, my_logposterior1,'grey',label='p(theta|my_data1)')plt.plot(possible_pW1, my_logposterior2,'purple',label='p(theta|my_data2)')plt.plot(possible_pW1, my_logposterior3,'red',label='p(theta|my_data3)')plt.xlabel("theta")plt.ylabel("logp(theta|data)")plt.legend()```The numbers might look a bit different, because the y axis is logprobs, but the posterior should be highest where you expext it to be and low where you expect it to be low.  Iterated learningAt last, we have all the bits we need to do iterated learning: we can have a learner infer a value of `pW1` given some observed data, using the `loglearn` function, then we can have that individual produce data which another individual can learn from, using the `logproduce` function.<jupyter_code>def iterate(alpha, n_productions, starting_count_w1, generations):
    prior = calculate_logprior(alpha)
    pW1_accumulator = [] #we build up lists of our inferred pW1 at each generation
    data_accumulator = [] #and number of produced word1 at each generation
    data = [1] * starting_count_w1 + [0] * (n_productions - starting_count_w1)
    for generation in range(generations):
        logpW1 = loglearn(data, prior)
        data = logproduce(logpW1, n_productions)
        pW1_accumulator.append(exp(logpW1)) #will convert the inferred pW1 to a normal probability for you!
        data_accumulator.append(sum(data))
    return pW1_accumulator, data_accumulator<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06111804892382282, 'fraction_numerical': 0.014030705942848748, 'mean_word_length': 4.489516129032258, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '658230', 'license': 'MIT', 'n_tokens_mistral': 3872, 'n_tokens_neox': 3653, 'n_words': 1773, 'path': 'lab4.ipynb', 'repo_name': 'ponrawee/simlang2022', 'size': '13613'}