To answer the second question first, it\u2019s someone who supported chartism, which was \u201ca working-class movement for political reform in the United Kingdom that erupted from 1838 to 1857\u201d, quoth Wikipedia<\/a>. The name comes from the People\u2019s Charter of 1838, which called for voting rights for unpropertied men, among other reforms.<\/p>\n\n\n\n To answer the first question, we\u2019ll do some Bayesian statistics. My solution is based on a model that\u2019s not very realistic, so we should not take the result too seriously, but it demonstrates some interesting methods, I think. And as you\u2019ll see, there is a connection to Zipf\u2019s law, which I wrote about last week<\/a>.<\/p>\n\n\n\n Since last week\u2019s post was at the beginner level, I should warn you that this one is more advanced \u2013 in rapid succession, it involves the beta distribution, the t<\/em> distribution, the negative binomial, and the binomial.<\/p>\n\n\n\n This post is based on Think Bayes 2e<\/em>, which is available from Bookshop.org<\/a> and Amazon<\/a>.<\/p>\n\n\n\n Click here to run this notebook on Colab<\/a>. <\/p>\n\n\n\n If you don\u2019t hear a word for more than 50 years, that suggests it is not a common word. We can use Bayes\u2019s theorem to quantify this intuition. First we\u2019ll compute the posterior distribution of the word\u2019s frequency, then the posterior predictive distribution of hearing it again within a day.<\/p>\n\n\n\n Because we have only one piece of data \u2013 the time until first appearance \u2013 we\u2019ll need a good prior distribution. Which means we\u2019ll need a large, good quality sample of English text. For that, I\u2019ll use a free sample of the COCA dataset from CorpusData.org<\/a>. The following cells download and read the data.<\/p>\n\n\n\n We\u2019ll use a The dataset includes about 188,000 unique strings, but not all of them are what we would consider words.<\/p>\n\n\n\n To narrow it down, I\u2019ll remove anything that starts or ends with a non-alphabetical character \u2013 so hyphens and apostrophes are allowed in the middle of a word.<\/p>\n\n\n\n This filter reduces the number of unique words to about 151,000.<\/p>\n\n\n\n Of the 50 most common words, all of them have one syllable except number 38. Before you look at the list, can you guess the most common two-syllable word? Here\u2019s a theory about why common words are short<\/a>.<\/p>\n\n\n\n There are about 72,000 words that only appear once in the corpus, technically known as hapax legomena<\/a>.<\/p>\n\n\n\n Here\u2019s a random selection of them. Many are proper names, typos, or other non-words, but some are legitimate but rare words.<\/p>\n\n\n\n Now let\u2019s see what the distribution of word frequencies looks like.<\/p>\n\n\n\n One way to visualize the distribution is a Zipf plot, which shows the ranks on the x-axis and the frequencies on the y-axis.<\/p>\n\n\n\n Here\u2019s what it looks like on a log-log scale.<\/p>\n\n\n\n Zipf\u2019s law suggest that the result should be a straight line with slope close to -1. It\u2019s not exactly a straight line, but it\u2019s close, and the slope is about -1.1.<\/p>\n\n\n\n The Zipf plot is a well-known visual representation of the distribution of frequencies, but for the current problem, we\u2019ll switch to a different representation.<\/p>\n\n\n\n Given the number of times each word appear in the corpus, we can compute the rates, which is the number of times we expect each word to appear in a sample of a given size, and the inverse rates, which are the number of words we need to see before we expect a given word to appear.<\/p>\n\n\n\n We will find it most convenient to work with the distribution of inverse rates on a log scale. The first step is to use the observed frequencies to estimate word rates \u2013 we\u2019ll estimate the rate at which each word would appear in a random sample.<\/p>\n\n\n\n We\u2019ll do that by creating a beta distribution that represents the posterior distribution of word rates, given the observed frequencies (see this section of Think Bayes<\/em><\/a>) \u2013 and then drawing a random sample from the posterior. So words that have the same frequency will not generally have the same inferred rate.<\/p>\n\n\n\n Now we can compute the inverse rates, which are the number of words we have to sample before we expect to see each word once.<\/p>\n\n\n\n And here are their magnitudes, expressed as logarithms base 10.<\/p>\n\n\n\n To represent the distribution of these magnitudes, we\u2019ll use a Here\u2019s how we make the survival function.<\/p>\n\n\n\n And here\u2019s what it looks like.<\/p>\n\n\n\n The tail distribution has the sigmoid shape that is characteristic of normal distributions and t<\/em> distributions, although it is notably asymmetric.<\/p>\n\n\n\n And here\u2019s what the tail probabilities look like on a log-y scale.<\/p>\n\n\n\n If this distribution were normal, we would expect this curve to drop off with increasing slope. But for the words with the lowest frequencies \u2013 that is, the highest inverse rates \u2013 it is almost a straight line. And that suggests that a<\/p>\n\n\n\n distribution might be a good model for this data.<\/p>\n\n\n\n To estimate the frequency of rare words, we will need to model the tail behavior of this distribution and extrapolate it beyond the data. So let\u2019s fit a t<\/em> distribution and see how it looks. I\u2019ll use code from Chapter 8 of Probably Overthinking It<\/em><\/a>, which is all about these long-tailed distributions.<\/p>\n\n\n\n The following function makes a If we are given the But since we are not given Here\u2019s the With the y-axis on a linear scale, we can see that the model fits the data reasonably well, except for a range between 5 and 6 \u2013 that is for words that appear about 1 time in a million.<\/p>\n\n\n\n Here\u2019s what the model looks like on a log-y scale.<\/p>\n\n\n\n The model fits the data well in the extreme tail, which is exactly where we need it. And we can use the model to extrapolate a little beyond the data, to make sure we cover the range that will turn out to be likely in the scenario where we hear a word for this first time after 50 years.<\/p>\n\n\n\n The model we\u2019ve developed is the distribution of inverse rates for the words that appear in the corpus and, by extrapolation, for additional rare words that didn\u2019t appear in the corpus. This distribution will be the prior for the Bayesian update. We just have to convert it from a survival function to a PMF (remembering that these are equivalent representations of the same distribution).<\/p>\n\n\n\n To compute the likelihood of the observation, we have to transform the inverse rates to probabilities.<\/p>\n\n\n\n Now suppose that in a given day, you read or hear 10,000 words in a context where you would notice if you heard a word for the first time. Here\u2019s the number of words you would hear in 50 years.<\/p>\n\n\n\n Now, what\u2019s the probability that you fail to encounter a word in With this likelihood and the prior, we can compute the posterior distribution in the usual way.<\/p>\n\n\n\n And here\u2019s what it looks like.<\/p>\n\n\n\n If you go 50 years without hearing a word, that suggests that it is a rare word, and the posterior distribution reflects that logic.<\/p>\n\n\n\n The posterior distribution represents a range of possible values for the inverse rate of the word you heard. Now we can use it to answer the question we started with: what is the probability of hearing the same word again on the same day \u2013 that is, within the next 10,000 words you hear?<\/p>\n\n\n\n To answer that, we can use the survival function of the binomial distribution<\/a> to compute the probability of more than 0 successes in the next And we can use the probabilities in the posterior to compute the expected value \u2013 by the law of total probability, the result is the probability of hearing the same word again within a day.<\/p>\n\n\n\n The result is about 1 in 6000.<\/p>\n\n\n\n With all of the assumptions we made in this calculation, there\u2019s no reason to be more precise than that. And as I mentioned at the beginning, we should probably not take this conclusion too seriously. For one thing, it’s likely that my experience is an example of the frequency illusion<\/a>, which is “a cognitive bias in which a person notices a specific concept, word, or product more frequently after recently becoming aware of it.” Also, if you hear a word for the first time after 50 years, there\u2019s a good chance the word is having a moment, which greatly increases the chance you\u2019ll hear it again. I can\u2019t think of why chartism might be in the news at the moment, but maybe this post will go viral and make it happen.<\/p>\n","protected":false},"excerpt":{"rendered":" Recently I heard the word \u201cchartist\u201d for the first time in my life (that I recall). And then later the same day, I heard it again. So that raises two questions: To answer the second question first, it\u2019s someone who supported chartism, which was \u201ca working-class movement for political reform in the United Kingdom that erupted from 1838 to 1857\u201d, quoth Wikipedia. The name comes from the People\u2019s Charter of 1838, which called for voting rights for unpropertied men, among…<\/p>\nWord Frequencies<\/h2>\n\n\n\n
download(\"https:\/\/www.corpusdata.org\/coca\/samples\/coca-samples-text.zip\")\n<\/pre>\n\n\n\n
import zipfile\n\n\ndef generate_lines(zip_path=\"coca-samples-text.zip\"):\n with zipfile.ZipFile(zip_path, \"r\") as zip_file:\n file_list = zip_file.namelist()\n for file_name in file_list:\n with zip_file.open(file_name) as file:\n lines = file.readlines()\n for line in lines:\n yield (line.decode(\"utf-8\"))\n<\/pre>\n\n\n\n
Counter<\/code> to count the number of times each word appears.<\/p>\n\n\n\nimport re\nfrom collections import Counter\n\npattern = r\"[ \/\\n]+|--\"\n\ncounter = Counter()\nfor line in generate_lines():\n words = re.split(pattern, line)[1:]\n counter.update(word.lower() for word in words if word)\n<\/pre>\n\n\n\n
len(counter), counter.total()\n<\/pre>\n\n\n\n
(188086, 11503819)\n<\/pre>\n\n\n\n
for s in list(counter.keys()):\n if not s[0].isalpha() or not s[-1].isalpha():\n del counter[s]\n<\/pre>\n\n\n\n
num_words = counter.total()\nlen(counter), num_words\n<\/pre>\n\n\n\n
(151414, 8889694)\n<\/pre>\n\n\n\n
for i, (word, freq) in enumerate(counter.most_common(50)):\n print(f'{i+1}\\t{word}\\t{freq}')\n<\/pre>\n\n\n\n1\tthe\t461991\n2\tto\t237929\n3\tand\t231459\n4\tof\t217363\n5\ta\t203302\n6\tin\t153323\n7\ti\t137931\n8\tthat\t123818\n9\tyou\t109635\n10\tit\t103712\n11\tis\t93996\n12\tfor\t78755\n13\ton\t64869\n14\twas\t64388\n15\twith\t59724\n16\the\t57684\n17\tthis\t51879\n18\tas\t51202\n19\tn't\t49291\n20\twe\t47694\n21\tare\t47192\n22\thave\t46963\n23\tbe\t46563\n24\tnot\t43872\n25\tbut\t42434\n26\tthey\t42411\n27\tat\t42017\n28\tdo\t41568\n29\twhat\t35637\n30\tfrom\t34557\n31\this\t33578\n32\tby\t32583\n33\tor\t32146\n34\tshe\t29945\n35\tall\t29391\n36\tmy\t29390\n37\tan\t28580\n38\tabout\t27804\n39\tthere\t27291\n40\tso\t27081\n41\ther\t26363\n42\tone\t26022\n43\thad\t25656\n44\tif\t25373\n45\tyour\t24641\n46\tme\t24551\n47\twho\t23500\n48\tcan\t23311\n49\ttheir\t23221\n50\tout\t22902\n<\/pre>\n\n\n\n
singletons = [word for (word, freq) in counter.items() if freq == 1]\nlen(singletons), len(singletons) \/ counter.total() * 100\n<\/pre>\n\n\n\n
(72159, 0.811715228893143)\n<\/pre>\n\n\n\n
np.random.choice(singletons, 100)\n<\/pre>\n\n\n\n
array(['laneer', 'emc', 'literature-like', 'tomyworld', 'roald',\n 'unreleased', 'basemen', 'kielhau', 'clobber', 'feydeau',\n 'symptomless', 'channelmaster', 'v-i', 'tipsha', 'mjlkdroppen',\n 'harlots', 'phaetons', 'grlinger', 'naniwa', 'dadian',\n 'banafionen', 'ceramaseal', 'vine-covered', 'terrafirmahome.com',\n 'hesten', 'undertheorized', 'fantastycznie', 'kaido', 'noughts',\n 'hannelie', 'cacoa', 'subelement', 'mestothelioma', 'gut-level',\n 'abis', 'potterville', 'quarter-to-quarter', 'lokkii', 'telemed',\n 'whitewood', 'dualmode', 'plebiscites', 'loubrutton',\n 'off-loading', 'abbot-t-t', 'whackaloons', 'tuinal', 'guyi',\n 'samanthalaughs', 'editor-sponsored', 'neurosciences', 'lunched',\n 'chicken-and-brisket', 'korekane', 'ruby-colored',\n 'double-elimination', 'cornhusker', 'wjounds', 'mendy', 'red.ooh',\n 'delighters', 'tuviera', 'spot-lit', 'tuskarr', 'easy-many',\n 'timepoint', 'mouthfuls', 'catchy-titled', 'b.l', 'four-ply',\n \"sa'ud\", 'millenarianism', 'gelder', 'cinnam',\n 'documentary-filmmaking', 'huviesen', 'by-gone', 'boy-friend',\n 'heartlight', 'farecompare.com', 'nurya', 'overstaying',\n 'johnny-turn', 'rashness', 'mestier', 'trivedi', 'koshanska',\n 'tremulousness', 'movies-another', 'womenfolks', 'bawdy',\n 'all-her-life', 'lakhani', 'screeeeaming', 'marketings', 'girthy',\n 'non-discriminatory', 'chumpy', 'resque', 'lysing'], dtype='<U24')\n<\/pre>\n\n\n\n
Zipf\u2019s Law<\/h2>\n\n\n\n
freqs = np.array(sorted(counter.values(), reverse=True))\n<\/pre>\n\n\n\n
n = len(freqs)\nranks = range(1, n + 1)\n<\/pre>\n\n\n\n
plt.plot(ranks, freqs)\n\ndecorate(\n title=\"Zipf plot\", xlabel=\"Rank\", ylabel=\"Frequency\", xscale=\"log\", yscale=\"log\"\n)\n<\/pre>\n\n\n\n
![\"_images\/204f4aae3fe537fefdbe43abadd4be2a854bc627c7d5e064d9efebd9cc6a58df.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/204f4aae3fe537fefdbe43abadd4be2a854bc627c7d5e064d9efebd9cc6a58df.png\")
<\/figure>\n\n\n\nrise = np.log10(freqs[-1]) - np.log10(freqs[0])\nrise\n<\/pre>\n\n\n\n
-5.664633515191604\n<\/pre>\n\n\n\n
run = np.log10(ranks[-1]) - np.log10(ranks[0])\nrun\n<\/pre>\n\n\n\n
5.180166032638616\n<\/pre>\n\n\n\n
rise \/ run\n<\/pre>\n\n\n\n
-1.0935235433575892\n<\/pre>\n\n\n\n
Tail Distribution<\/h2>\n\n\n\n
from scipy.stats import beta\n\nnp.random.seed(17)\nalphas = freqs + 1\nbetas = num_words - freqs + 1\ninferred_rates = beta(alphas, betas).rvs()\n<\/pre>\n\n\n\n
inverse_rates = 1 \/ inferred_rates\n<\/pre>\n\n\n\n
mags = np.log10(inverse_rates)\n<\/pre>\n\n\n\n
Surv<\/code> object, which represents survival functions, but we\u2019ll use a variation of the survival function which is the probability that a randomly-chosen value is greater than or equal to a given quantity. The following function computes this version of a survival function, which is called a tail probability.<\/p>\n\n\n\nfrom empiricaldist import Surv\n\n\ndef make_surv(seq):\n \"\"\"Make a non-standard survival function, P(X>=x)\"\"\"\n pmf = Pmf.from_seq(seq)\n surv = pmf.make_surv() + pmf\n\n # correct for numerical error\n surv.iloc[0] = 1\n return Surv(surv)\n<\/pre>\n\n\n\n
surv = make_surv(mags)\n<\/pre>\n\n\n\n
options = dict(marker=\".\", ms=2, lw=0.5, label=\"data\")\nsurv.plot(**options)\ndecorate(xlabel=\"Inverse rate (log10 words per appearance)\", ylabel=\"Tail probability\")\n<\/pre>\n\n\n\n
![\"_images\/a85267e340ee9e92988ee9ce9ec80c2edcf7c526f04c09a28b53b4ca43f0fa1a.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/a85267e340ee9e92988ee9ce9ec80c2edcf7c526f04c09a28b53b4ca43f0fa1a.png\")
<\/figure>\n\n\n\nsurv.plot(**options)\ndecorate(xlabel=\"Inverse rate (words per appearance)\", yscale=\"log\")\n<\/pre>\n\n\n\n
![\"_images\/ad17e7857447f99903d3a718f91bcaae09fdcf067df4d4f48f5127d3e8151c5d.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/ad17e7857447f99903d3a718f91bcaae09fdcf067df4d4f48f5127d3e8151c5d.png\")
<\/figure>\n\n\n\nFitting a Model<\/h2>\n\n\n\n
Surv<\/code> object that represents a t<\/em> distribution with the given parameters.<\/p>\n\n\n\nfrom scipy.stats import t as t_dist\n\n\ndef truncated_t_sf(qs, df, mu, sigma):\n \"\"\"Makes Surv object for a t distribution.\n \n Truncated on the left, assuming all values are greater than min(qs)\n \"\"\"\n ps = t_dist.sf(qs, df, mu, sigma)\n surv_model = Surv(ps \/ ps[0], qs)\n return surv_model\n<\/pre>\n\n\n\n
df<\/code> parameter, we can use the following function to find the values of mu<\/code> and sigma<\/code> that best fit the data, focusing on the central part of the distribution.<\/p>\n\n\n\nfrom scipy.optimize import least_squares\n\n\ndef fit_truncated_t(df, surv):\n \"\"\"Given df, find the best values of mu and sigma.\"\"\"\n low, high = surv.qs.min(), surv.qs.max()\n qs_model = np.linspace(low, high, 2000)\n ps = np.linspace(0.1, 0.8, 20)\n qs = surv.inverse(ps)\n\n def error_func_t(params, df, surv):\n mu, sigma = params\n surv_model = truncated_t_sf(qs_model, df, mu, sigma)\n\n error = surv(qs) - surv_model(qs)\n return error\n\n pmf = surv.make_pmf()\n pmf.normalize()\n params = pmf.mean(), pmf.std()\n res = least_squares(error_func_t, x0=params, args=(df, surv), xtol=1e-3)\n assert res.success\n return res.x\n<\/pre>\n\n\n\n
df<\/code>, we can use the following function to search for the value that best fits the tail of the distribution.<\/p>\n\n\n\nfrom scipy.optimize import minimize
def minimize_df(df0, surv, bounds=[(1, 1e3)], ps=None):
low, high = surv.qs.min(), surv.qs.max()
qs_model = np.linspace(low, high * 1.2, 2000)
if ps is None:
t = surv.ps[0], surv.ps[-5]
low, high = np.log10(t)
ps = np.logspace(low, high, 30, endpoint=False)
qs = surv.inverse(ps)
def error_func_tail(params):
(df,) = params
# print(df)
mu, sigma = fit_truncated_t(df, surv)
surv_model = truncated_t_sf(qs_model, df, mu, sigma)
errors = np.log10(surv(qs)) - np.log10(surv_model(qs))
return np.sum(errors**2)
params = (df0,)
res = minimize(error_func_tail, x0=params, bounds=bounds, method=\"Powell\")
assert res.success
return res.x
<\/pre>\n\n\n\ndf = minimize_df(25, surv)\ndf\n<\/pre>\n\n\n\n
array([22.52401171])\n<\/pre>\n\n\n\n
mu, sigma = fit_truncated_t(df, surv)\ndf, mu, sigma\n<\/pre>\n\n\n\n
(array([22.52401171]), 6.433323515095857, 0.49070837962997577)\n<\/pre>\n\n\n\n
t<\/code> distribution that best fits the data.<\/p>\n\n\n\nlow, high = surv.qs.min(), surv.qs.max()\nqs = np.linspace(low, 10, 2000)\nsurv_model = truncated_t_sf(qs, df, mu, sigma)\n<\/pre>\n\n\n\n
surv_model.plot(color=\"gray\", alpha=0.4, label=\"model\")\nsurv.plot(**options)\ndecorate(xlabel=\"Inverse rate (log10 words per appearance)\", ylabel=\"Tail probability\")\n<\/pre>\n\n\n\n
![\"_images\/6a78db0b7207d7492f70ad6ec717b9441d528c6bae9ad44e8e5673a7982f4777.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/6a78db0b7207d7492f70ad6ec717b9441d528c6bae9ad44e8e5673a7982f4777.png\")
<\/figure>\n\n\n\nsurv_model.plot(color=\"gray\", alpha=0.4, label=\"model\")\nsurv.plot(**options)\ndecorate(\n xlabel=\"Inverse rate (log10 words per appearance)\",\n ylabel=\"Tail probability\",\n yscale=\"log\",\n)\n<\/pre>\n\n\n\n
![\"_images\/8486cbced76150d86f5639eb3320a68b34fd31325bd8d20f9caf3d994dd0c669.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/8486cbced76150d86f5639eb3320a68b34fd31325bd8d20f9caf3d994dd0c669.png\")
<\/figure>\n\n\n\nThe Update<\/h2>\n\n\n\n
prior = surv_model.make_pmf()\nprior.plot(label=\"prior\")\ndecorate(\n xlabel=\"Inverse rate (log10 words per appearance)\",\n ylabel=\"Density\",\n)\n<\/pre>\n\n\n\n
![\"_images\/baed114476768c8511a77423e7217807b462ce1d867fe098a75626df9ed59042.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/baed114476768c8511a77423e7217807b462ce1d867fe098a75626df9ed59042.png\")
<\/figure>\n\n\n\nps = 1 \/ np.power(10, prior.qs)\n<\/pre>\n\n\n\n
words_per_day = 10_000\ndays = 50 * 365\nk = days * words_per_day\nk\n<\/pre>\n\n\n\n
182500000\n<\/pre>\n\n\n\n
k<\/code> attempts and then encounter it on the next attempt? We can answer that with the negative binomial distribution, which computes the probability of getting the n<\/code>th success after k<\/code> failures, for a given probability \u2013 or in this case, for a sequence of possible probabilities.<\/p>\n\n\n\nfrom scipy.stats import nbinom\n\nn = 1\nlikelihood = nbinom.pmf(k, n, ps)\n<\/pre>\n\n\n\n
posterior = prior * likelihood\nposterior.normalize()\n<\/pre>\n\n\n\n
1.368245917258196e-11\n<\/pre>\n\n\n\n
prior.plot(alpha=0.5, label=\"prior\")\nposterior.plot(label=\"posterior\")\ndecorate(\n xlabel=\"Inverse rate (log10 words per appearance)\",\n ylabel=\"Density\",\n)\n<\/pre>\n\n\n\n
![\"_images\/9ee20cb48e5268c6716d40b69fd8eede6f4a244ae418222ecfaaa83a632dd371.png\"\/](\"https:\/\/allendowney.github.io\/ThinkBayes2\/_images\/9ee20cb48e5268c6716d40b69fd8eede6f4a244ae418222ecfaaa83a632dd371.png\")
<\/figure>\n\n\n\nn_pred<\/code> attempts. We\u2019ll compute this probability for each of the ps<\/code> that correspond to the inverse rates in the posterior.<\/p>\n\n\n\nfrom scipy.stats import binom\n\nn_pred = words_per_day\nps_pred = binom.sf(0, n_pred, ps)\n<\/pre>\n\n\n\n
p = np.sum(posterior * ps_pred)\np, 1 \/ p\n<\/pre>\n\n\n\n
(0.00016019406802217392, 6242.42840166579)\n<\/pre>\n\n\n\n