Probability distribution for different probabilities If I wanted to get the probability of 9 successes in 16 trials with each trial having a probability of 0.6 I could use a binomial distribution. What could I use if each of the 16 trials has a different probability of success?
 A: The (in general intractable) pmf is
$$
  \Pr(S=k) = \sum_{\substack{A\subset\{1,\dots,n\}\\ |A|=k}} \left( \prod_{i\in A} p_i \right)\left(\prod_{j\in \{1,\dots,n\}\setminus A} (1-p_j) \right) \, .
$$
R code:
p <- seq(1, 16) / 17
cat(p, "\n")
n <- length(p)
k <- 9
S <- seq(1, n)
A <- combn(S, k)
pr <- 0
for (i in 1:choose(n, k)) {
    pr <- pr + exp(sum(log(p[A[,i]])) + sum(log(1 - p[setdiff(S, A[,i])])))
}
cat("Pr(S = ", k, ") = ", pr, "\n", sep = "")

For the $p_i$'s used in wolfies answer, we have:
Pr(S = 9) = 0.1982677

When $n$ grows, use a convolution.
A: This is the sum of 16 (presumably independent) Binomial trials.  The assumption of independence allows us to multiply probabilities.  Whence, after two trials with probabilities $p_1$ and $p_2$ of success the chance of success on both trials is $p_1 p_2$, the chance of no successes is $(1-p_1)(1-p_2)$, and the chance of one success is $p_1(1-p_2) + (1-p_1)p_2$.  That last expression owes its validity to the fact that the two ways of getting exactly one success are mutually exclusive: at most one of them can actually happen.  That means their probabilities add.
By means of these two rules--independent probabilities multiply and mutually exclusive ones add--you can work out the answers for, say, 16 trials with probabilities $p_1, \ldots, p_{16}$.  To do so, you need to account for all the ways of obtaining each given number of successes (such as 9).  There are $\binom{16}{9} = 11440$ ways to achieve 9 successes.  One of them, for example, occurs when trials 1, 2, 4, 5, 6, 11, 12, 14, and 15 are successes and the others are failures.  The successes had probabilities $p_1, p_2, p_4, p_5, p_6, p_{11}, p_{12}, p_{14},$ and $p_{15}$ and the failures had probabilities $1-p_3, 1-p_7, \ldots, 1-p_{13}, 1-p_{16}$.  Multiplying these 16 numbers gives the chance of this particular sequence of outcomes.  Summing this number along with the 11,439 remaining such numbers gives the answer.
Of course you would use a computer.
With many more than 16 trials, there is a need to approximate the distribution.  Provided none of the probabilities $p_i$ and $1-p_i$ get too small, a Normal approximation tends to work well.  With this method you note that the expectation of the sum of $n$ trials is $\mu = p_1 + p_2 + \cdots + p_n$ and (because the trials are independent) the variance is $\sigma^2 = p_1(1-p_1) + p_2(1-p_2) + \cdots + p_n(1-p_n)$.  You then pretend the distribution of sums is Normal with mean $\mu$ and standard deviation $\sigma$.  The answers tend to be good for computing probabilities corresponding to a proportion of successes that differs from $\mu$ by no more than a few multiples of $\sigma$.  As $n$ grows large this approximation gets ever more accurate and works for even larger multiples of $\sigma$ away from $\mu$.
A: One alternative to @whuber's normal approximation is to use "mixing" probabilities, or a hierarchical model.  This would apply when the $p_i$ are similar in some way, and you can model this by a probability distribution $p_i\sim Dist(\theta)$ with a density function of $g(p|\theta)$ indexed by some parameter $\theta$.  you get a integral equation:
$$Pr(s=9|n=16,\theta)={16 \choose 9}\int_{0}^{1} p^{9}(1-p)^{7}g(p|\theta)dp $$
The binomial probability comes from setting $g(p|\theta)=\delta(p-\theta)$, the normal approximation comes from (I think) setting $g(p|\theta)=g(p|\mu,\sigma)=\frac{1}{\sigma}\phi\left(\frac{p-\mu}{\sigma}\right)$ (with $\mu$ and $\sigma$ as defined in @whuber's answer) and then noting the "tails" of this PDF fall off sharply around the peak.
You could also use a beta distribution, which would lead to a simple analytic form, and which need not suffer from the "small p" problem that the normal approximation does - as beta is quite flexible.  Using a $beta(\alpha,\beta)$ distribution with $\alpha,\beta$ set by the solutions to the following equations (this is the "mimimum KL divergence" estimates):
$$\psi(\alpha)-\psi(\alpha+\beta)=\frac{1}{n}\sum_{i=1}^{n}log[p_{i}]$$
$$\psi(\beta)-\psi(\alpha+\beta)=\frac{1}{n}\sum_{i=1}^{n}log[1-p_{i}]$$
Where $\psi(.)$ is the digamma function - closely related to harmonic series.
We get the "beta-binomial" compound distribution:
$${16 \choose 9}\frac{1}{B(\alpha,\beta)}\int_{0}^{1} p^{9+\alpha-1}(1-p)^{7+\beta-1}dp ={16 \choose 9}\frac{B(\alpha+9,\beta+7)}{B(\alpha,\beta)}$$
This distribution converges towards a normal distribution in the case that @whuber points out - but should give reasonable answers for small $n$ and skewed $p_i$ - but not for multimodal $p_i$, as beta distribution only has one peak.  But you can easily fix this, by simply using $M$ beta distributions for the $M$ modes.  You break up the integral from $0<p<1$ into $M$ pieces so that each piece has a unique mode (and enough data to estimate parameters), and fit a beta distribution within each piece. then add up the results, noting that making the change of variables $p=\frac{x-L}{U-L}$ for $L<x<U$ the beta integral transforms to:
$$B(\alpha,\beta)=\int_{L}^{U}\frac{(x-L)^{\alpha-1}(U-x)^{\beta-1}}{(U-L)^{\alpha+\beta-1}}dx$$
A: Let $X_i$ ~ $Bernoulli(p_i)$ with probability generating function (pgf):
$$\text{pgf} = E[t^{X_i}] = 1 - p_i (1-t)$$
Let $S = \sum_{i=1}^n X_i$ denote the sum of $n$ such independent random variables. Then, the pgf for the sum $S$ of $n=16$ such variables is:
$$\begin{align*}\displaystyle \text{pgfS} &=  E[t^S] 
\\&=  E[t^{X_1}]  E[t^{X_2}]  \dots E[t^{X_{16}}] \text{          (... by independence)}
\\ &= \prod _{i=1}^{16} \left(1-p_i(1-t) \right)\end{align*}$$
We seek $P(S=9)$, which is:
$$\frac{1}{9!}\frac{d^9 \text{pgfS}}{dt^9}|_{t=0}$$
ALL DONE. This produces the exact symbolic solution as a function of the $p_i$. The answer is rather long to print on screen, but it is entirely tractable, and takes less than $\frac{1}{100}$th of a second to evaluate using Mathematica on my computer.
Examples
If $p_i = \frac{i}{17}, i= 1 \text{ to } 16$, then:  $P(S=9) = \frac{9647941854334808184}{48661191875666868481} = 0.198268 \dots$
If $p_i = \frac{\sqrt{i}}{17}, i= 1 \text{ to } 16$, then:  $P(S=9) = 0.000228613 \dots$
More than 16 trials?
With more than 16 trials, there is no need to approximate the distribution. The above exact method works just as easily for examples with say $n = 50$ or $n = 100$. For instance, when $n = 50$, it takes less than $\frac{1}{10}$th of second to evaluate the entire pmf (i.e. at every value $s = 0, 1, \dots, 50$) using the code below. 
Mathematica code
Given a vector of $p_i$ values, say:
n = 16;   pvals = Table[Subscript[p, i] -> i/(n+1), {i, n}];

... here is some Mathematica code to do everything required:
pgfS = Expand[ Product[1-(1-t)Subscript[p,i], {i, n}] /. pvals];
D[pgfS, {t, 9}]/9! /. t -> 0  // N


0.198268

To derive the entire pmf:
Table[D[pgfS, {t,s}]/s! /. t -> 0 // N, {s, 0, n}]

... or use the even neater and faster (thanks to a suggestion from Ray Koopman below):
CoefficientList[pgfS, t] // N

For an example with $n = 1000$, it takes just 1 second to calculate pgfS, and then 0.002 seconds to derive the entire pmf using CoefficientList, so it is extremely efficient.
