As a question in an assignment, I have been asked to write an Octave function that simulates 1000 experiments of finding a Random Variable X with alphabet {0, 1, 2, 3} and pmf:
Px(0) = 1/8
Px(1) = 1/4
Px(2) = 1/2
Px(3) = 1/8
by asking a sequence of binary, "yes" or "no" questions.
I have determined that the optimal sequence of binary questions asked to find the value of X is to simply ask "Is X = p?" where p are the possible values, in order of decreasing probabilty.
So the optimal sequence would be:
Is X = 2?
If not:
Is X = 1?
If not:
Is X = 0?
If not then X = 3
This is the function I have written:
function x = guessing_experiment(probabilities, n)
% generates n simulations of finding a random number in an alphabet by asking binary questions,
% where 'probabilities' is a list of the probabilities per number in the order the questions will be asked
num_Qs = zeros(1,n); % allocate array of size n for number of questions asked per experiment
[num_col, alphabet_size] = size(probabilities); % get size of alphabet
for i = 1:n % generate n experiments
Qs = 0; % number of questions asked in this experiment
for j = 1:alphabet_size - 1 % iterate through questions
question = rand; % generate random number in range [0, 1]
Qs++; % incremenet number of questions asked
if (question <= probabilities(j)) % if question produces a "yes" answer
break;
endif
endfor
num_Qs(i) = Qs; % store number of questions asked for this experiment
endfor
x = mean(num_Qs); % calculate mean number of questions asked over the n experiments
end
Which is called as guessing_experiment([1/2, 1/4, 1/8, 1/8], 1000)
where the array is the probability of each question producing a "yes" answer, in order of how they are to be asked, and n is the number of experiments.
Asking these questions should produce an average number of questions of 1.75, but my program is always producing an average of ~1.87. Where is my script in error?
I am assuming it has something to do with generating a new random number to simulate each of the 3 questions being asked.
I deleted my previous wrong answer which was stated the your script is right and your calculation is wrong. I think about it once more and your calculation is right. I tried it myself using the following MATLAB script:
% probabilities for each number
p = [1/8,1/4,1/2,1/8];
% sort them from higher to lower
p = sort(p,'descend');
% number of questions per probability
nq = 1:length(p)-1;
% the last question can distinguish between two variables
nq(end+1) = nq(end);
% number of trials
n = 100000;
% random sample number of questions
q = randsample(nq,n,true,p);
% mean number of questions
avgQ = mean(q)
and the obtained avg. was 1.75 - as you calculated. I will try to take another look on your code to see what's wrong
EDIT
the problem with your script is that you ignored the conditional probability, i.e., when asking questions about the variable you ignored the information you already know about it. for example, when you're asking your third question, the probability that the value will be 0 is not p=1/8 but p=1/2 because you already know that it is not 1 or 2.
The fix you need to do is as dividing the probability by the possible events probabilities probabilities(j)/sum(probabilities(j:end)):
n = 10000;
p = [1/8,1/4,1/2,1/8];
% sort them from higher to lower
probabilities = sort(p,'descend');
probabilities(end-1) = probabilities(end-1) + probabilities(end);
probabilities(end) = [];
alphabet_size = numel(probabilities);
num_Qs = zeros(1,n); % allocate array of size n for number of questions asked per experiment
for i = 1:n % generate n experiments
Qs = 0; % number of questions asked in this experiment
for j = 1:alphabet_size % iterate through questions
question = rand; % generate random number in range [0, 1]
Qs = Qs + 1; % incremenet number of questions asked
if question < probabilities(j)/sum(probabilities(j:end)) % if question produces a "yes" answer
break;
end
end
num_Qs(i) = Qs; % store number of questions asked for this experiment
end
x = mean(num_Qs)
x ~ 1.75
the vector of conditional probabilities in this scenario is:
p = [1/8,1/4,1/2,1/8];
p = sort(p,'descend');
cond_p = p./cumsum(p,'reverse')
cond_p =
0.5000 0.5000 0.5000 1.0000