CODEWITHSUNDEEP
pythonnumpy
Ana Branco
Ana Branco

Reputation: 11

Inverse method for discrete variables in Python

I am trying to do the inverse method for a discrete variable, but I don't get anything (if I print the sample I get only one number, not a sample). How do I make this work correctly?

import numpy as np
import numpy.random as rd
import matplotlib.pyplot as plt
def inversao_discreto(p,M):
    amostra=np.zeros(M)
    acum=np.zeros(len(p))
    acum[:]=p
    for i in range(1, len(acum)):
        acum[i]+=acum[i-1]
    r=rd.random_sample(M)
    for i in range(M):
        i=0
        k=0
        while(r[k]>acum[i]):
            i+=1
            amostra=i
    return amostra
model=np.array([0.1,0.2,0.1,0.2,0.4])
sample=inversao_discreto(model,10000)

Upvotes: 1

Views: 592

Answers (1)

Matheus Portela
Matheus Portela

Reputation: 2460

As far as I could understand, you want to implement the Inverse transform sampling for discrete variables, which works like this:

Given a probability distribution `p`
Given a number of samples `M`
Calculate the cumulative probability distribution function `acum` from `p`
Generate `M` uniformly distributed samples and store into `r`
For each sample `r[i]`
    Get the index where the cumulative distribution `acum` is larger or equal to the sample `r[i]`
    Store this value into `amostra[i]`

If my understanding is correct, then your code is almost there. Your mistake is only in the last for loop. Your variable i is tracking the position of your samples stored in r. Then, k will keep track of where in the accumulator acum you are comparing r to. Here is the pseudo-code:

For each sample `r[i]`
    Start with `k = 0`
    While `r[i] > acum[k]`
        Increment `k`
    Now that `r[i] <= acum[k]`, store `k` into `amostra[i]`

Translating it to Python:

for i in range(M):
    k = 0

    while (r[i] > acum[k]):
        k += 1

    amostra[i] = k

So here is the code fixed:

import numpy as np
import numpy.random as rd

def inversao_discreto(p, M):
    amostra = np.zeros(M)
    acum = np.zeros(len(p))
    acum[:] = p

    for i in range(1, len(acum)):
        acum[i] += acum[i - 1]

    r = rd.random_sample(M)

    for i in range(M):
        k = 0

        while (r[i] > acum[k]):
            k += 1

        amostra[i] = k

    return amostra

model = np.array([0.1, 0.2, 0.1, 0.2, 0.4])
sample = inversao_discreto(model, 10000)

I tried to change as little code as possible to make it work since you are relatively new to Python. Most of the changes I implemented were based on this Style Guide for Python Code, which I recommend you to take a look at since will improve your code, visually speaking.

Upvotes: 1

Related Questions