0
$\begingroup$

Following is the pdf from which I want to sample so, I used roulette wheel samplinguser defined pdf

Code to generate pdf

pdf=np.array([1,5,1,10,1,5,1])
pdf=pdf/np.sum(pdf)
x=np.arange(0,np.size(pdf),1)
plt.bar(x,pdf)

Code to generate cdf

def pdf2cdf(tpdf):
tcdf=np.zeros_like(tpdf)
itr=0
for i in tpdf:
    tcdf[itr]=tcdf[itr-1]+tpdf[itr]
    itr=itr+1
return tcdf
cdf=pdf2cdf(pdf)
plt.bar(x,cdf)

CDF of the given pdfCDF

Code to sample from pdf using roulette wheel

lst=[]
uni=[]
for j in range(100000):
    tmp=np.random.uniform(0,1)
    uni.append(tmp)
    tmp2=np.argsort(np.abs(cdf-tmp))
    lst.append(tmp2[0])

histogram of uni

uni

histogram of samples drawn

samples

I don't understand what's going wrong. The pdf from which I sampled(first plot) and the histogram of samples (last plot) are not the same. Can anybody please suggest the correction of any mistake I made above.


Thanks @jsk

for k in range(np.size(cdf)):
    if cdf[k]>tmp:
        break
tmp2=k

plt.hist(lst,bins=range(0,8), align='left', rwidth=1, normed=True)

enter image description here

$\endgroup$
0
$\begingroup$

The problem is with the command

np.argsort(np.abs(cdf-tmp)).

The command will find which element of the cdf is closest to tmp, regardless of whether or not tmp is greater than cdf.

According to your problem,

$P(0) = 1/24 \approx 0.04167$

$P(1) = 5/24 \approx 0.2083 $

so $ CDF(1) = 6/24 = 0.25$.

In your code, you will end up sampling a 0 if tmp is closer to 0.04167 than 0.25, for example tmp = 0.1 will return a 0, which means that you are assigning 0 a much higher proportion of the time than 0.04167. What range of values for tmp will ensure that you assign a 0 with probability 0.04167? What range of values for tmp will ensure that you will assign a 1 with the appropriate probability?

I will leave it as an exercise for you to figure out how to adjust your code in light of the problem I have identified with the current code.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.