biased.m

一些常被用于教学或者参考的概率论的实例的源代码

%biased.m/created by PJNahin for "Duelling Idiots"(12/27/98)
%This m-file simulates 1,000 games of odd-person-out, for N people as
%defined by the user, with N-1 persons flipping a fair coin and a single person
%person flipping a biased coin (with the probability of heads equal to
%to q). The elements of the row vector duration are the number of games of length
%i, i.e., duration(i)=# of games that require i flips to complete, where i=1,2,3...
%
%
duration=zeros(1,50);         %initialize to zero the number of games
                              %that have durations from 1 to 50 flips
N=input('Enter number of people playing the game:')
q=input('Enter probability of heads for the biased coin:')                              
rand('state',100*sum(clock))  %set new seed for the random number
                              %generator
for game=1:1000               %simulate 1,000 games
   gameover=0;                %gameover=1 is the flag that the
                              %current game is done
   flip=1;                    %initialize duration of current game
  while gameover==0           %keep playing current game
   s=0;                       %initialize number of heads before
                              %each flipping of N coins
   fair=N-1;
   for n=1:fair               %flip a fair coin for each of N-1
      coin=rand;              %people, and increment s by 1 for
      if coin<1/2             %each head
         s=s+1;
      end
   end
   coin=rand;                 %flip the biased coin;
   if coin<q                  %increment s if biased
      s=s+1;                  %coin shows heads;
   end
   if s==1 | s==N-1           %is there an odd-person-out?
         gameover=1;          %yes, set game-over flag
      else
         flip=flip+1;         %no, so flip again
      end                     %check while statement for end of game
   end
      duration(flip)=duration(flip)+1; %game is over, up-date duration vector
end                           %play new game
   average=0;
   for i=1:length(duration)
      average=average+i*duration(i);
   end
average=average/1000;
disp(['The average number of flips/game for ',num2str(N),' players is ',num2str(average)])