GIBBS SAMPLING FOR N DISCRETE VARIABLES IN A N-SPACE

Dear all,

I would like to ask how can be perfomed a Gibbs sampling for a space constituted by n discrete distributions. As an example, I provide 5 distributions with different data, different size, and different p(data), which represent the posterior probability:

- data1 = [50 70 100 130 170 230 300 400]; p(data1) = [0.20 0.05 0.05 0.10 0.09 0.01 0.28 0.22];

- data2 = [1 2 3 4 5 6 7 8]; p(data2) = [0.05 0.1 0.25 0.1 0.05 0.05 0.2 0.2];

- data3 = [1 2 3 4]; p(data3) = [0.3 0.5 0.1 0.1];

- data4 = [1 2]; p(data4) = [0.85 0.15];

- data5 = [1 2]; p(data5) = [0.9 0.1];

[Star Strider]

If you have the Econometrics Toolbox, the sampleroptions function (introduced in R2017b) may be helpful.

This is far from my areas of expertise, so I cannot help with it. (I did not know it even existed before reading your post.) I just searched to see what I could find.

[Sergio RUGGIERI]

Dear Star Strider, thanks for your answer.

Unfortunately, I did not find any example about this sampling technique with discrete distributions, and even less if in the space I consider also continous distributions.

[Star Strider]

My pleasure. I did an Interweb search, and unfortunately found nothing with respect to your question, although you may be able to find something since you know what you’re looking for. What I posted was as close as I can get. You may have to write your own code for this project.

Answer 1

[Vinayak]

Hi Sergio,

It seems like you want to distribute 'x' data points according to their probabilities into 10,000 samples. This process should be repeated for all 'n' distributions, generating a sample object where each object includes all 'n' data points.

You may use the “randsample” function, which takes data, sample size, and probability as inputs to output the sample dataset with the probabilities provided. This function can be used for all independent 'n' distributions separately.

For the example data you provided, the following code should work:

% Your declaration of data as matrices. 
samples = zeros(10000, 5); 
for i = 1:10000 
    samples(i, 1) = randsample(data1, 1, true, pdata1); 
    samples(i, 2) = randsample(data2, 1, true, pdata2); 
    samples(i, 3) = randsample(data3, 1, true, pdata3); 
    samples(i, 4) = randsample(data4, 1, true, pdata4); 
    samples(i, 5) = randsample(data5, 1, true, pdata5); 
end 

If you want to learn more about the "randsample" function, refer to this documentation link - https://www.mathworks.com/help/stats/randsample.html

I hope this helps!

[Sergio RUGGIERI]

Dear Vinayak,

thanks for your answer. This is the unique way I found to have numerically correct results, but I did not take it into account because in this way, the information about the space (i.e., the Joint PMF) should be lost.