【Programming using PyMC3】Doing Bayesian Data Analysis by John K. Kruschke

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

本帖隐藏的内容

Doing Bayesian Data Analysis-master.zip (11.59 MB, 需要: 5 个论坛币)

2016-12-11 11:00:54 上传

需要: 5 个论坛币  [购买]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Programming Bayesian Analysis Kruschke Program

1. """
   
2. Goal: Toss a coin N times and compute the running proportion of heads.
   
3. """
   
4. import matplotlib.pyplot as plt
   
5. import numpy as np
   
6. 
   
7. # Specify the total number of flips, denoted N.
   
8. N = 500
   
9. # Generate a random sample of N flips for a fair coin (heads=1, tails=0);
   
10. np.random.seed(47405)
    
11. flip_sequence = np.random.choice(a=(0, 1), p=(.5, .5), size=N, replace=True)
    
12. # Compute the running proportion of heads:
    
13. r = np.cumsum(flip_sequence)
    
14. n = np.linspace(1, N, N)  # n is a vector.
    
15. run_prop = r/n  # component by component division.
    
16. 
    
17. # Graph the running proportion:
    
18. plt.plot(n, run_prop, '-o', )
    
19. plt.xscale('log')  # an alternative to plot() and xscale() is semilogx()
    
20. plt.xlim(1, N)
    
21. plt.ylim(0, 1)
    
22. plt.xlabel('Flip Number')
    
23. plt.ylabel('Proportion Heads')
    
24. plt.title('Running Proportion of Heads')
    
25. # Plot a dotted horizontal line at y=.5, just as a reference line:
    
26. plt.axhline(y=.5, ls='dashed')
    
27. 
    
28. # Display the beginning of the flip sequence.
    
29. flipletters = ''
    
30. for i in flip_sequence[:10]:
    
31.     if i == 1:
    
32.         flipletters += 'H'
    
33.     else:
    
34.         flipletters += 'T'
    
35. 
    
36. plt.text(10, 0.8, 'Flip Sequence = %s...' % flipletters)
    
37. # Display the relative frequency at the end of the sequence.
    
38. plt.text(25, 0.2, 'End Proportion = %s' % run_prop[-1])
    
39. 
    
40. plt.savefig('Figure_3.1.png')

复制代码

1. """
   
2. Inferring a binomial proportion via exact mathematical analysis.
   
3. """
   
4. import sys
   
5. import numpy as np
   
6. from scipy.stats import beta
   
7. from scipy.special import beta as beta_func
   
8. import matplotlib.pyplot as plt
   
9. from HDIofICDF import *
   
10. 
    
11. 
    
12. def bern_beta(prior_shape, data_vec, cred_mass=0.95):
    
13.     """Bayesian updating for Bernoulli likelihood and beta prior.
    
14.      Input arguments:
    
15.        prior_shape
    
16.          vector of parameter values for the prior beta distribution.
    
17.        data_vec
    
18.          vector of 1's and 0's.
    
19.        cred_mass
    
20.          the probability mass of the HDI.
    
21.      Output:
    
22.        post_shape
    
23.          vector of parameter values for the posterior beta distribution.
    
24.      Graphics:
    
25.        Creates a three-panel graph of prior, likelihood, and posterior
    
26.        with highest posterior density interval.
    
27.      Example of use:
    
28.      post_shape = bern_beta(prior_shape=[1,1] , data_vec=[1,0,0,1,1])"""
    
29. 
    
30.     # Check for errors in input arguments:
    
31.     if len(prior_shape) != 2:
    
32.         sys.exit('prior_shape must have two components.')
    
33.     if any([i < 0 for i in prior_shape]):
    
34.         sys.exit('prior_shape components must be positive.')
    
35.     if any([i != 0 and i != 1 for i in data_vec]):
    
36.         sys.exit('data_vec must be a vector of 1s and 0s.')
    
37.     if cred_mass <= 0 or cred_mass >= 1.0:
    
38.         sys.exit('cred_mass must be between 0 and 1.')
    
39. 
    
40.     # Rename the prior shape parameters, for convenience:
    
41.     a = prior_shape[0]
    
42.     b = prior_shape[1]
    
43.     # Create summary values of the data:
    
44.     z = sum(data_vec[data_vec == 1])  # number of 1's in data_vec
    
45.     N = len(data_vec)   # number of flips in data_vec
    
46.     # Compute the posterior shape parameters:
    
47.     post_shape = [a+z, b+N-z]
    
48.     # Compute the evidence, p(D):
    
49.     p_data = beta_func(z+a, N-z+b)/beta_func(a, b)
    
50.     # Construct grid of theta values, used for graphing.
    
51.     bin_width = 0.005  # Arbitrary small value for comb on theta.
    
52.     theta = np.arange(bin_width/2, 1-(bin_width/2)+bin_width, bin_width)
    
53.     # Compute the prior at each value of theta.
    
54.     p_theta = beta.pdf(theta, a, b)
    
55.     # Compute the likelihood of the data at each value of theta.
    
56.     p_data_given_theta = theta**z * (1-theta)**(N-z)
    
57.     # Compute the posterior at each value of theta.
    
58.     post_a = a + z
    
59.     post_b = b+N-z
    
60.     p_theta_given_data = beta.pdf(theta, a+z, b+N-z)
    
61.     # Determine the limits of the highest density interval
    
62.     intervals = HDIofICDF(beta, cred_mass, a=post_shape[0], b=post_shape[1])
    
63. 
    
64.     # Plot the results.
    
65.     plt.figure(figsize=(12, 12))
    
66.     plt.subplots_adjust(hspace=0.7)
    
67. 
    
68.     # Plot the prior.
    
69.     locx = 0.05
    
70.     plt.subplot(3, 1, 1)
    
71.     plt.plot(theta, p_theta)
    
72.     plt.xlim(0, 1)
    
73.     plt.ylim(0, np.max(p_theta)*1.2)
    
74.     plt.xlabel(r'$\theta$')
    
75.     plt.ylabel(r'$P(\theta)$')
    
76.     plt.title('Prior')
    
77.     plt.text(locx, np.max(p_theta)/2, r'beta($\theta$;%s,%s)' % (a, b))
    
78.     # Plot the likelihood:
    
79.     plt.subplot(3, 1, 2)
    
80.     plt.plot(theta, p_data_given_theta)
    
81.     plt.xlim(0, 1)
    
82.     plt.ylim(0, np.max(p_data_given_theta)*1.2)
    
83.     plt.xlabel(r'$\theta$')
    
84.     plt.ylabel(r'$P(D|\theta)$')
    
85.     plt.title('Likelihood')
    
86.     plt.text(locx, np.max(p_data_given_theta)/2, 'Data: z=%s, N=%s' % (z, N))
    
87.     # Plot the posterior:
    
88.     plt.subplot(3, 1, 3)
    
89.     plt.plot(theta, p_theta_given_data)
    
90.     plt.xlim(0, 1)
    
91.     plt.ylim(0, np.max(p_theta_given_data)*1.2)
    
92.     plt.xlabel(r'$\theta$')
    
93.     plt.ylabel(r'$P(\theta|D)$')
    
94.     plt.title('Posterior')
    
95.     locy = np.linspace(0, np.max(p_theta_given_data), 5)
    
96.     plt.text(locx, locy[1], r'beta($\theta$;%s,%s)' % (post_a, post_b))
    
97.     plt.text(locx, locy[2], 'P(D) = %g' % p_data)
    
98.     # Plot the HDI
    
99.     plt.text(locx, locy[3],
    
100.              'Intervals = %.3f - %.3f' % (intervals[0], intervals[1]))
     
101.     plt.fill_between(theta, 0, p_theta_given_data,
     
102.                     where=np.logical_and(theta > intervals[0],
     
103.                     theta < intervals[1]),
     
104.                         color='blue', alpha=0.3)
     
105.     return intervals
     
106. 
     
107. data_vec = np.repeat([1, 0], [11, 3])  # 11 heads, 3 tail
     
108. intervals = bern_beta(prior_shape=[100, 100], data_vec=data_vec)
     
109. plt.savefig('Figure_5.2.png')
     
110. plt.show()

复制代码

1. """
   
2. Posterior predictive check. Examine the veracity of the winning model by
   
3. simulating data sampled from the winning model and see if the simulated data
   
4. 'look like' the actual data.
   
5. """
   
6. import numpy as np
   
7. from scipy.stats import beta
   
8. import matplotlib.pyplot as plt
   
9. 
   
10. # Specify known values of prior and actual data.
    
11. prior_a = 100
    
12. prior_b = 1
    
13. actual_data_Z = 8
    
14. actual_data_N = 12
    
15. # Compute posterior parameter values.
    
16. post_a = prior_a + actual_data_Z
    
17. post_b = prior_b + actual_data_N - actual_data_Z
    
18. # Number of flips in a simulated sample should match the actual sample size:
    
19. sim_sample_size = actual_data_N
    
20. # Designate an arbitrarily large number of simulated samples.
    
21. n_sim_samples = 1000
    
22. # Set aside a vector in which to store the simulation results.
    
23. sim_sample_Z_record = np.zeros(n_sim_samples)
    
24. # Now generate samples from the posterior.
    
25. for sample_idx in range(0, n_sim_samples):
    
26.     # Generate a theta value for the new sample from the posterior.
    
27.     sample_theta = beta.rvs(post_a, post_b)
    
28.     # Generate a sample, using sample_theta.
    
29.     sample_data = np.random.choice([0, 1], p=[1-sample_theta, sample_theta],
    
30.                                   size=sim_sample_size, replace=True)
    
31.     sim_sample_Z_record[sample_idx] = sum(sample_data)
    
32. 
    
33. 
    
34. ## Make a histogram of the number of heads in the samples.
    
35. plt.hist(sim_sample_Z_record)
    
36. plt.show()

复制代码

1. """
   
2. Use this program as a template for experimenting with the Metropolis algorithm
   
3. applied to 2 parameters called theta1,theta2 defined on the domain [0,1]x[0,1].
   
4. """
   
5. from __future__ import division
   
6. import numpy as np
   
7. from scipy.stats import beta
   
8. import matplotlib.pyplot as plt
   
9. 
   
10. 
    
11. # Define the likelihood function.
    
12. # The input argument is a vector: theta = [theta1 , theta2]
    
13. 
    
14. def likelihood(theta):
    
15.     # Data are constants, specified here:
    
16.     z1, N1, z2, N2 = 5, 7, 2, 7
    
17.     likelihood = (theta[0]**z1 * (1-theta[0])**(N1-z1)
    
18.                  * theta[1]**z2 * (1-theta[1])**(N2-z2))
    
19.     return likelihood
    
20. 
    
21. 
    
22. # Define the prior density function.
    
23. # The input argument is a vector: theta = [theta1 , theta2]
    
24. def prior(theta):
    
25.     # Here's a beta-beta prior:
    
26.     a1, b1, a2, b2 = 3, 3, 3, 3
    
27.     prior = beta.pdf(theta[0], a1, b1) * beta.pdf(theta[1], a2, b2)
    
28.     return prior
    
29. 
    
30. 
    
31. # Define the relative probability of the target distribution, as a function
    
32. # of theta.  The input argument is a vector: theta = [theta1 , theta2].
    
33. # For our purposes, the value returned is the UNnormalized posterior prob.
    
34. def target_rel_prob(theta):
    
35.     if ((theta >= 0.0).all() & (theta <= 1.0).all()):
    
36.         target_rel_probVal =  likelihood(theta) * prior(theta)
    
37.     else:
    
38.         # This part is important so that the Metropolis algorithm
    
39.         # never accepts a jump to an invalid parameter value.
    
40.         target_rel_probVal = 0.0
    
41.     return target_rel_probVal
    
42. #    if ( all( theta >= 0.0 ) & all( theta <= 1.0 ) ) {
    
43. #        target_rel_probVal =  likelihood( theta ) * prior( theta )
    
44. 
    
45. 
    
46. # Specify the length of the trajectory, i.e., the number of jumps to try:
    
47. traj_length = 5000 # arbitrary large number
    
48. # Initialize the vector that will store the results.
    
49. trajectory = np.zeros((traj_length, 2))
    
50. # Specify where to start the trajectory
    
51. trajectory[0, ] = [0.50, 0.50] # arbitrary start values of the two param's
    
52. # Specify the burn-in period.
    
53. burn_in = np.ceil(.1 * traj_length) # arbitrary number
    
54. # Initialize accepted, rejected counters, just to monitor performance.
    
55. n_accepted = 0
    
56. n_rejected = 0
    
57. # Specify the seed, so the trajectory can be reproduced.
    
58. np.random.seed(47405)
    
59. # Specify the covariance matrix for multivariate normal proposal distribution.
    
60. n_dim, sd1, sd2 = 2, 0.2, 0.2
    
61. covar_mat = [[sd1**2, 0], [0, sd2**2]]
    
62. 
    
63. # Now generate the random walk. step is the step in the walk.
    
64. for step in range(traj_length-1):
    
65.     current_position = trajectory[step, ]
    
66.     # Use the proposal distribution to generate a proposed jump.
    
67.     # The shape and variance of the proposal distribution can be changed
    
68.     # to whatever you think is appropriate for the target distribution.
    
69.     proposed_jump = np.random.multivariate_normal(mean=np.zeros((n_dim)),
    
70.                                                  cov=covar_mat)
    
71.     # Compute the probability of accepting the proposed jump.
    
72.     prob_accept = np.minimum(1, target_rel_prob(current_position + proposed_jump)
    
73.                             / target_rel_prob(current_position))
    
74.     # Generate a random uniform value from the interval [0,1] to
    
75.     # decide whether or not to accept the proposed jump.
    
76.     if np.random.rand() < prob_accept:
    
77.         # accept the proposed jump
    
78.         trajectory[step+1, ] = current_position + proposed_jump
    
79.         # increment the accepted counter, just to monitor performance
    
80.         if step > burn_in:
    
81.             n_accepted += 1
    
82.     else:
    
83.         # reject the proposed jump, stay at current position
    
84.         trajectory[step+1, ] = current_position
    
85.         # increment the rejected counter, just to monitor performance
    
86.         if step > burn_in:
    
87.             n_rejected += 1
    
88. 
    
89. # End of Metropolis algorithm.
    
90. 
    
91. #-----------------------------------------------------------------------
    
92. # Begin making inferences by using the sample generated by the
    
93. # Metropolis algorithm.
    
94. 
    
95. # Extract just the post-burnIn portion of the trajectory.
    
96. accepted_traj = trajectory[burn_in:]
    
97. 
    
98. # Compute the means of the accepted points.
    
99. mean_traj =  np.mean(accepted_traj, axis=0)
    
100. # Compute the standard deviations of the accepted points.
     
101. stdTraj =  np.std(accepted_traj, axis=0)
     
102. 
     
103. # Plot the trajectory of the last 500 sampled values.
     
104. plt.plot(accepted_traj[:,0], accepted_traj[:,1], marker='o', alpha=0.3)
     
105. plt.xlim(0, 1)
     
106. plt.ylim(0, 1)
     
107. plt.xlabel(r'$\theta1$')
     
108. plt.ylabel(r'$\theta2$')
     
109. 
     
110. # Display means in plot.
     
111. plt.plot(0, label='M = %.3f, %.3f' % (mean_traj[0], mean_traj[1]), alpha=0.0)
     
112. # Display rejected/accepted ratio in the plot.
     
113. plt.plot(0, label=r'$N_{pro}=%s$ $\frac{N_{acc}}{N_{pro}} = %.3f$' % (len(accepted_traj), (n_accepted/len(accepted_traj))), alpha=0)
     
114. 
     
115. # Evidence for model, p(D).
     
116. # Compute a,b parameters for beta distribution that has the same mean
     
117. # and stdev as the sample from the posterior. This is a useful choice
     
118. # when the likelihood function is binomial.
     
119. a =   mean_traj * ((mean_traj*(1-mean_traj)/stdTraj**2) - np.ones(n_dim))
     
120. b = (1-mean_traj) * ( (mean_traj*(1-mean_traj)/stdTraj**2) - np.ones(n_dim))
     
121. # For every theta value in the posterior sample, compute
     
122. # beta.pdf(theta, a, b) / likelihood(theta) * prior(theta)
     
123. # This computation assumes that likelihood and prior are properly normalized,
     
124. # i.e., not just relative probabilities.
     
125. 
     
126. wtd_evid = np.zeros(np.shape(accepted_traj)[0])
     
127. for idx in range(np.shape(accepted_traj)[0]):
     
128.     wtd_evid[idx] = (beta.pdf(accepted_traj[idx,0],a[0],b[0] )
     
129.         * beta.pdf(accepted_traj[idx,1],a[1],b[1]) /
     
130.         (likelihood(accepted_traj[idx,]) * prior(accepted_traj[idx,])))
     
131. 
     
132. p_data = 1 / np.mean(wtd_evid)
     
133. # Display p(D) in the graph
     
134. plt.plot(0, label='p(D) = %.3e' % p_data, alpha=0)
     
135. plt.legend(loc='upper left')
     
136. plt.savefig('Figure_8.3.png')
     
137. 
     
138. # Estimate highest density region by evaluating posterior at each point.
     
139. accepted_traj = trajectory[burn_in:]
     
140. npts = np.shape(accepted_traj)[0]
     
141. post_prob = np.zeros((npts))
     
142. for ptIdx in range(npts):
     
143.     post_prob[ptIdx] = target_rel_prob(accepted_traj[ptIdx,])
     
144. 
     
145. # Determine the level at which credmass points are above:
     
146. credmass = 0.95
     
147. waterline = np.percentile(post_prob, (credmass))
     
148. 
     
149. HDI_points = accepted_traj[post_prob > waterline, ]
     
150. 
     
151. plt.figure()
     
152. plt.plot(HDI_points[:,0], HDI_points[:,1], 'ro')
     
153. plt.xlim(0,1)
     
154. plt.ylim(0,1)
     
155. plt.xlabel(r'$\theta1$')
     
156. plt.ylabel(r'$\theta2$')
     
157. 
     
158. # Display means in plot.
     
159. plt.plot(0, label='M = %.3f, %.3f' % (mean_traj[0], mean_traj[1]), alpha=0.0)
     
160. # Display rejected/accepted ratio in the plot.
     
161. plt.plot(0, label=r'$N_{pro}=%s$ $\frac{N_{acc}}{N_{pro}} = %.3f$' % (len(accepted_traj), (n_accepted/len(accepted_traj))), alpha=0)
     
162. # Display p(D) in the graph
     
163. plt.plot(0, label='p(D) = %.3e' % p_data, alpha=0)
     
164. plt.legend(loc='upper left')
     
165. 
     
166. plt.savefig('Figure_8.3_HDI.png')
     
167. 
     
168. plt.show()

复制代码

......

好东西，可以看一看

好东西，可以看一看。

回复过了