how to do simulation for branching process?

From: lucy (losemind_at_yahoo.com)
Date: 01/25/05 

Date: Mon, 24 Jan 2005 21:47:03 -0800

Hi all,

I am interested in finding the extinction probability for a branching
process. At time 0, there is 1 red cell. At the end of each time period, a
red cell produces

2 red cells with probability 1/4;
1 red cell + 1 white cell with probability 2/3;
2 white cells with probability 1/12.

The life span of both red cell and white cell is one time period. A white
cell does not produce at all.

Now I have already find the extinction probability for the whole population
starting from time 0 with that 1 red cell. I found the extinction
probability = 1/3.

I want to try some Matlab simulation to verify my results.

I found this simulation is not easy to do, because how to decide the number
of runs "R" for this simulation? In each run, how to decide the number of
generations "G"?

If "R" and "G" are too small, the result is not representative: for example,
when R=G=100, I got 25 zeros out of 100 runs, so the probability of
extinction is 1/4. In fact, I have just simulated the probalility of
extinction by 100th generation, instead of the probablity of extinction
eventually, which is the one we truely wanted...

But if I choose larger "R" and "G", the simulation time is huge. The above
R=G=100 already requires about several hours to run.

How to improve this simulation? Please throw me some lights. Thanks a lot!

-----------------------------------------------

%do 100 runs.
R=100;
%number of generations=100.
for kk=1:R;
    % do 100 generations.
    G=100;
    %using 1 for red cells and 0 for white cells.
    %generation 0 has 1 red cells.
    %only count the number of red cells.
    c(1)=1;
    for ii=2:G+1;
        cc=c(ii-1);
        for jj=1:c(ii-1);
            pp=rand(1);
            if pp<(1/4)
                cc=cc+1;
            elseif pp<(1/4+2/3)
                cc=cc;
            else
                cc=cc-1;
            end;
        end;
        c(ii)=cc;
    end;
    r(kk)=cc;
    % r(kk) holds the final number of red cells in each run after 100
generations...
end;

-------------------------------------

Results:

>> r

r =

  Columns 1 through 6

     2593271 677674 14366112 0 11180607 0

  Columns 7 through 12

           0 21666187 8045838 1349831 0 1323461

  Columns 13 through 18

     5991805 1233180 6837545 5374192 8542731 30977858

  Columns 19 through 24

     9941138 9686152 11542565 5264693 0 9034561

  Columns 25 through 30

     5810981 4311641 0 4301531 2394947 0

  Columns 31 through 36

     5341814 4326037 8229412 0 802612 358343

  Columns 37 through 42

     9254916 1965937 14193663 7576754 12425946 11298911

  Columns 43 through 48

      155572 25635151 0 0 210559 27023402

  Columns 49 through 54

     2723847 0 812472 13527307 2051866 1892091

  Columns 55 through 60

     1114763 1880878 4468172 0 5088758 0

  Columns 61 through 66

           0 11706819 0 11960996 0 1319624

  Columns 67 through 72

    14287190 1706401 11244218 3295325 5002472 3607422

  Columns 73 through 78

           0 2181729 6089518 8182148 0 6573386

  Columns 79 through 84

      976495 0 2611333 7735778 0 4880631

  Columns 85 through 90

    10126430 4570258 11645403 0 6071208 7322621

  Columns 91 through 96

           0 6907213 10800563 0 16968393 814406

  Columns 97 through 100

           0 18954344 3892277 0