Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import random
- import matplotlib.pyplot as plt
- import numpy as np
- import time
- import scipy.stats as ss
- a=0 # xmin
- b=1 # xmax
- m=3/2 # ymax
- variables = [] #list for variables
- def f(x):
- return 3/2 * (1 - x**2) #probability density function
- reject = 0 # number of rejections
- start = time.time()
- while len(variables) < 100000: #I want to generate 100 000 variables
- u1 = random.uniform(a,b)
- u2 = random.uniform(0,m)
- if u2 <= f(u1):
- variables.append(u1)
- else:
- reject +=1
- end = time.time()
- print("Time: ", end-start)
- print("Rejection: ", reject)
- x = np.linspace(a,b,1000)
- plt.hist(variables,50, density=1)
- plt.plot(x, f(x))
- plt.show()
- ss.probplot(variables, plot=plt)
- plt.show()
- import numpy as np
- import scipy.stats as ss
- import matplotlib.pyplot as plt
- import time
- a = 0. # xmin
- b = 1. # xmax
- m = 3.0/2.0 # ymax
- def f(x):
- return 1.5 * (1.0 - x*x) # probability density function
- start = time.time()
- N = 100000
- u1 = np.random.uniform(a, b, N)
- u2 = np.random.uniform(0.0, m, N)
- negs = np.empty(N)
- negs.fill(-1)
- variables = np.where(u2 <= f(u1), u1, negs) # accepted samples are positive or 0, rejected are -1
- end = time.time()
- accept = np.extract(variables>=0.0, variables)
- reject = N - len(accept)
- print("Time: ", end-start)
- print("Rejection: ", reject)
- x = np.linspace(a, b, 1000)
- plt.hist(accept, 50, density=True)
- plt.plot(x, f(x))
- plt.show()
- ss.probplot(accept, plot=plt) # against normal distribution
- plt.show()
- class my_pdf(ss.rv_continuous):
- def _pdf(self, x):
- return 1.5 * (1.0 - x*x)
- ss.probplot(accept, dist=my_pdf(a=a, b=b, name='my_pdf'), plot=plt)
- n_before_accept_reject = 150000
- u1 = np.random.uniform(a, b, size=n_before_accept_reject)
- u2 = np.random.uniform(0, m, size=n_before_accept_reject)
- variables = u1[u2 <= f(u1)]
- reject = n_before_accept_reject - len(variables)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement