Bad coin factory: Do coins in a set tend to be unfair? The question I'm trying to answer is (fundamentally) this: 
I have a bag of coins that I suspect are weighted, some towards heads, some towards tails. I toss each coin 4 times and record the outcomes (e.g., 3H1T). As a group, do the coins tend to be unfair? 
I can't figure out what an appropriate test would be, though it seems like there ought to be one. Here are some relevant thoughts and options I've considered.
(1) Binomial test - Appropriate way to test EACH COIN's fairness, but (a) 4 tosses isn't enough for statistical significance ($\alpha$ = .05) at the level of the individual coin and (b) since I suspect different coins may be weighted in opposite directions, lumping all the data together would make these coins cancel each other out. (See similar comments on this question)
(2) Chi-square goodness-of-fit or multinomial test over counts - This will tell me if my observed counts for each outcome (4H0T, 3H1T, 2H2T...) differ from the expected counts (they do), but not how. It will return a high test statistic whether my coins are all magically fair (all 2H2T results) or if they are all weighted (all either 4H0T, or 0H4T). It also ignores the underlying binomial nature of the data.
(3) Regression seems like overkill for this data, and linear regression/linear mixed models wouldn't answer the right question anyway: coins with opposite weighting would cancel each other out.
For reference, my actual counts are as follows, of a total of 56 "coins". Unfortunately, since they're not real coins, and the experiment is over, I can't just go flip each one a few more times!
16 4H, 10 3H, 9 2H, 8 1H, 13 0H
 A: The Multinomial Test seems appropriate. You can do this through a significance test, or by using Bayesian inference with a Dirichlet prior. For a frequentist significance test, your null hypothesis is:
$$
\begin{split}
P(4H) &= 1/16 \\
P(3H) &= 1/4 \\
P(2H) &= 3/8 \\
P(1H) &= 1/4 \\
P(0H) &= 1/16.
\end{split}
$$
That's the probability distribution of an outcome of four flips of a fair coin. What you're trying to find is if the distribution of "head count" outcomes represents the distribution given by a set of fair coins. While this won't answer the question of whether any one coin is unfair (which, for the reasons you stated in your question, you just cannot do) this does help answer whether your set of coins are all fair. If the probability of getting an outcome that's as likely or less likely than yours given that they're all fair coins, then that's a reason to reject the null hypothesis.
As mentioned in the article I linked, the Multinomial test is intractable for a large number of samples and categories. That's because, for $N$ samples and $k$ multinomial categories, there are $N + k - 1 \choose k - 1$ possible outcomes, and you have to figure out which outcomes, from that set, are less likely than the one you got. So an approximation scheme like the likelihood ratio test might be appropriate. In your case ($N=56$, $k=5$) there are 490 thousand possible outcomes. If you don't mind the computational load, you can brute force calculate the probability of all those outcomes and search for which ones are lower than, or equal to, the probability of your outcome. The sum of those probabilities is your $p$-value. Otherwise you can use the LR test as illustrated in the wiki article, which is simple enough.
A: What is unknown is about the factory. You want to test a property of the factory (not the coins). 
While a coin can be entirely described by a single number (the probability of "heads"), a factory is naturally described by a distribution : the distribution of the random variable $P$ that is the probability of heads for a coin produced by this factory. The parameter $\theta$ is thus a distribution (not necessarily continuous) over $[0;1]$.
Now you can't observe $P$ directly : $P$ is latent. You can only observe a variable $X$ whose distribution given $P$ is $B(4,P)$. First let's calculate the distribution of $X$ (given $\theta$) :
$$P_\theta(X=k)=E_\theta\left(\binom{4}{k}P^k(1-P)^{4-k}\right)=\int_0^1\binom{4}{k}p^k(1-p)^{4-k}d\theta(p)$$ for $k=0,1,2,3,4$
This distribution $d(\theta)$ is an element $[0;1]^5$. $d(\theta)$ is not always a binomial distribution. It is a binomial distribution when $\theta$ is a Dirac. $(0.5;0;0;0;0;0.5)$ is a possible value for $d(\theta)$ that is not binomial. Yet $X$ cannot have ANY possible distribution. For example, the distribution $(0,0,1,0,0)$ is not a possible value for $d$.
Let's look at the set of all possible distributions :
$$D=\{d(\theta)/\theta \text{ is a probability distribution over [0;1]}\}$$
Finding $D$ is maybe difficult. I won't try to do it. But you will have to I guess. It is convex : image of the convex space of distributions by a linear map. I think it must be a kind of polyhedron (like a simplex).
Now the idea is to use $d\in D$ as the parameter instead of $\theta$. The null hypothesis is $d=(1/16;1/4;3/8;1/4;1/16)$. 
You can calculate the likelihood of a vector of 56 independent variables having $d$ as a distribution (multinomial stuff). Then you can use a likelihood ratio test as suggested by Bridgeburner. But the range of the parameters won't be every possible distribution, only those in $D$. See : https://en.wikipedia.org/wiki/Likelihood-ratio_test
I think this solves the problem of rejecting if all values are 2 : the distribution $(0;0;1;0;0)$ is not an element of $D$.
The difficult point seems to find $D$ and to calculate the maximum of the likelihood over $D$.
