it=100000
llik = vector(mode='numeric',length = it)
Nm = vector(mode='numeric',length = it)
size = vector(mode='numeric',length = it)
extant = vector(mode='numeric',length = it)
for(i in 1:it){
  s = sim_phyl()
  llik[i] = -llik(pars=c(0.8,0.1,40),n = s$tree$n,E = s$tree$E,t = s$tree$wt)
  Nm[i] = number_missing(s)
  size[i] = length(s$tree$wt)  
  extant[i] = size[i] - Nm[i]
}


qplot(size,llik)

qplot(Nm,llik)

qplot(extant,llik)

a=lm(llik~size)
summary(a)
## 
## Call:
## lm(formula = llik ~ size)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.435  -7.241  -0.505   6.757  46.421 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.73628    0.21204    69.5   <2e-16 ***
## size        -2.73431    0.00187 -1461.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.32 on 99998 degrees of freedom
## Multiple R-squared:  0.9553, Adjusted R-squared:  0.9553 
## F-statistic: 2.137e+06 on 1 and 99998 DF,  p-value: < 2.2e-16
a$coefficients[1]
## (Intercept) 
##    14.73628
a$coefficients[2]
##      size 
## -2.734306
it=100000
llik = vector(mode='numeric',length = it)
Nm = vector(mode='numeric',length = it)
size = vector(mode='numeric',length = it)
extant = vector(mode='numeric',length = it)
for(i in 1:it){
  s = sim_phyl()
  llik[i] = -llik(pars=c(1.8,0.4,40),n = s$tree$n,E = s$tree$E,t = s$tree$wt)
  Nm[i] = number_missing(s)
  size[i] = length(s$tree$wt)  
  extant[i] = size[i] - Nm[i]
}

qplot(size,llik)

qplot(Nm,llik)

qplot(extant,llik)

a=lm(llik~size)
summary(a)
## 
## Call:
## lm(formula = llik ~ size)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -123.056  -21.730   -0.921   20.875  131.071 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -46.433814   0.644939   -72.0   <2e-16 ***
## size         -3.295959   0.005686  -579.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.34 on 99998 degrees of freedom
## Multiple R-squared:  0.7706, Adjusted R-squared:  0.7706 
## F-statistic: 3.36e+05 on 1 and 99998 DF,  p-value: < 2.2e-16
a$coefficients[1]
## (Intercept) 
##   -46.43381
a$coefficients[2]
##      size 
## -3.295959

dimension in differenth simulations

time = proc.time()
it = 100000
size2 = vector(mode='numeric',length = it)
for(i in 1:it){
  s = sim_phyl()
  size2[i] = length(s$tree$wt)
}
print(proc.time()-time)
##      user    system   elapsed 
## 17543.392     0.224 17542.347
s=sim_phyl(seed=1)
pars=mle_tree(s)
length(s$tree$wt)
## [1] 120
size3 = vector(mode='numeric',length = it)
time=proc.time()
for(i in 1:it){
  rec=rec_tree(wt=s$tree.extant$wt,pars=pars)
  size3[i] = length(rec$tree$wt)
}
print(proc.time()-time)
##    user  system elapsed 
## 227.376   0.000 227.352
time=proc.time()
for(i in 1:it){
  rec=rec_tree(wt=s$tree.extant$wt,pars=pars)
  size3[i] = length(rec$wt)
}
print(proc.time()-time)
##    user  system elapsed 
## 252.656   0.004 252.641
length(s$tree$wt)
## [1] 120
summary(size2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     3.0   102.0   113.0   111.9   124.0   177.0
qplot(size2,geom='histogram')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(size3)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      80     108     116     116     122     168
qplot(size3,geom='histogram')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.