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`.
