Can Bernoulli random variables be used to approximate more than just the normal distribution? Most statistics students are familiar with the normal approximation of the binomial distribution. And since binomial distributions are created from sums of Bernoulli random variables, it would follow that with a linear combination of enough Bernoulli random variables you could approximate any normal distribution.
So I'm curious: can you use combinations of Bernoulli random variables to make arbitrarily precise approximations of more/most common families of random variables? For example, if you used Bernoulli distributions to randomly turn on and off a string of bits that was the seed for a linear congruential generator, couldn't you then use inverse transform sampling to create all sorts of distributions? 
Alternatively, if you could map a normal distribution to another distribution, wouldn't you by extension have a way of creating that distribution from Bernoulli random variables?
 A: Bernoulli random variables can approximate (almost) any distribution to arbitrary accuracy: A sequence of Bernoulli values gives us a binary sequence, which can be interpreted as the binary representation of a real number.  (This is not surprising, given that a real number is essentially just an infinite sequence of discrete digits.)  Hence, by an appropriate mapping, we can transform a Bernoulli sequence into a continuous uniform random variable.  Once we have done this, we can use the standard technique of inverse transformation sampling to get a random variable from an arbitrary distribution.
Now, there is some limitation on this, because in practice we never have an infinite sequence of Bernoulli values, but we can generate a large finite sequence.  This allows us to approximate a uniform random variable arbitrarily well, and so we can then approximate any distribution that can be approximated arbitrarily well by a mapping of a random variable that is arbitrarily close to a uniform random variable.

Mathematical details for generation with an infinite sequence: Suppose you want to generate a scalar random variable with distribution function $F$.  To do this, we consider an exchangeable binary sequence $X_1, X_2, X_3, ... \sim \text{IID Bern}(\tfrac{1}{2})$ and define the corresponding random variables:
$$A = A(\boldsymbol{X}) = \inf \Big\{ r \in \mathbb{R} \Big| F(r) \geqslant U \Big\} \quad \quad \quad U = U(\boldsymbol{X}) = \sum_{i=1}^\infty \frac{X_i}{2^i} \sim \text{U}(0,1).$$
(This function is well-defined, by the completeness of the real numbers.)  Now, since $F$ is a non-decreasing function, we have:
$$\mathbb{P}(A \leqslant a) = \mathbb{P} \Big( \inf \Big\{ r \in \mathbb{R} \Big| F(r) \geqslant U \Big\} \leqslant a \Big) = \mathbb{P} \Big( U \leqslant F(a) \Big) = F(a).$$
(Note that this result does not require continuity of $F$, so it works for general distributions, not just continuous distributions.)

Mathematical details for generation with a finite sequence: The above case is an idealised case where we can generate an infinite sequence of Bernoulli random variables.  We now consider the more realistic case where we can generate some arbitrarily large finite sequence with $k \in \mathbb{N}$ terms.  Hence, we now have the finite sequence $X_1, X_2, ..., X_k \sim \text{IID Bern}(\tfrac{1}{2})$ and we define the corresponding random variables:
$$A_k = \inf \Big\{ r \in \mathbb{R} \Big| F(r) \geqslant U_k \Big\} \quad \quad \quad U_k = \sum_{i=1}^k \frac{X_i}{2^i} + \frac{1}{2^{k+1}}.$$
(We have included an additional "continuity correction" term in $U_k$ so that its distribution is still symmetric around the value $\mathbb{E}(U_k) = \tfrac{1}{2}$.)  We now have:
$$\mathbb{P}(A_k \leqslant a) = \mathbb{P} \Big( \inf \Big\{ r \in \mathbb{R} \Big| F(r) \geqslant U_k \Big\} \leqslant a \Big) = \mathbb{P} \Big( U_k \leqslant F(a) \Big).$$
For large $k$, we then have:
$$\mathbb{P}(A_k \leqslant a) = \mathbb{P} \Big( U_k \leqslant F(a) \Big) \approx \mathbb{P} \Big( U \leqslant F(a) \Big) = F(a).$$
As you can see, this approximation relies on our ability to approximate the event $U \leqslant F(a)$ by the event $U_k \leqslant F(a)$, for large $k$.  For all but pathological distributions this approximation can be made arbitrarily close by taking $k$ to be sufficiently large.  There are some pathological distributions where this is not the case (e.g., any distribution with non-zero probability on the irrational numbers, or more broadly, on real values that cannot be represented as a finite binary number), but this is quite a narrow class of distributions.  Hence, this technique will approximate a random variable with (almost) any distribution to an arbitrary degree of accuracy.
A: Yes, depending on circumstances. Binomial and negative binomial (including geometric) distributions are defined in terms of Bernoulli random variables so the relationship of Bernoulli random variables to those distributions is not an approximation. Especially for large $n$ and small $p$,
$\mathsf{Binom}(n, p)$ is approximately $\mathsf{Pois}(\lambda = np).$
Your purpose is not clear, but once you have approximated normal distributions in terms of Bernoulli random variables, you could exploit the relationship of normal distributions to chi-squared, t (extended to Cauchy), and F distributions to approximate those by simulation as well.
In the same spirit, you might exploit the relationship of Poisson to exponential random variables, and go further to approximate Laplace, gamma (at least with integer shape parameters), and some Weibull distributions. [This is a partial list.]
Simulations generally begin from pseudorandom realizations of $\mathsf{Unif}(0,1),$ and other distributions can then be simulated in terms of uniform random variables. Simulations in terms of Bernoulli random variables might be somewhat more
awkward and limited than simulation in terms of uniform random variables, but there are, nevertheless, many possibilities. M/M Queueing processes (defined in continuous time) can be simulated approximately by 'discretizing time' and using Bernoulli instead of exponential random variables in the simulation.
A: This will work in most cases, at least if the Bernoulli variables involved are independent and identically distributed, or at least exchangeable and identically distributed (Peres 1992), and an unlimited sequence of them is available. This answer will include discussion on algorithms to approximate a distribution using Bernoulli random variables.
The rest of this answer assumes the Bernoulli random variables have a mean (bias) of 1/2 (they represent unbiased random bits). If the random variables (bits) have an unknown bias, a randomness extraction method such as the von Neumann or Peres extractor can be used to generate unbiased random bits from them (see my note on randomness extraction).
There are many ways to approximate a distribution to an arbitrary accuracy, and they are discussed in Devroye and Gravel 2015/2020:

*

*For discrete distributions, there are many possibilities. One is to take the binary expansions of the distribution's probabilities and perform a random walk on them, driven by unbiased random bits. This is the essence of the DDG tree algorithm by Knuth and Yao. In fact any practical algorithm for discrete distributions can be recast to a DDG tree algorithm. See also Gryszka 2020.
The binary entropy of a discrete distribution is a lower bound on the average number of random bits needed to produce a variate of that distribution. Unfortunately, some discrete distributions have unbounded entropy; Devroye and Gravel give as an example the zeta Dirichlet distribution (in certain parameterizations).


*For continuous distributions, Devroye and Gravel discuss inversion, discretization, convolution, and rejection sampling to produce a variate from a given distribution using random bits, and give lower bounds on the number of random bits needed to generate a variate with the desired accuracy. Different approaches work well for different distributions. One notable example of generating a continuous variate using unbiased bits is the exponential distribution, such as C.F.F. Karney's improvement in "Sampling exactly from a normal distribution" of von Neumann's time-honored exponential generator. Also, I specify algorithms to do arithmetic on random variates represented by sequences of bits or digits. On the inversion side, if a distribution's inverse CDF is known, enough bits of a uniform(0,1) random variate can be generated so that the inverse CDF can be calculated from that variate with a low enough variance.
REFERENCES:

*

*Devroye, L., Gravel, C., "Random variate generation using only finitely many unbiased, independently and identically distributed random bits", arXiv:1502.02539v6 [cs.IT], 2020.

*Peres, Y., "Iterating von Neumann's procedure for extracting random bits", Annals of Statistics 1992,20,1, p. 590-597.

*Knuth, Donald E. And Andrew Chi-Chih Yao. "The complexity of nonuniform random number generation", in Algorithms and Complexity: New Directions and Recent Results, 1976.

*Gryszka, K. From biased coin to any discrete distribution. Period Math Hung (2020). https://doi.org/10.1007/s10998-020-00363-w
