<jupyter_start><jupyter_text>Question: How does Zipf's Exponent change with sample size?Our entire thesis is that Zipf's Law and Heap's Law are (excellent) approximations for describing text, but ultimately wrong. Let's face it, if they're wrong, it's not by much. So, if these polylogarithmic functions are correct, then what's the "apparent" Zipf exponent $\alpha(x)$, and is there a compelling derivation that matches the empirical curve?<jupyter_code># bloody dependencies
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = 16,8
from nltk.corpus import gutenberg
import numpy as np
import pandas as pd
from scipy.special import zeta, zetac
from scipy.optimize import fsolve

# custom classes
from legomena import Corpus, LogModel
def zetainv(h):
    func = lambda z: zeta(z) - h
    z0 = 1.1
    z = fsolve(func, z0)[0]
    return z

def hapax(x):
    if abs(x-1) < 0.001:
        return 0.5
    return 1./np.log(x) - 1./(x-1)<jupyter_output><empty_output><jupyter_text>Choosing a BookSelect a book for use as our fitting corpus.<jupyter_code># moby dick
words = gutenberg.words("melville-moby_dick.txt")
corpus = Corpus(words)
corpus.resolution = 2000
corpus.dimension = 20
TTR = corpus.TTR
model = LogModel().fit(TTR.m_tokens, TTR.n_types)
print("Optimum (M_z,N_z) =", model.params)

# plot TTR curve
plt.plot(TTR.m_tokens, TTR.n_types)
plt.title("Type-Token Relation")
plt.xlabel("tokens")
plt.ylabel("types")
plt.show()

# plot WFD
xmax = 50
i = np.arange(1, xmax)
df = corpus.fdist.head(xmax)
alpha = corpus.alpha
plt.bar(df.index, df.freq)
plt.plot(i, corpus.fdist.freq[1]/i**alpha, color = 'r')
plt.title("Word Frequency Distribution (alpha=%0.4f)" % alpha)
plt.xlabel("rank")
plt.ylabel("frequency")
plt.show()

# plot the k-vector
gamma = corpus.gamma
plt.bar(i, corpus.k[i])
plt.plot(i, corpus.k[1]/i**gamma, color = 'r')
plt.title("n-Legomena Distribution (gamma=%0.4f)" % gamma)
plt.xlabel("n = number of repeats")
plt.ylabel("words repeating n times")
plt.show()<jupyter_output>Optimum (M_z,N_z) = LogParams(M_z=173158, N_z=15848)<jupyter_text>How do alpha & gamma vary over sample size?<jupyter_code># get data
data = []
for _ in range(999):
    m_tokens = np.random.choice(range(corpus.M))
    minicorpus = corpus.sample(m_tokens)
    alpha, gamma = minicorpus.alpha, minicorpus.gamma
    data.append((m_tokens, alpha, gamma))
data = pd.DataFrame(data, columns = ["m_tokens", "alpha", "gamma"]).sort_values("m_tokens")

# plot it!
plt.plot(data.m_tokens, data.alpha, label="alpha")
plt.plot(data.m_tokens, data.gamma, label="gamma")
plt.axvline(x = model.M_z, color = 'r', label = "optimum")
plt.title("Alpha & Gamma Evolution")
plt.xlabel("tokens")
plt.legend()
plt.show()

# gamma ~ alpha relationship
plt.scatter(data.alpha, data.gamma, label="observed")
plt.scatter(data.alpha, 1 + 1/data.alpha, label="theory")
plt.title("Gamma ~ Alpha")
plt.xlabel("alpha")
plt.ylabel("gamma")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>What do Alpha & Gamma look like at the optimum sample size?<jupyter_code># get data
data = []
for _ in range(999):
    optimum = corpus.sample(model.M_z)
    alpha, gamma = optimum.alpha, optimum.gamma
    data.append((alpha, gamma))
data = pd.DataFrame(data, columns = ["alpha", "gamma"])

# alpha
plt.hist(data.alpha, bins=100)
plt.title("Optimum Alpha")
plt.show()

# gamma
plt.hist(data.gamma, bins=100)
plt.title("Optimum Gamma")
plt.show()

# gamma ~ alpha relationship
plt.scatter(data.alpha, data.gamma)
plt.title("Gamma ~ Alpha")
plt.xlabel("alpha")
plt.ylabel("gamma")
plt.show()
# word-frequency distribution
alpha = optimum.alpha
df = optimum.fdist
logrank = np.log(df.index.values)
logfreq = np.log(df.freq.values)
plt.plot(logrank, logfreq, label="observed")
plt.plot(logrank, max(logfreq) -alpha * logrank, label="theory")
plt.plot(logrank, max(logfreq) -1 * logrank, label="perfect")
plt.title("Word Frequency Distribution of Optimum, alpha = %0.4f" % alpha)
plt.xlabel("log(rank)") 
plt.ylabel("log(freq)")
plt.legend()
plt.show()

# n-legomena distribution
gamma = optimum.gamma
df = optimum.kdf
df.head()
logrank = np.log(df.index.values)
logfreq = np.log(df.freq.values)
plt.plot(logrank, logfreq, label="observed")
plt.plot(logrank, max(logfreq) -gamma * logrank, label="theory")
plt.plot(logrank, max(logfreq) -2 * logrank, label="perfect")
plt.title("n-Legomena Distribution of Optimum, gamma = %0.4f" % gamma)
plt.xlabel("log(rank)")
plt.ylabel("log(freq)")
plt.legend()
plt.show()
# manual
df = optimum.kdf
x = df.index.values
df["norm"] = df.freq / df.freq.sum()
df["poly"] = 1.0 / x / (x+1)
df.head()

y = df.norm.values

# method 1: curve_fit
def func_powerlaw(x, a, b):
    return a * x**b
from scipy.optimize import curve_fit
[a,b], _ = curve_fit(func_powerlaw, x, y)
curve_fit_ = func_powerlaw(x, a, b)
df["curve_fit"] = curve_fit_
print("curve_fit=", b)

# method 2: log linear regression
alpha = corpus._powerlaw(x, y)
print("linregress =", alpha)
linregress_ = y[0] * x**alpha
df["linregress"] = linregress_

# method 3: clauset
import powerlaw
dist = powerlaw.Fit(y).power_law
alpha = dist.alpha
print("clauset=", -alpha)
clauset_ = y[0] * x**-alpha
df["clauset"] = clauset_

# pretty picture
plt.scatter(x, y, label="observed")
plt.scatter(x, df["poly"], label="poly 1/n/(n+1)")
plt.loglog(x, curve_fit_, label="curve_fit")
plt.loglog(x, linregress_, label="linregress")
plt.loglog(x, clauset_, label="clauset")
plt.legend()
plt.show()

# table
display(df.head())
display(df.tail())
# calculate apparent gamma
TTR = corpus.TTR
k_x = model.predict_k(TTR.m_tokens, corpus.dimension)
E_x = model.predict(TTR.m_tokens)

# enrich TTR data
df = TTR
df["x"] = df.m_tokens / corpus.M
df["y_obs"] = df.n_types / corpus.N
df["y_exp"] = E_x / corpus.N
df["h_obs"] = df.lego_1 / df.n_types
df["h_exp"] = k_x[:,1] / E_x
df["g_exp"] = [ zetainv(1/h) for h in df["h_exp"] ]

for col in df.columns:
    if "lego_" in col:
        n = float(col.split("_")[1])
        df[col + "_p_obs"] = df[col] / df.n_types
        df[col + "_p_exp"] = df["h_exp"] / n**df["g_exp"]
        if n < 6:
            df[col + "_p_mod"] = k_x[:,int(n)] / E_x
display(df.head())

# draw a pretty plot
plt.plot(df["x"], df["y_exp"], label = "y_exp")
plt.plot(df["x"], df["h_exp"], label = "h_exp")
plt.plot(df["x"], df["g_exp"], label = "g_exp")
plt.scatter(df["x"], df["y_obs"], label = "y_obs")
plt.scatter(df["x"], df["h_obs"], label = "h_obs")
plt.scatter(df["x"], df["gamma"], label = "g_obs")
plt.axvline(x = 1, color = 'r', label = "x=1")
plt.xlabel("x = m/M_z")
plt.ylabel("γ = γ(x)")
plt.legend()
plt.title("Apparent Gamma")
plt.show()

# draw more pretty plots
for n in range(1, 6):
    plt.scatter(df["x"], df[f"lego_{n}_p_obs"], label = f"lego_{n}_p_obs")
    plt.plot(df["x"], df[f"lego_{n}_p_exp"], label = "zeta")
    plt.plot(df["x"], df[f"lego_{n}_p_mod"], label = "model")
    plt.axvline(x = 1, color = 'r', label = "x=1")
    plt.xlabel("x = m/M_z")
    plt.ylabel("P(freq = n) = p_n(x)")
    plt.legend()
    plt.title(f"N-Legomena Proportions (n = {n})")
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1794113505140121, 'fraction_numerical': 0.01197014504999296, 'mean_word_length': 4.85973597359736, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '438138', 'license': 'MIT', 'n_tokens_mistral': 2994, 'n_tokens_neox': 2787, 'n_words': 636, 'path': 'notebooks/pseudo-powerlaw.ipynb', 'repo_name': 'VictorDavis/legomena', 'size': '7101'}