# Illustrate Bayesian inference. GOAL: Illustrate Bayesian inference. Probability is the codification of logical thought in the presence of uncertainty. It the ONLY consistent system for logical inference with uncertainty (TODO: Cox's Theorem, Fisher's Information, Cramer-Rao bound, KL co-ordinate free invariant information geometry) When interpreting the world you get observations $O$ and you want to decide what model of truth $T$ best explains those observations. This inference of model from data sometimes allows you to understand and predict things about the world (The Scientific Method). To do this inference, you do experiments to determine what you observe given you have fixed and know the truth, $P(O|T)$, the likelihood / density (interpreting as function of T / O). Bayes rule then turns that around to compute $P(T|O)$ the posterior (after-data) of the true hypothesis given what you observe and your prior about what might be true before seeing any data: $P(T)$: $$P(T|O) = P(O|T)*P(T) / P(O)$$ This is inspired from MacKay's book in which he shows the simple intuitive beauty of Bayes rule: <a href="http://www.inference.org.uk/mackay/itprnn/ps/47.59.pdf">http://www.inference.org.uk/mackay/itprnn/ps/47.59.pdf</a> For example, say I want to predict the species for a given tree based on its height. <img src="redmaple2__58290.1538684018__32799.1539281421.jpg"> <img width="300px" src="Kousa_Dogwood_1_540x.jpg"> First I construct some probability distributions. I go out and look at only maple trees and get their heights. I then I look at only dogwood trees and get their heights. I capture my observations in functions $P(height|maple), P(height|dogwood)$. Maples are generally taller (70 feet) than dogwoods (30 feet). Now I have a tree that is 69 feet; what species is it? I want to compute the posterior probability of it being each species of tree: $P(maple | 69 feet), P(dogwood | 69 feet)$. I might also have a prior that I might see more maples because I'm farther north. I look at which posterior has greater probability and decide that species of tree. Skipping all the details, this procedure would likely call Maple because 69 feet is much closer to the average 70 of Maple than the average of 30 feet of Dogwood. The details make this all very precise and quantitative. Probability distributions forces you to place rational bets across all possibilities. This allows you precisely weigh how often different events happen. Bayesian probability mean you write down the probability of everything. After that it is _mechanical_ to infer different quantities. ================================ # Simple Binomial Example Let's work through a simple example: <a href="https://en.wikipedia.org/wiki/Binomial_distribution">https://en.wikipedia.org/wiki/Binomial_distribution</a> Flip a coin $N$ times where the probability of heads is $p$. How often do you see $k$ heads where $k$ is between $[0,1,...,N]$? The binomial distribution tells you this: $$ P(k|N,p) = choose(N, k) p^k (1-p)^{(N-k)} $$ The three terms are 1) probability of heads to the power of the number of observed heads $(p^k)$, 2) probability of tails to the power of number of observed tails $(1-p)^{(N-k)}$, 3) how many ways are there are to choose number of observed heads out of the total $(choose(N,k) = N!/(k! * (N-k)!))$ or the binomial coefficient. These distributions are carefully constructed to sum to 1.0. The powers of probability are simply because probabilities multiply for independent events. This is mathematical truth. If you independently flip coins, binomial describes exactly the outcomes. ================================ # Generate Some Data Let's generate some data from the model using R. Binomial data is generated using "rbinom". Set the model parameters to $N=10$ flips, $p=0.3$ or 30% for heads. Let's run 100 experimental trials. ``` N=10 p=0.3 trials= 100 obs = rbinom(trials, size=N, prob=p) png("tableobs.png") plot(table(observed)) title("counts of observed values") dev.off() ``` <img src="tableobs.png"> The count of about 30 at observed=3 says that out of the 100 experiments you ran, ~30 of those experiments showed 3 heads and 7 tails when you flipped the coin ten times. ================================ # Inference OK, we have some data. Let's forget that we know the probability of heads that generated the data (0.3). Let's infer is value from the data itself where we know we flipped a coin 10 times for 100 trials. To make things simple, let's assume the true probability is quantized so that it can only take on 1/10'th values. There are nine different hypotheses to test: ``` hyp = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9) ``` We want to know the posterior probability of each hypothesis given the data. This says after you are given the data, what is the probability of each of the hypotheses. You can use this and do things like take the maximum one as the best bet as to which one is the correct one. How do we do that? Write down Bayes rule: $$P(T|O) = P(O|T)*P(T) / P(O)$$ Here $T$ or truth is over the possible hypotheses and $O$ is the observed data. Let's fill out the terms. We assume a model. Here it is binomial but it could be any model that assigns a probability to a data point given some parameters of the model. (Note: here the model "reflects" reality. What happens when the model or class of models is limited and doesn't exactly reflect reality? In machine learning, is your hypothesis concept class capable of representing the problem? ). R binomial probability is "dbinom". $$P(O|T) = P(k|N,p) = choose(N,k) p^k (1-p)^{(N-k)} = dbinom(k,size=N,prob=p)$$ The observed data is $k$ the number of heads and the truth is the setting of the model parameters $N,p$. When you sum over all possibilities of the observed data (0:N) you get 1.0, or that you've placed rational bets on all possibilities. $P(T)$ is the prior. This is your a-priori belief that each of the hypotheses is true. In our case we don't have any real information and therefore assign uniform probability over the nine hypotheses. $P(O)$ is a normalizing factor to make sure the posterior sums to 1.0 over the different hypotheses. It is $\sum_i P(O|T_i)*P(T_i)$. We have all the ingredients! ================================ # Inference Calculation: Likelihoods Compute some numbers! Both the likelihood and the posterior are 2d arrays $[O,T]$. In this case, $O$ can take on 11 possible values (0:10 the number of heads out of 10 flips) and $T$ can take on 9 values (0.1:0.9, the 9 hypothesis values that we decided to consider). Let's use R to compute these things. ``` # storage [hpy,obj] likelihood = matrix(nrow=9, ncol=11) # step through all hyptheses for (hypii in 1:length(hyp)){ # step through all possibilties for running one trial, # the number of observed heads obspossible = seq(0,N) for (obsii in 1:length(obspossible)){ # compute the likelihood for this observed, hypothesis combination likelihood[hypii,obsii] = dbinom(obspossible[obsii], size=N, prob=hyp[hypii]) } } write.table( likelihood, file="likelihood.tsv", sep="\t", col.names=F,row.names=F) write.table( hyp, file="likelihood_rowval.tsv", sep="\t", col.names=F,row.names=F) write.table( obspossible, file="likelihood_colval.tsv", sep="\t", col.names=F,row.names=F) ``` OK, now we have the likelihood for each possible hypothesis over every possible observed outcome using a simple mechanical computation. Let's make sure that our distributions properly sum to 1.0. The sum across the columns should be 1.0 for each of the 9 hypotheses: ``` apply(likelihood, 1, sum) [1] 1 1 1 1 1 1 1 1 1 ``` GOOD! Each hypothesis places a distribution over all observed possibilities that sum to 1.0 Let's look at all the hypotheses: ``` for (hh in 1:length(hyp)){ fn = paste("hypoth.",hh,".png",sep="") png(fn) plot(obspossible, likelihood[hh,],type="b") title(paste("hypothesis p=",hyp[hh])) dev.off() } mycol=rainbow(9) png("hypoth.all.png") plot(obspossible, likelihood[1,],col=mycol[1],type="b",ylim=c(0,0.4)) for (hh in 1:length(hyp)){ points(obspossible, likelihood[hh,],col=mycol[hh],type="b",ylim=c(0,1)) } dev.off() ``` These are plots of the likelihood over the observations 0:10 given the hypothesis probability. <img src="hypoth.all.png"> SlideIt: <input id="val1" type="range" min="1" max="9" value="1" step="1" oninput="showVal1(this.value)" onchange="showVal1(this.value)" /> <span id="range">1</span> <img id="img1" src="hypoth.1.png"> A plot shows that given a particular hypothesis probability of heads say 0.1 and flip a coin 10 times, the probability distribution over the outcomes of a total of 0 heads, 1 heads, ... 10 heads given my probabilistic model. For binomial and p=0.1 it says there is about 40% probability of seeing 1 heads out of the 10 flips. ================================ # Likelihood Dual / Posterior Let's look at the dual or along the other axis of the matrix we constructed. This is the value of the likelihood at each hypothesis=(0.1,...0.9) given a single observation of 0:11. First note that they do not sum to 1.0: ``` apply(likelihood, 2, sum) [1] 0.4914342 0.8287079 0.9061899 0.9099927 0.9091388 0.9090730 0.9091388 [8] 0.9099927 0.9061899 0.8287079 0.4914342 ``` ``` for (oo in 1:length(obspossible)){ fn = paste("obs.",oo,".png",sep="") png(fn) plot(hyp, likelihood[,oo],type="b") title(paste("obs=",obspossible[oo])) dev.off() } mycol=rainbow(11) png("obs.all.png") plot(hyp, likelihood[,1],col=mycol[1],type="b",ylim=c(0,0.4)) for (oo in 1:length(obspossible)){ points(hyp, likelihood[,oo],col=mycol[oo],type="b",ylim=c(0,0.4)) } dev.off() ``` These are the un-normalized distributions of the different hypotheses given a single observation of obs: <img src="obs.all.png"> SlideIt: <input id="val2" type="range" min="1" max="11" value="1" step="1" oninput="showVal2(this.value)" onchange="showVal2(this.value)" /> <span id="range2">1</span> <img id="img2" src="obs.1.png"> Note that these are very similar to the likelihood plots from before. There is a dual symmetry. The dual conjugate to the <a href="https://en.wikipedia.org/wiki/Binomial_distribution">binomial distribution</a> is the beta distribution: <a href="https://en.wikipedia.org/wiki/Beta_distribution">beta distribution</a>. The distributions are EXACTLY the same except for the normalizing factor and interpretation! (Which might be hard to tell from the equations at first glance) A plot shows that given a particular observation from my experiment, say 0 heads out of 10 flips, show me the probabilities over the different hypotheses (0.1, 0.2, ..., 0.9). This is the dual of the likelihood plots shown first. We are taking slices through the probability distributions: <img style="box-shadow: 0 4px 8px 0 rgba(0, 0, 0, 0.2), 0 6px 20px 0 rgba(0, 0, 0, 0.19);" width=400px src="bayesSlice.svg.png"> Note these "dual" distributions don't sum to 1.0! They are still values between 0 and 1 and can be compared against each other because we've taken effort to make sure our models return proper probability distributions (viz: examine the maximum value). However, we can't interpret them as probability _distributions_ yet. ***This is very important!:*** How can we interpret these numbers as probability distributions? Bayes rule says to simply multiply by prior and normalize. - Going one way, we can propose a multiple hypothesis models to make predictions about how the outcome of an experiment will turn out (the likelihood) - We can also quantify how much we believe each model is likely to be true (the prior). - We can then run experiments and let Nature tell us the results (the observations). - We can then mechanically plug these resulting observations into our models and make _inferences_ about what how likely each of the different models is to be ***true*** (the posterior). - This ***is*** the scientific method and codifies logical thought under uncertainty. Note that for uniform prior, this multiplying by prior and normalizing is the same as dividing by the sum, as the prior cancels out, and does not change the shape. For non-uniform prior, more probability will be put on those with higher a-prior probability. We do this in the next section. Rather than cuts parallel to the two axes, you can also plot the likelihood as a 2d image: ``` png("like.2d.png") image(likelihood) dev.off() ``` <img src="like.2d.png"> Or in 3d: <a href="three.html">three.html</a> <iframe src="three.html" width=800 height=800></iframe> ================================ # Uniform Prior Say we see a single of observation of 3 out of 10. First let's pull out the probability of observing a single 3 (out of 10) under the different hypotheses (obspossible[4]= 3 because of the 0 obs) ``` Pobs3 = likelihood[,4] [1] 0.057395628 0.201326592 0.266827932 0.214990848 0.117187500 0.042467328 [7] 0.009001692 0.000786432 0.000008748 sum(Pobs3) [1] 0.9099927 which.max(Pobs3) 3 png("Pobs3.png") plot(hyp, Pobs3,type="b") dev.off() ``` <img src="Pobs3.png"> Good! If you observe a single observation of 3 (out of 10) then the hypothesis with the highest likelihood is p=0.3 with probability 26.7% . This is the maximum likelihood choice. The next most likely is p=0.4 at 21.5%. Here we can compare probabilities through things like maximum even though they do not sum to one and are not a distribution Let's look at the true posterior with uniform prior properly normalized: ``` sum(Pobs3) [1] 0.9099927 # doesn't sum to 1 because over hypotheses # uniform prior and normalize prior = rep(1/9,9) posterior = prior*Pobs3 posterior = posterior/sum(posterior) posterior [1] 6.307262e-02 2.212398e-01 2.932199e-01 2.362556e-01 1.287785e-01 [6] 4.666777e-02 9.892049e-03 8.642179e-04 9.613264e-06 png("Pobs3post.png") plot(hyp, posterior,type="b") dev.off() ``` <img src="Pobs3post.png"> As you can see, the uniform prior posterior has changed the scale, but the shape remains unchanged. But now summing over all hypotheses gives 1.0 so it is a nice distribution. Note that getting the "partition function" normalizations are correct so it sums to 1.0 is key if you want probabilities. With more complicated models, this can sometimes be tricky. So much so that fields like Bayes decision theory construct functions that look at ratios so the these normalizations cancel out. With a real probability distribution, you can also also talk about the maximum a posterior choice (MAP), which for uniform prior is the maximum-likelihood choice. There is the mean posterior estimate as alternative to the maximum likelihood estimate: ``` sum(hyp*posterior) [1] 0.3330378 ``` The mean posterior estimate is 0.333 which is close the maximum likelihood estimate. Note it doesn't lie right on a hypothesis 0.3, 0.4. This is the Bayesian way. You keep around the distributions of everything carrying full information and then possibly take reductions like max or mean if you want a single estimate. ================================ # Informative Prior Let's say we believe a-priori that the coin is biased and will come up heads 30% of the time. Let's put a prior probability 76% on the 0.3 hypothesis and uniform over others. ``` hypPrior= c( 0.03, 0.03, 0.76, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03 ) posterior = likelihood # set up storage for (hypii in 1:length(hyp)){ for (obsii in 1:length(obspossible)){ posterior[hypii,obsii] = likelihood[hypii,obsii]*hypPrior[hypii] } } for (obsii in 1:length(obspossible)){ # normalize posterior[,obsii] = posterior[,obsii]/sum(posterior[,obsii]) } for (oo in 1:length(obspossible)){ fn = paste("post3.",oo,".png",sep="") png(fn) plot(hyp, posterior[,oo],type="b") title(paste("prior=3 76% obs=",obspossible[oo])) dev.off() } apply(posterior,2,sum) [1] 1 1 1 1 1 1 1 1 1 1 1 ``` SlideIt: <input id="val3" type="range" min="1" max="11" value="1" step="1" oninput="showVal3(this.value)" onchange="showVal3(this.value)" /> <span id="range3">1</span> <img id="img3" src="post3.1.png"> The sums over hypotheses are to 1. You can see how the prior really favors the 0.3 hypothesis. Only when you seen an observation of 7 or greater does the likelihood overcome the prior and change the most likely call. ================================ # Multiple Observations Of Same Outcome Before we had the probability given a single observation: $P(obs | hyp)$ or likelihood[hyp,obs] For multiple observations if you assume independence then you can multiply together: $$P({obs} | hyp ) = \prod_N P(obs_N | hyp)$$ ( Multiplying makes probabilities smaller because they are less than or equal to 1.0. With large number of observations, these probabilities can get very small (Think 10^-100 for 100 observations each at 10% probability). Are they still a distribution because they are so small? Yes, because the input space is much larger for the multiple observations as you have to sum over all possibilties of _all_ observations. See appendix. ) Bayes rule applies in the same way to multiple observations. Let's say you see 5 observations of 3. What is the likelihood under the 9 different hypotheses for these multiple observations? This is just probability of each hypothesis given the single observation of 3 (Prob3 = $P(3|H_i)$) to the 5th power because observations are independent. For example, for just one hypothesis $$ P( (3,3,3,3,3) |H) = P(3|H) P(3|H) P(3|H) P(3|H) P(3|H) = (P(3|H))^5 $$ ``` Pobs3^5 [1] 6.228652e-07 3.307545e-04 1.352560e-03 4.593036e-04 2.210072e-05 [6] 1.381258e-07 5.910453e-11 3.008194e-16 5.123230e-26 # With uniform prior, shape is UNCHANGED prior = rep(1/9,9) posterior = prior*Pobs3^5 posterior = posterior/sum(posterior) [1] 2.876338e-04 1.527396e-01 6.246005e-01 2.121025e-01 1.020592e-02 [6] 6.378530e-05 2.729397e-08 1.389158e-13 2.365864e-23 png("Pobs3pow5.png") plot(hyp, posterior,type="b") points(hyp, Pobs3/sum(Pobs3) ,type="b",col="red") title("Black: posterior 5 observations of 3/10.\nRed: posterior of single 3/10 observation") dev.off() ``` <img src="Pobs3pow5.png"> The maximum likelihood choice is again p=0.3 with probability 1.3e-03. But you can see that the likelihood is peaking up on 0.3 and falling off rapidly on other hypotheses because of the multiplication. The posterior is 62.5% at the 0.3 hypothesis With a large number of 3-observations 100, almost all the probability falls on 0.3. Note that with such small probabilities you must be careful of floating point underflow. In this case computing the posterior puts machine 100% probability on 0.3 with "rounding error" at the other hypotheses. As shown later, you usually solve this problem by storing log(probability) where multiplication is addition and addition involves moving in and out of log using exponential and log functions (slower than multiplying probabilities). ``` Pobs3^100 [1] 7.724693e-125 2.455354e-70 4.198837e-58 1.745811e-67 7.729100e-94 [6] 6.389894e-138 2.706543e-205 3.682287e-311 0.000000e+00 prior = rep(1/9,9) posterior100 = prior*Pobs3^100 posterior100 = posterior100/sum(posterior100) [1] 1.839722e-67 5.847701e-13 1.000000e+00 4.157844e-10 1.840772e-36 [6] 1.521825e-80 6.445934e-148 8.769779e-254 0.000000e+00 sum(posterior) [1] 1 png("Pobs3pow5c100.png") plot(hyp, posterior100,type="b",col="blue") points(hyp, Pobs3/sum(Pobs3) ,type="b",col="red") points(hyp, posterior ,type="b") title("Black: posterior 5 observations of 3/10.\nRed: posterior of single 3/10 observation\nBlue: 100 observations") dev.off() ``` <img src="Pobs3pow5c100.png"> In general the likelihood under these multinomial distributions is the probability to the number of observations: $$P( obs | hyp) = \prod_i (hypP_i) ^ {obs_i}$$ ================================ # Multiple observations Of Different Outcomes We've seen if you observe the same value several times then the likelihood peaks up on the more probable hypotheses looking at single observation. But what if you observe multiple observations at different points, where does the likelihood peak up? Say you observed the following observations after running 102 different experiments: ``` myobs = c(3,11,24,27,20,10,4,2,1,0,0) sum(myobs) [1] 102 cbind(obspossible,myobs) obspossible myobs [1,] 0 3 [2,] 1 11 [3,] 2 24 [4,] 3 27 [5,] 4 20 [6,] 5 10 [7,] 6 4 [8,] 7 2 [9,] 8 1 [10,] 9 0 [11,] 10 0 ``` What is the likelihood of this ensemble of data under the different hypotheses? For each hypothesis take a product of the probability estimated by that hypothesis to the number of times it is seen in the data. ``` mylike = sapply(1:9, function(hh){ prod(likelihood[hh,]^myobs)}) [1] 1.447318e-156 1.882205e-97 1.239293e-82 2.775077e-90 2.104950e-115 [6] 1.317852e-158 2.080566e-225 0.000000e+00 0.000000e+00 mylike/sum(mylike) [1] 1.167858e-74 1.518773e-15 1.000000e+00 2.239241e-08 1.698509e-33 [6] 1.063389e-76 1.678832e-143 0.000000e+00 0.000000e+00 -log2(mylike) [1] 517.6874 321.3146 272.0886 297.5010 380.9479 524.4665 746.3768 Inf [9] Inf ``` The maximum likelihood is p=0.3 that is 8 orders of magnitude greater than p=0.4 given the data (The posterior gets machine 100% due to underflow). Because probabilities get so small with larger number of observations, taking -log2 not only makes numerics easier but also can be interpreted as bits. Even here, you run into underflow issues that make normalization tricky. Here p=0.3 has (- 297.5010 272.0886)=24.5 bits in support of it over the p=0.4 hypothesis. Roughly seeing this as a log odds ratio, the probability of making a classification error is $2^{-24.5}=4.2E-08$. Thinking in bits is useful! When there are two hypotheses, examining $LR= log_2 P1/P2$, the log probability ratio, has nice properties in Bayes Decision Theory. You decide P1 is the winner if LR is positive otherwise P2. $LR$ gives the bits in support over the alternative. ================================ # Maximum Likelihood and KL From before we saw the likelihood from multiple independent observations: $$P( obs | hyp) = \prod_i (hypP_i) ^ {obs_i}$$ Say $hypP$ is a set of hypotheses ${hypA, hypB, ...}$. Maximum likelihood is: $$\max_X \prod_i (hypX_i) ^ {obs_i}$$ Take the negative log and now minimizing as log in monotonic and negative flips $$\min_X - \sum_i obs_i \log (hypX_i)$$ Divide by total number of observations $N$, a constant $$\min_X - \sum_i (obs_i/N) \log (hypX_i)$$ Introduce Kullback-Liebler between distributions P and Q $$KL(P|Q) = \sum_i P_i * \log( P_i/Q_i )$$ $$KL(P|Q) = \sum_i P_i \log( P_i) - \sum_i P_i \log( Q_i )$$ Fix $P$ to be the truth and find min KL to $hypX$, the different hypotheses to truth $$min_X KL(P|hypX) = \sum_i P_i \log( P_i) - \sum_i P_i \log( hypX_i )$$ Note the first term the entropy of $P$ is constant $$min_X - \sum_i P_i \log( hypX_i )$$ This is exactly like maximum likelihood except $obs_i/N$ is replaced by $P_i$. As $N$ gets big $(obs_i/N) \rightarrow P_i$ by the strong law of large numbers. Finding the maximum likelihood hypothesis given the observed set of data is the same as computing KL between the hypothesis' distributions and the distribution from the observed set of data. The nice thing about thinking in KL minimization is that you can apply a bound from Hoeffding related to Bayes decision theory: $$ P(\mbox{Err}) = P( \log \frac {p_1(O)}{p_2(O)} > 0 | O \sim p_2 ) < \exp(-2N KL(p_2, p_1)^2) $$ This says the that the probability of making a misclassification error by choosing the maximum posterior through the log posterior ratio, goes down exponentially fast with the number of data points $N$ and the square KL divergence in the competing hypotheses $p1,p2$. So with large data and hypotheses that don't have the same shape you will make very few errors. So taking your observations and estimating a distribution and comparing that distribution through KL to a set of competing distributions is the same as maximum likelihood! KL is looking for the same shape of distribution. Let's compute KL for the data in the previous: ``` kl = function(pp, qq){ sum(pp*log2(pp/qq)) } # estimate distribution from observed data with 0.01 Laplace pseudo count myobsP = (myobs+0.01)/sum(myobs+0.01) # compute KL from observed to hypothesis distributions sapply(1:length(hyp), function(hh){ kl(myobsP, likelihood[hh,])}) [1] 2.43565192 0.50803132 0.02383006 0.27166534 1.08857409 2.49441972 [7] 4.66870578 8.10767875 14.48080813 png("kl.png") plot(obspossible,likelihood[2,],type="b",col="red") points(obspossible,likelihood[3,],type="b",col="green") points(obspossible,likelihood[4,],type="b",col="blue") points(obspossible,myobsP,type="b",lwd=2) title("data = black. hypothesis=0.3=green 0.4=blue 0.2=red\nKL(black,green)=0.0238 ,blue)=0.27 ,red)=0.5") dev.off() ``` <img src="kl.png"> Yes the minimum KL is the maximum likelihood answer! You can see that green and black have almost the exactly same shape. So with multiple observations over several values, the likelihood peaks up on the hypothesis that has the most similar shape as measured by KL! TODO: what if the data has nothing to do with your models and the "fit" is horrible? KL is the number of bits required to communicate the differences between two distributions. This example has 102 data points. Let's construct the number of bits between the top two hypotheses: ``` 102*(0.27166534-0.02383006) [1] 25.2792 ``` The top hypothesis has 25.2792 bit in support of it over the 2nd alternative as measured by KL. This is very close to the maximum likelihood difference in bits! Min KL is max likelihood. TODO: why not exactly the same. ================================ # Appendix: Multiplying Independent Distributions is Still a Distribution Let's consider multiple Bernoulli trials of flipping a biased coin: $ P(H)=90 \%, P(T)=10 \% $ For a single flip of a single coin it is obviously a distribution: $\sum_{c1} P(C) = 0.9 + 0.1 = 1.0$ For two independent flips we want the joint which is a 2d matrix covering possibilities [[HH, HT],[TH,TT]] where we multiply to Bernoulli together $$P(c1,c2) = P(c1) P(c2) = [[0.9 \times 0.9, 0.9 \times 0.1], [0.9 \times 0.1, 0.1 \times 0.1]]$$ The probabilities get smaller because of the multiplying but we have a larger number of terms in the sum to account for the different coins. It is a distribution: $\sum_{c1,c2} P(c1,c2)= (0.81 + 0.09 + 0.09 + 0.01) = 1.0$ For 3 coins, the joint is a 3d matrix. For $N$ coins, the joint is a N-d matrix which can be very large. For 100 coins there would be 2^100 terms in the sum to evaluate! However, note that the probabilities simply depend on the number of heads and tails for Bernoulli: $$\sum_{c1,c2} P(c1,c2)= 0.9^2 0.1^0+ 0.9^1 0.1^1 + 0.9^1 0.1^1 + 0.9^0 0.1^2 $$ $$= choose(2,2) 0.9^2 0.1^0 + choose(2,1) (0.9 0.1) + choose(2,0) 0.9^0 0.1^2$$ This is exactly the binomial distribution! It sums to 1.0: $$ 1 = 1^N = (p + 1-p)^N = \sum_{k=0}^N choose(N, k) p^k (1-p)^{(N-k)} $$ This is motivation why multiplying independent distributions is still a distribution and sums to 1 over all possibilities.