| # Install on Terminal of MacOS #pip3 install -U numpy #pip3 install -U scipy #pip3 install -U matplotlib |
1_MacOS_Terminal.txt
| ########## Run Terminal on MacOS and execute ### TO UPDATE cd "YOUR_WORKING_DIRECTORY" |
Python files
bayesian01.py
| ########## Bayesian Statistics Example for Beginners in Python ########## # ##### Reference #https://www.quantstart.com/articles/Bayesian-Statistics-A-Beginners-Guide/ ##### Bayes' Rule for Bayesian Inference ''' P(θ|D) * P(D) = P(D|θ) * P(θ) P(θ|D) = P(D|θ) * P(θ) / P(D) where: P(θ) is the prior (belief). This is the strength in our belief of θ without considering the evidence D. Our prior view on the probability of how fair the coin is. P(θ|D) is the (updated) posterior (belief). This is the (refined) strength of our belief of θ once the evidence D has been taken into account. After seeing 4 heads out of 8 flips, say, this is our updated view on the fairness of the coin. P(D|θ) is the likelihood. This is the probability of seeing the data D as generated by a model with parameter θ. If we knew the coin was fair, this tells us the probability of seeing a number of heads in a particular number of flips. P(D) is the evidence. This is the probability of the data as determined by summing (or integrating) across all possible values of θ, weighted by how strongly we believe in those particular values of θ. If we had multiple views of what the fairness of the coin is (but didn't know for sure), then this tells us the probability of seeing a certain sequence of flips for all possibilities of our belief in the coin's fairness. ''' ##### import import numpy as np from scipy import stats from matplotlib import pyplot as plt ##### Coin-Flipping Example # #The coin will actually be fair, but we won't learn this until the trials are carried out. #At the start we have no prior belief on the fairness of the coin, that is, we can say that any level of fairness is equally likely. #We are going to use a Bayesian updating procedure to go from our prior beliefs to posterior beliefs as we observe new coin flips. This is carried out using a particularly mathematically succinct procedure using conjugate priors. We won't go into any detail on conjugate priors within this article, as it will form the basis of the next article on Bayesian inference. It will however provide us with the means of explaining how the coin flip example is carried out in practice.# # #See the figure: #In the following figure we can see 6 particular points at which we have carried out a number of Bernoulli trials (coin flips). # #1. (top left) #In the first sub-plot we have carried out no trials and hence our probability density function (in this case our prior density) is the uniform distribution. #It states that we have equal belief in all values of θ representing the fairness of the coin. # #2. (top right) #The next panel shows 2 trials carried out and they both come up heads. #(This is the case of 2 heads. Note that this result can be different, such as, 1 head and 1 tail or 2 tails, rather than 2 heads.) #Our Bayesian procedure using the conjugate Beta distributions now allows us to update to a posterior density. #Notice how the weight of the density is now shifted to the right hand side of the chart. #This indicates that our prior belief of equal likelihood of fairness of the coin, coupled with 2 new data points, #leads us to believe that the coin is more likely to be unfair (biased towards heads) than it is tails. # # #The following two panels show 10 and 20 trials respectively. # #3. (middle left) #Notice that even though we have seen 2 tails in 10 trials we are still of the belief that the coin is likely to be unfair and biased towards heads. # #4. (middle right) #After 20 trials, we have seen a few more tails appear. The density of the probability has now shifted closer to #θ = P(H) = 0.5. Hence we are now starting to believe that the coin is possibly fair. # # #After 50 and 500 trials respectively, we are now beginning to believe that the fairness of the coin is very likely to be around #θ = 0.5 # #5. (bottom left) # #6. (bottom right) #This is indicated by the shrinking width of the probability density, which is now clustered tightly around θ=0.46 in the final panel. #Were we to carry out another 500 trials (since the coin is actually fair) #we would see this probability density become even tighter and centred closer to #θ=0.5. if __name__ == "__main__": # Create a list of the number of coin tosses ("Bernoulli trials") number_of_trials = [0, 2, 10, 20, 50, 500] # Conduct 500 coin tosses and output into a list of 0s and 1s # where 0 represents a tail and 1 represents a head data = stats.bernoulli.rvs(0.5, size=number_of_trials[-1]) # #If you would like to think about an unfair coin with a bias heavily towards heads, say, 0.80 instead of 0.50, then it goes like this: #data = stats.bernoulli.rvs(0.80, size=number_of_trials[-1]) # Discretise the x-axis into 100 separate plotting points x = np.linspace(0, 1, 100) # Loops over the number_of_trials list to continually add # more coin toss data. For each new set of data, we update # our (current) prior belief to be a new posterior. This is # carried out using what is known as the Beta-Binomial model. # For the time being, we won't worry about this too much. It # will be the subject of a later article! for i, N in enumerate(number_of_trials): # Accumulate the total number of heads for this # particular Bayesian update heads = data[:N].sum() # Create an axes subplot for each update ax = plt.subplot(len(number_of_trials) / 2, 2, i + 1) ax.set_title("%s trials, %s heads" % (N, heads)) # Add labels to both axes and hide labels on y-axis plt.xlabel("$P(H)$, Probability of Heads") plt.ylabel("Density") if i == 0: plt.ylim([0.0, 2.0]) plt.setp(ax.get_yticklabels(), visible=False) # Create and plot a Beta distribution to represent the # posterior belief in fairness of the coin. y = stats.beta.pdf(x, 1 + heads, 1 + N - heads) plt.plot(x, y, label="observe %d tosses,\n %d heads" % (N, heads)) plt.fill_between(x, 0, y, color="#aaaadd", alpha=0.5) # Expand plot to cover full width/height and show it plt.tight_layout() plt.savefig("Figure_Bayesian_Statistics_Example_for_Beginners.png") plt.show() |
Figures
Figure_Bayesian_Statistics_Example_for_Beginners.png
qs-bayes-bernoulli.png
(an image file from the following Reference)
No comments:
Post a Comment