# Lognormal distribution and Galton Watson process

I have been studying a Galton Watson process that creates random binary trees with probability of survival $$p_{s}$$ and offspring distribution $$p(k)=p_{s}\delta(k-2)+(1-p_{s})\delta(k)$$. I'm particularly interested in the subcritical regime ($$p_{s}<0.5$$) where the tree eventually stops growing with probability 1. I'm not focusing on the standard variables, i.e., the number of elements at generation $$n$$, $$X_n$$, but on the total size of the tree $$T=\sum_n^{\inf}X_n$$. From the numerical simulations attached below, I find that a Lognormal distribution emerges as a good approximation of $$T$$ across several iterations of the same stochastic process, as $$X_0=a$$ increases, especially around $$p_{s}\approx 0.45$$. It is unclear to me why this happens since I would have a better argument for Lognormality in the supercritical regime ($$p_{s}<0.5$$) with noisy $$p_{s}$$, emerging due to the dominance of the last $$X_n$$ in $$T$$, and driven by the multiplicative central limit. I am very new to this literature, and I would appreciate any help or relevant references.

a=100;
c=2;

loop_population=100000;
pvec=linspace(0.1,0.5);
dp=length(pvec);
max_iter=100;
time2extinction=zeros(loop_population, dp);
random_var_vec=zeros(loop_population, dp);

for jnd=1:dp
selp=pvec(jnd);
for ind=1:loop_population
noutcome=a;
counter=1;
veciter=zeros(1, max_iter);
veciter(1)=noutcome;
while(noutcome>0)
noutcome=c*binornd(noutcome,selp);
counter=counter+1;
veciter(counter)=noutcome;
if(counter==max_iter)
break
end
end
time2extinction(ind,jnd)=counter;
random_var_vec(ind,jnd)=sum(veciter);

end
end

ksLvector=zeros(1,length(pvec));
ksWvector=zeros(1,length(pvec));
ksGvector=zeros(1,length(pvec));
ksNvector=zeros(1,length(pvec));

ksLvectorTE=zeros(1,length(pvec));
ksWvectorTE=zeros(1,length(pvec));
ksGvectorTE=zeros(1,length(pvec));
ksNvectorTE=zeros(1,length(pvec));

for ind=1:length(pvec)
vector2test=random_var_vec(:,ind);

[pHat,pCI] = lognfit(vector2test);
test_cdf = makedist('Lognormal','mu',pHat(1),'sigma',pHat(2));
[~,~,ksLvector(ind),~]= kstest(vector2test,'CDF',test_cdf);

[pHat,pCI]  = wblfit(vector2test);
test_cdf = makedist('Weibull','A',pHat(1),'B',pHat(2));
[~,~,ksWvector(ind),~]= kstest(vector2test,'CDF',test_cdf);

[pHat,pCI] = gamfit(vector2test);
test_cdf = makedist('Gamma','a',pHat(1),'b',pHat(2));
[~,~,ksGvector(ind),~]= kstest(vector2test,'CDF',test_cdf);

[muHat,sigmaHat]  = normfit(vector2test);
test_cdf = makedist('Normal','mu',muHat,'sigma',sigmaHat);
[~,~,ksNvector(ind),~]= kstest(vector2test,'CDF',test_cdf);

vector2test=time2extinction(:,ind);

[pHat,pCI] = lognfit(vector2test);
test_cdf = makedist('Lognormal','mu',pHat(1),'sigma',pHat(2));
[~,~,ksLvectorTE(ind),~]= kstest(vector2test,'CDF',test_cdf);

[pHat,pCI]  = wblfit(vector2test);
test_cdf = makedist('Weibull','A',pHat(1),'B',pHat(2));
[~,~,ksWvectorTE(ind),~]= kstest(vector2test,'CDF',test_cdf);

[pHat,pCI] = gamfit(vector2test);
test_cdf = makedist('Gamma','a',pHat(1),'b',pHat(2));
[~,~,ksGvectorTE(ind),~]= kstest(vector2test,'CDF',test_cdf);

[muHat,sigmaHat]  = normfit(vector2test);
test_cdf = makedist('Normal','mu',muHat,'sigma',sigmaHat);
[~,~,ksNvectorTE(ind),~]= kstest(vector2test,'CDF',test_cdf);
end

figure,
loglog(pvec, mean(random_var_vec))
hold on
loglog(pvec, std(random_var_vec))
xlabel('p')
legend('mean tree size', 'std tree size')

figure,
plot(pvec, ksLvector)
hold on
plot(pvec, ksWvector)
hold on
plot(pvec, ksGvector)
hold on
plot(pvec, ksNvector)
legend('Lognormal', 'Weibull', 'Gamma', 'Normal')
xlabel('p')
ylabel('kstat')
title (['tree size a=', num2str(a)])

figure,
plot(pvec, ksLvectorTE)
hold on
plot(pvec, ksWvectorTE)
hold on
plot(pvec, ksGvectorTE)
hold on
plot(pvec, ksNvectorTE)
legend('Lognormal', 'Weibull', 'Gamma', 'Normal')
xlabel('p')
ylabel('kstat')
title (['time to extinction a=', num2str(a)])

figure,
histogram(log(random_var_vec(:,88)))
title(['p=', num2str(pvec(88))])
xlabel('log(T)')
ylabel('PDF')

figure,
histogram(log(time2extinction(:,88)), 20)
title(['p=', num2str(pvec(88))])
xlabel('log(time to extinction)')
ylabel('PDF')