Untitled

def Gibbs(nr_samples, P_cond,j = 0, samples_shape = (11,11)):

    samples = np.zeros(samples_shape)
    for _ in range(nr_samples):

        X = P_cond(j)
        i = np.random.choice([i for i in range(len(X))],p=X)
        X = P_cond(i)
        j = np.random.choice([i for i in range(len(X))],p=X)
        samples[i,j] += 1

    return samples

m = 10
def P_cond(j):
    X = []
    denominator = 0

    for x in range(m - j + 1):
        denominator += truncatePoiss(4,x)

    for i in range(m - j + 1):
        X += [truncatePoiss(4,i)/denominator]

    return X

samples_gibbs = Gibbs(100000,P_cond)
P_empir_gibbs = samples_gibbs/samples.sum()