simulating an unbias coin with a bias coin flip

Hello, I have a question about achieving an unbias coin with a bias coin flip for N = 1000 tosses, for a set of I have written code for unknown bias probabilities p = [ 0.5 0.4 0.3 0.2 0.1].I would like to choose these arbitrarily to achieve an unbias coin for each. I visited this site to read about the process http://www.pcoder.net/change-an-unfairly-biased-coin-to-a-fair-coin-a-classical-randomized-algorithm/#axzz308a9TpBx but im confuse about the implementation in matlab. thanks in advance for any feedback. thanks

if true
% calling the function
p = [ 0.5 0.4 0.3 0.2 0.1];
N = 1000;
for i = 1:length(p)

   function outcome = mysim(p, N)
P = cumsum(p); 
u = rand(1, N); 
outcome = zeros(1, N); % A blank array for holding the outcomes 
for n=100:100:N, 
   h = find(u(n)<P, 1 );
    outcome(n) = h;
end % code
  end

[Andrew Newell]

Let's see if you can get some code that actually runs first:

1. if true is pointless - you can leave it off
2. The function should be in a separate file, not inside a for loop.
3. In your for loop, you need to call your function (e.g.:

for i = 1:length(p)
    outcome(i) = mysim(p,N);
end

[Genus]

I implemented your recommendation however there is an error at outcome(i) = mysim(p, N); it states the variable appears to change size on every loop iteration (within a script) consider preallocating for speed. What are your reccomendations for rectifying this issue, thanks

[Andrew Newell]

That's not an error, just a suggestion. You can preallocate outcome by adding the line

outcome = zeros(size(p));

before the loop.

Answer 1

[Andrew Newell]

I think that the simplest approach is to write a little function to simulate a single toss, and then build up your simulation from that. Here is that function:

function side = simulateOneToss
% Make two tosses (outcomes 0 and 1 could stand for heads and tails)
twoTosses = round(rand(1,2));

 % If outcome is not HT or TH (both of which sum to 1), try again.
while sum(twoTosses) ~= 1
    twoTosses = round(rand(1,2));
end

 % Take first of the two tosses as the answer.
side = twoTosses(1);

Based on the comments below, here is a version that allows a probability other than 0.5:

function side = simulateOneToss(p)
% Make two tosses (outcomes 0 and 1 could stand for heads and tails)
twoTosses = rand(1,2)<p;

% If outcome is not HT or TH (both of which sum to 1), try again.
while sum(twoTosses) ~= 1
  twoTosses = rand(1,2)<p;
end

% Take first of the two tosses as the answer.
side = twoTosses(1);

Note that I don't round before comparing with the probability (otherwise you'll always get p=0.5 in effect). Now this can be used in a simulator:

function outcomes = simulateTosses(p,nTosses,nTrials)

tosses = zeros(nTosses,nTrials);
for i=1:numel(tosses)
  tosses(i) = simulateOneToss(p);
end

outcomes = sum(tosses>p);

Finally, you can run it and plot the results as follows:

nTosses = 100;
nTrials = 10;
p = 0.6;
outcomes = simulateTosses(p,nTosses,nTrials);
bar(1:nTrials,outcomes)

[Genus]

yes Thanks this works beautifully, however, I was seeking to incorporate my for loop from the first function to implement trials from 100:100:1000. I tried to plot this but my plot chart exhibited a blank plot.

    function outcome1 = simulateToss(p,N)
% Make two tosses (outcomes 0 and 1 could stand for heads and tails)
 P = cumsum(p); 
 twoTosses = round(rand(1,100));
 outcome1 = zeros(1, 100);
 % If outcome is not HT or TH (both of which sum to 1), try again.
 for n=100:100:N
 while sum(twoTosses) ~= 1
     twoTosses = (round(rand(1,100))< P);
end
 % Take first of the two tosses as the answer.
outcome1(n) = twoTosses(1);

 end
  end

[Andrew Newell]

I'm not sure what you're trying to do - simulate 100 tosses at a time?

[Genus]

yes exactly, I'm trying to simulate 100 tosses from 100 to 1000 in steps of 100, I would like to arbitrarily pick a probability for this purpose.

[Andrew Newell]

Why are you using cumsum on p? Is the probability changing?

[Andrew Newell]

For now, I'll assume it doesn't change.

[Genus]

yes, the probablilties should be abritrarily picked form the set of p = [ 0.5 0.4 0.3 0.2 0.1];so they should be changing,I assumed that will be the right statement to abritrarly select probabilities from the set. the statement should be 1000 trials from 100 to 1000 in steps of 100.

[Andrew Newell]

By "arbitrarily" you mean randomly? In your initial attempt at this, you had a for loop where you went through the choices systematically. Do you want to do 1000 trials for each probability and plot the results?

[Genus]

yes i apologise about the confusion about the arbitrarily terminology. I made an error in my initial description. It should be out of an observation of coins i.e 5 or 10 coins. we should be be able to randomly select a coin. and determine if the coin is unbiased. One thing I initially forgot to emphasize is that a coin is described by its bias and defined by the probability of coming up as head. Then a simulation should be done in (two different ways of 1000 trials) in steps of 100 starting from 100 for each random selected coin and a combined 1000 trials as a whole. We should be able to select each probability in the set to do this. Thanks

[Andrew Newell]

Well, if I understand your problem, you can just put the last part of the code in a for loop:

p = [ 0.5 0.4 0.3 0.2 0.1];
nTosses = 100;
nTrials = 10;
for i = 1:length(p)
  outcomes = simulateTosses(p(i),nTosses,nTrials);
  figure; bar(1:nTrials,outcomes)
end

[Genus]

excellent, it works. I have a question about descision or hypothesis testing. suppose I have a conditional probability defined by the binomial probability

    function [p] = calculateInputProbability(N, T)
     inputP = [];

    for i=1:N+1
        combo = nchoosek(N,i-1);
        temp  = combo * power(T , i-1);
        temp1  = temp * power((1-T), N-(i-1));
        inputP(i) = temp1;          
    end

    p  = inputP;

end

im hoping to do hypothesis testing with H0 vs H1 the same specifications as before. but now with theta is in[0,1];. Ho having the conditional probability and theta = 1/2 and for H1 defined by the same conditional probability but theta =[ 0.5 0.4 0.3 0.2 0.1];

1) suppose we have 5 bias coins corresponding to each bias element theta =[ 0.5 0.4 0.3 0.2 0.1] and unbias coin corresponding to theta = 1/2. 2) we would like to choose an arbitrary bias single coin from the 5 coins and flip it 100 times and determine if the coin is bias or unbias. if its unbias probability is 1/2 otherwise bias. 3) we would like for each selected coin to simulate it from 100 to 1000 in steps of 100. . how do i rate you on this site. you have helped me greatly. here is my code so far. thanks

%% Input Values n0 = 0.5; n1 = [0.4 0.3 0.2 0.1]; nTosses = [100 200 300 400 500 600 700 800 900 1000]; nTrials = 1; for i = 1:length(n1) for j = 1: length(nToses) outcomes = simulateTosses(n1(i),nTosses(j),nTrials); end %% Calculate the prbability [p] = calculateInputProbability(nTosses(j), n1(i)); % code end

[Andrew Newell]

Genus, now you're on to a question on hypothesis testing, which should be posted separately. But first, I would recommend reading up on the relevant statistical methods. For example, there is a Wikipedia page Checking whether a coin is fair.

I'm glad I have been able to help you. We don't have ratings as such, but you can click on the "Accept this answer" button.