# Sampling uniformly from the language generated by a grammar

Suppose I have the following formal grammar:

$$S \rightarrow \varepsilon + a S + b S + c S d S$$

where $S$ is a nonterminal symbol and $a$, $b$, $c$, $d$ are terminal symbols. How can I sample uniformly from the set of $n$-length strings that are generated by this grammar, without having to generate all $n$-length strings first?

For example, if the grammar were just

$$S \rightarrow \varepsilon + a S + bS$$

we could generate a random $n$-length string by choosing $\varepsilon$ if $n = 0$, and otherwise choosing between $a$ and $b$ with equal probability followed by a random $(n-1)$-length string.

• At least superficially, skimming title and abstract: relevant? www.sciencedirect.com/science/article/pii/S0890540197926213 – Clement C. Oct 23 '17 at 1:41

I managed to solve the problem with an exact procedure. First, we find the generating function for the grammar:

\begin{align} y &= 1 + x y + x y + x y x y \\ y &= \frac{1 - 2x - \sqrt{1 - 4x}}{2x^2} \\ &= 1 + 2x + 5x^2 + 14x^3 + 42x^4 + 132x^5 + 429x^6 + 1430x^7 + \ldots \end{align}

The coefficient $y_n$ of $x^n$ tells us how many strings of length $n$ there are in the grammar. In this case, the coefficients are exactly the Catalan numbers, shifted by one. They satisfy the recurrence relation

$$y_n = \begin{cases} 1 & n = 0 \\ \frac{2(2n+1)}{n+2} y_{n-1} & \text{otherwise} \end{cases}$$

We can pre-compute $y_n$ for all $n$ up to the maximum length we will sample. Next, we use the technique for random generation of combinatorial structures described in section 3 ("basic generation schemes") of A calculus for the random generation of labelled combinatorial structures by Flajolet et al:

Adapted to Python code:

coef = [1]
for n in range(1, 20):
coef.append(int(2 * (2 * n + 1) / (n + 2) * coef[n - 1]))

def sample_grammar(n):
if n == 0:
return ''
elif np.random.uniform() < 2 * coef[n - 1] / coef[n]:
return np.random.choice(['a', 'b']) + sample_grammar(n - 1)
else:
u = np.random.uniform()
k = 0
s = coef[0] * coef[n - 2] / (coef[n] - 2 * coef[n - 1])
while s < u:
k += 1
s += coef[k] * coef[n - 2 - k] / (coef[n] - 2 * coef[n - 1])
return 'c' + sample_grammar(k) + 'd' + sample_grammar(n - 2 - k)

counter = Counter([sample_grammar(4) for i in range(100000)])
# Print the counts, normalized so they should all be approximately 1
print('\n'.join(
'{}\t{:.4f}'.format(string, count / sum(counter.values()) * len(counter.values()))
for string, count in sorted(counter.items()))
)

This outputs, for instance,

aaaa    0.9949
aaab    1.0006
aaba    0.9936
aabb    0.9926
aacd    0.9942
abaa    1.0091
abab    0.9979
abba    1.0022
abbb    1.0055
abcd    1.0007
acbd    1.0120
acda    0.9951
acdb    0.9947
baaa    0.9924
baab    0.9992
baba    0.9991
babb    1.0020
bacd    0.9967
bbaa    0.9944
bbab    1.0046
bbba    0.9897
bbbb    0.9996
bbcd    1.0148
bcbd    0.9946
bcda    1.0080
bcdb    1.0028
cabd    0.9960