R: How do I approximate distributions of AICs with simulated data?

Question

I have likelihood functions LFi,i=1,2 with parameters [lambda,p] and [lambda]:

LF1<-lambda^sp^(n-n0)(1-p)^n0/((1-1/exp(lambda))^(n-n0)(exp(lambda))^(n-n0));
LF2<-lambda^s/(exp(lambda))^n

where s is the sum of observed counts, n is the number of observations, and n0 is the number of observed zeros. As usual, maximised log(LFi)s and AICs can be computed for each i. Now I'd like to

1) simulate a sequence of n Poisson(lambda_true) variates N times;

2) given lambda_true, count how often (out of N) the min AIC was due to LF1;

3) approximate the distribution of the two AICs (to compare them).

Anyone can help me learn how do I do that properly in R? Thanks.