Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I wanna simulate a Galton-Watson Tree to a maximum of n generation given a reproduction law P. I use Maple but I am unable to create the edges of the tree whenever there are more than two vertices in the tree. I think that the problem is that I am unable to store the number of children that each node has. Here is the code:

I use the law below to test, but I believe it should work with any discrete law on posint. P:=RandomVariable(EmpiricalDistribution(Array([0,1,2]))): Then the procedure to generate the tree is

GWT:=proc(n,P) #Renvoie un Arbre de Galton-Watson avec au plus n génération sous loi P;
local G,g,p,s,N,k;
if n=0 then
if p[1]=1 then s=0; else s:=((p[1])^n-1)/(p[1]-1);
if p[1]>0 then
G:=AddEdge(G,{seq({1, 2+k*s},k = 0 .. p[1]-1)});
end proc;

Just to compare, I was able to code the construction of a m-ary tree up to n generation, given n and m with a pretty similar code:

  Arbre:=proc(n,m) # m-ary arbre à n génération
local k,T,t,s,N;
if n=0 then 
if m=1 then s:=0; else s:=(m^n-1)/(m-1);
T:=AddEdge(T,{seq({1, 2+k*s},k = 0 .. m-1)});
end proc;

I was wondering if anyone could help with the coding? Maybe another way of storing the children. It may also be that it cannot be done with recursion.


Edit#1 : Pseudo code would look like:

Input: n:nonneg integer, P: probability law
if n=0 then return 1 vertice
create 1 vertice; i=0; while i<=n 
for each vertices in generation i generate a random number p from Law P and connect it 
with p other vertices (in generation i+1); i=i+1;
Return the graph
share|cite|improve this question
Perhaps you should write the algorithm in pseudo code, in case the problem is there and other people can help you too. – Patrick Da Silva May 15 '12 at 21:25
up vote 0 down vote accepted

Here's my version, which seems to work well in Maple 16.

uses Statistics, GraphTheory;
local v,v1,v2,gen,S,T1,T2,E,i;
  v:= 1; E:= {};
  v1:=1; v2:= 1;
  for gen from 1 to n do
    if v1 > v2 then break end if;
    S:= map(round, Sample(P,v2-v1+1));
    T2[v1-1]:= v2; 
    for i from v1 to v2 do
      T1[i]:= T2[i-1]+1;
      T2[i]:= T1[i]+S[i-v1+1]-1;
      E:= E union {seq({i,j},j=T1[i] .. T2[i])};
    end do;
    v1:= T1[v1];
    v2:= T2[v2];
  end do;
end proc;  

To try it out (with your P):

GraphTheory:-DrawGraph(GWT(5,P), style=tree, root=1);

enter image description here

share|cite|improve this answer
Thanks Robert! I'll look it up tomorrow and try to break down what you've done. I tested it as well and got it to work with Probability laws that take only a finite numbers of values, but I couldnt get it to process with, say a Poisson or a Geometric distribution. – Jean-Sébastien May 16 '12 at 3:45
Forget what I said about the laws, it works with everything I tried. Thanks again – Jean-Sébastien May 16 '12 at 3:49

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.