Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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
G:=Graph(1);
else
g:=GWT(n-1,P);
p:=Sample(P,1);
G:=DisjointUnion(Graph(1),`$`(g,p[1]));
if p[1]=1 then s=0; else s:=((p[1])^n-1)/(p[1]-1);
fi;
N:=nops(Vertices(G));
G:=RelabelVertices(G,[seq(k,k=1..N)]);
if p[1]>0 then
G:=AddEdge(G,{seq({1, 2+k*s},k = 0 .. p[1]-1)});
fi;
G;
fi;
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 
T:=Graph(1);
else
t:=Arbre(n-1,m);
T:=DisjointUnion(Graph(1),`$`(t,m));
if m=1 then s:=0; else s:=(m^n-1)/(m-1);
fi;
N:=nops(Vertices(T));
T:=RelabelVertices(T,[seq(k,k=1..N)]);
T:=AddEdge(T,{seq({1, 2+k*s},k = 0 .. m-1)});
T;
fi;
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.

Thanks

Edit#1 : Pseudo code would look like:

Input: n:nonneg integer, P: probability law
if n=0 then return 1 vertice
else
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|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
add comment

1 Answer

up vote 0 down vote accepted

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

GWT:=proc(n,P) 
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;
  Graph([$1..v2],E);
end proc;  

To try it out (with your P):

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

enter image description here

share|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
add comment

Your Answer

 
discard

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.