1 import numpy as np
 2 import matplotlib.pyplot as plt
 3
 4 # probability distribution we're trying to calculate
 5 p = lambda x: np.exp(-x)
 6
 7 # CDF of p
 8 CDF = lambda x: 1-np.exp(-x)
 9
10 # invert the CDF
11 invCDF = lambda x: -np.log(1-x)
12
13 # domain limits
14 xmin = 0 # the lower limit of our domain
15 xmax = 6 # the upper limit of our domain
16
17 # range limits
18 rmin = CDF(xmin)
19 rmax = CDF(xmax)
20
21 N = 10000 # the total of samples we wish to generate
22
23 # generate uniform samples in our range then invert the CDF
24 # to get samples of our target distribution
25 R = np.random.uniform(rmin, rmax, N)
26 X = invCDF(R)
27
28 # get the histogram info
29 hinfo = np.histogram(X,100)
30
31 # plot the histogram
32 plt.hist(X,bins=100, label=u'Samples');
33
34 # plot our (normalized) function
35 xvals=np.linspace(xmin, xmax, 1000)
36 plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')
37
38 # turn on the legend
39 plt.legend()
40 plt.show()

 1 # -*- coding: utf-8 -*-
 2 '''
 3 The following code produces samples that follow the distribution P(x)=e^−x
 4 for x=[0, 10] and generates a histogram of the sampled distribution.
 5 '''
 6 import numpy as np
 7 import matplotlib.pyplot as plt
 8
 9
10 P = lambda x: np.exp(-x)
11
12 # domain limits
13 xmin = 0  # the lower limit of our domain
14 xmax = 10 # the upper limit of our domain
15
16 # range limit (supremum) for y
17 ymax = 1
18
19 N = 10000    # the total of samples we wish to generate
20 accepted = 0 # the number of accepted samples
21 samples = np.zeros(N)
22 count = 0    # the total count of proposals
23
24 # generation loop
25 while (accepted < N):
26
27     # pick a uniform number on [xmin, xmax) (e.g. 0...10)
28     x = np.random.uniform(xmin, xmax)
29
30     # pick a uniform number on [0, ymax)
31     y = np.random.uniform(0,ymax)
32
33     # Do the accept/reject comparison
34     if y < P(x):
35         samples[accepted] = x
36         accepted += 1
37
38     count +=1
39
40 print count, accepted
41
42 # get the histogram info
43 # If bins is an int, it defines the number of equal-width bins in the given range
44 (n, bins)= np.histogram(samples, bins=30) # Returns: n-The values of the histogram,n是直方图中柱子的高度
45
46 # plot the histogram
47 plt.hist(samples,bins=30,label=u'Samples')   # bins=30即直方图中有30根柱子
48
49 # plot our (normalized) function
50 xvals=np.linspace(xmin, xmax, 1000)
51 plt.plot(xvals, n[0]*P(xvals), 'r', label=u'P(x)')
52
53 # turn on the legend
54 plt.legend()
55 plt.show()

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3
 4 p = lambda x: np.exp(-x)         # our distribution
 5 g = lambda x: 1/(x+1)            # our proposal pdf (we're choosing k to be 1)
 6 CDFg = lambda x: np.log(x +1)    # generates our proposal using inverse sampling
 7
 8 # domain limits
 9 xmin = 0  # the lower limit of our domain
10 xmax = 10 # the upper limit of our domain
11
12 # range limits for inverse sampling
13 umin = CDFg(xmin)
14 umax = CDFg(xmax)
15
16 N = 10000 # the total of samples we wish to generate
17 accepted = 0 # the number of accepted samples
18 samples = np.zeros(N)
19 count = 0 # the total count of proposals
20
21 # generation loop
22 while (accepted < N):
23
24     # Sample from g using inverse sampling
25     u = np.random.uniform(umin, umax)
26     xproposal = np.exp(u) - 1
27
28     # pick a uniform number on [0, 1)
29     y = np.random.uniform(0, 1)
30
31     # Do the accept/reject comparison
32     if y < p(xproposal)/g(xproposal):
33         samples[accepted] = xproposal
34         accepted += 1
35
36     count +=1
37
38 print count, accepted
39
40 # get the histogram info
41 hinfo = np.histogram(samples,50)
42
43 # plot the histogram
44 plt.hist(samples,bins=50, label=u'Samples');
45
46 # plot our (normalized) function
47 xvals=np.linspace(xmin, xmax, 1000)
48 plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')
49
50 # turn on the legend
51 plt.legend()
52 plt.show()