How to simulate a coin flip experiment and count the number of changeovers?

I have to create an experiment where a fair coin is flipped 20 times and X is the number of times it goes from Head to Tail or Tail to Head. I know the probability of a changeover is 0.5 and the maximum number of changeovers is 19 but I don't know to create the experiment. I don't know how to add a condition that there is a changeover or increment the count(no_Changeovers+1) by 1

[James Tursa]

Please post the code you have written so far and ask specific questions about it and then we can advise you .

[Kai Feng Gan]

This is the code so far:

N=20;
nreps=10000;
count=zeros(1,N);

for n=1:nreps
     no_Changeovers=0
        ntosses=N
      for i=1:N
        if (have to add condition that there is a changeover here)
        end
      end
     (Have to add code to increment count (no_Changeovers+1) by 1
end

Answer 1

[KSSV]

N = 20 ;
T = 0 ; H = 0 ;
for i = 1:N
    c = rand ;
    if c > 0.5
        T = T+1 ;
        fprintf('Tails\n')
    else
        H = H+1 ;
        fprintf('Heads\n')
    end
end

fprintf('In %d tosses,number of Heads=%d and number of Tails = %d\n',N,H,T)

[Kai Feng Gan]

Thank you for the help but unless I typed your code in wrongly, it only gave me the number of heads and tails. If it's not too much trouble, how do I get the number of changeovers(Lets call this X) and the mean and variance of X?

Answer 2

[John D'Errico]

Well, having posted your code, let me suggest the MATLAB way.

simsize = 20; % 20 samples per simulation
nsims = 10000; % number of simulated runs
% flip a coin
coins = rand(simsize,nsims) > 0.5;

So at this point, we have ALL of the necessary data in the array coins. Call a 0 tails, and a 1 heads.

When does a "changeover" occur? It happens at a transition between 0 and 1, or 1 and 0. How can we find them efficiently?

What would diff do here? THINK ABOUT IT! If the coin stays heads of tails, then the diff will be zero.

But if we go from 0 to 1, then the diff will be 1. If we go from 1 to 0, then the diff will be -1. How can we treat them both as the same thing?

D = abs(diff(coins,[],1));

Now, how can you count the number of changeover events in each sim?

nchanges = sum(D,1);

I suppose we might want to count the number of times there were 0,1,2,3,4,..,19 changes.

Change_freq = accumarray(1 + nchanges',ones(nsims,1),[20,1],@sum);

[(0:19)',Change_freq,Change_freq/nsims*100]
ans =
          0            0            0
          1            0            0
          2            4         0.04
          3           20          0.2
          4           75         0.75
          5          214         2.14
          6          514         5.14
          7          981         9.81
          8         1442        14.42
          9         1785        17.85
         10         1775        17.75
         11         1447        14.47
         12          932         9.32
         13          478         4.78
         14          232         2.32
         15           73         0.73
         16           24         0.24
         17            3         0.03
         18            1         0.01
         19            0            0

So the first column is the number of changeovers. The second column the count for that, and column 3 is the percentage of the time we saw that many changeovers.

Could this have been done differently? Well, yes. I might have used a Markov chain, where the state of the process is how many changeovers have been observed. Or, I might use some theory, viewing this properly as a binomial process.

Lets first try to compute the analytical distribution of the number of observed changeovers. I'll do it using a Markov chain.

Assume that the state of our Markov process is the number of changeovers we have seen after n coin tosses. The initial state of the system is

C = [1,zeros(1,19)];

So C(i) represents a probability that we have seen i-1 changeovers. At time 1, we have seen only one coin toss, so the initial state is 0 changeovers, with probability 1.

Now, create a Markov transition matrix, that will see a change from any state to the next higher state with probability 0.5, or you will stay in the current state with probability 0.5.

T = diag(repelem([.5,1],[19,1])) + diag(repelem(.5,19),1);

Look at T. think about what it means. Here are the first 5 rows and columns.

T(1:5,1:5)
ans =
        0.5          0.5            0            0            0
          0          0.5          0.5            0            0
          0            0          0.5          0.5            0
          0            0            0          0.5          0.5
          0            0            0            0          0.5

Now we can compute the state of our system after the second coin toss, just by multiplying by T.

C*T
ans =
  Columns 1 through 14
          0.5          0.5            0            0            0            0            0            0            0            0            0            0            0            0
  Columns 15 through 20
            0            0            0            0            0            0

After the third coin toss?

C*T^2
ans =
  Columns 1 through 14
         0.25          0.5         0.25            0            0            0            0            0            0            0            0            0            0            0
  Columns 15 through 20
            0            0            0            0            0            0

After the fourth coin toss?

C*T^3
ans =
  Columns 1 through 14
        0.125        0.375        0.375        0.125            0            0            0            0            0            0            0            0            0            0
  Columns 15 through 20
            0            0            0            0            0            0

At the 20'th coin toss? Remember, we don't count the first coin toss, since there can have been no changeovers after the very first coin toss. So we raise T to only the 19th power here.

C*T^19
ans =
  Columns 1 through 14
   1.9073e-06    3.624e-05   0.00032616    0.0018482    0.0073929     0.022179      0.05175     0.096107      0.14416       0.1762       0.1762      0.14416     0.096107      0.05175
  Columns 15 through 20
     0.022179    0.0073929    0.0018482   0.00032616    3.624e-05   1.9073e-06

Lets display this in the same form as what we had before. The second column here is the analytical fraction of time we expect to see that number of changes. Yes, I said analytical. This is not a Monte Carlo simulation as was done above.

[(0:19)',(C*T^19)'*100]
ans =
          0   0.00019073
          1     0.003624
          2     0.032616
          3      0.18482
          4      0.73929
          5       2.2179
          6        5.175
          7       9.6107
          8       14.416
          9        17.62
         10        17.62
         11       14.416
         12       9.6107
         13        5.175
         14       2.2179
         15      0.73929
         16      0.18482
         17     0.032616
         18     0.003624
         19   0.00019073

We can do more yet. That is indeed the fraction of times we will see any number of changeovers after 20 coin tosses, but what is the mean? Can we convert this to a statistical distribution? We might decide that this looks quite normally distributed. And while a normal distribution might make sense, perhaps there is a better choice of distribution.

Suppose we think about this not in terms of coin tosses, but purely in terms of did the coin flip sides from what it was last time we tossed it? If we think about it in those terms, we should recognize this as a binomial random process, of length 19. Each time the coin is flipped, we are not worrying about whether we see a head or tail, but did it change direction! So there are 19 such events we will observe.

This is indeed a binomial process. You just need to view it properly as such. What is the mean of a binomial random variable with p=0.5, and 19 samples?

https://en.wikipedia.org/wiki/Binomial_distribution

The mean of a binomial process with n samples and probability p is n*p.

n = 19;
p = 0.5;
n*p
ans =
        9.5

See that this looks absolutely consistent with the distribution we have observed, both from the Monte Carlo, as well as the Markov chain model we used.

A nice thing is we can use the above to predict the distribution of changeovers after 1000 tosses. Don't forget to subtract 1, because we cannot have a changeover until the second coin is tossed.

999*0.5
ans =
      499.5

So after the coin is tossed 1000 total times, we expect to see a mean number of changeovers of 499.5.

The standard deviation of that distribution will be

sqrt(n*p*(1-p))

And as the number of samples grows, it will asymptotically approach a Normal (Gaussian) distribution.

[Kai Feng Gan]

Thank you for the help, I really appreciate the amount of effort you put into your post. However, the question I'm doing requires that I enter the code using the format I posted above and if it's not too much trouble, I was wondering if you write a code using that format. I'm a complete beginner at Matlab so I apologise if it should be easy to write my code using what you posted.