n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl()
  pars = c(8,1,70)
  p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
}
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4548  0.7241  0.8248  0.8261  0.9044  1.2450
qplot(p[,2], geom="histogram") 

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.05803 0.09153 0.10210 0.10280 0.11210 0.15170
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 3.600e+01 3.900e+01 4.000e+01 1.268e+12 4.100e+01 1.268e+14
length(p[p[,3]<80,3])
## [1] 99

and if I start closer to the solution?

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl()
  pars = c(0.9,0.2,50)
  p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
}
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5430  0.7597  0.8330  0.8526  0.9688  1.2770
qplot(p[,2], geom="histogram") 

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.06448 0.09253 0.10170 0.10290 0.11140 0.15390
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   34.39   39.11   39.97   39.83   40.50   51.37

now with different parameters

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.2)
  pars = pars = c(8,1,70)
  p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
}
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4877  0.7659  0.8490  0.8642  0.9379  1.4440
qplot(p[,2], geom="histogram") 

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1524  0.1831  0.2024  0.2021  0.2186  0.2555
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 3.000e+01 3.800e+01 4.000e+01 3.633e+12 4.200e+01 3.633e+14
n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.4)
  pars = pars = c(8,1,70)
  p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
}
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5616  0.8014  0.8886  0.8978  0.9756  1.4760
qplot(p[,2], geom="histogram") 

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3136  0.3684  0.3931  0.4017  0.4321  0.5491
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 7.000e+00 3.500e+01 3.900e+01 4.050e+12 4.500e+01 2.851e+14
n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0)
  pars = pars = c(8,1,70)
  p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
}
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4797  0.6973  0.7894  0.8085  0.9273  1.2790
qplot(p[,2], geom="histogram") 

summary(p[,2])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.092e-18 8.668e-18 1.996e-17 2.441e-17 3.245e-17 1.079e-16
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   38.57   39.70   40.11   40.94   41.61   53.49

Crown age = 10. Note that we have to drop mall trees now

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
for(i in 1:n_sim){
  s = sim_phyl(ct = 10)
  if(s$newick$Nnode>5){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    i = i-1
    k=k+1}
}
k
## [1] 1
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2891  0.7026  0.8489  0.8187  0.9452  1.3160
qplot(p[,2], geom="histogram") 

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02791 0.08812 0.09947 0.10170 0.11410 0.20330
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 2.100e+01 3.800e+01 4.100e+01 2.783e+13 4.300e+01 9.057e+14
length(p[p[,3]<80,3])
## [1] 90

now with different parameters

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k= 0
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.2, ct = 10)
  if(s$newick$Nnode>5){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    i = i-1
    k=k+1}
}
k
## [1] 1
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1836  0.7647  0.8930  0.8771  0.9977  1.3360
qplot(p[,2], geom="histogram") 

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.1755  0.1995  0.2011  0.2342  0.3061
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 7.000e+00 3.500e+01 4.000e+01 3.087e+12 4.300e+01 1.056e+14
n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.4, ct = 10)
   if(s$newick$Nnode>15){
    pars = c(8,1,70)
    p[i-k,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    k=k+1}
}
k
## [1] 7
qplot(p[,1], geom="histogram") 
## Warning: Removed 7 rows containing non-finite values (stat_bin).

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.6291  0.9186  1.0160  1.0230  1.1410  1.8120       7
qplot(p[,2], geom="histogram") 
## Warning: Removed 7 rows containing non-finite values (stat_bin).

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.2911  0.3589  0.3899  0.3983  0.4241  0.6295       7
qplot(p[p[,3]<80,3], geom="histogram") 
## Warning: Removed 7 rows containing non-finite values (stat_bin).

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## 5.000e+00 2.300e+01 3.000e+01 1.347e+12 4.200e+01 9.584e+13         7
n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0,ct = 10)
   if(s$newick$Nnode>5){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    i = i-1
    k=k+1}
}
k
## [1] 0
qplot(p[,1], geom="histogram") 

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4239  0.7077  0.7982  0.8222  0.9371  1.3560
qplot(p[,2], geom="histogram") 

summary(p[,2])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 3.036e-20 6.801e-18 2.369e-17 2.638e-17 3.842e-17 1.075e-16
qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 3.000e+01 3.900e+01 4.000e+01 2.911e+10 4.200e+01 2.911e+12

Crown age = 4. Note that we have to drop mall trees now

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
for(i in 1:n_sim){
  s = sim_phyl(ct = 5)
  if(s$newick$Nnode>5){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    i = i-1
    k=k+1}
}
k
## [1] 14
qplot(p[,1], geom="histogram") 
## Warning: Removed 14 rows containing non-finite values (stat_bin).

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.4243  0.7945  0.9750  0.9876  1.1600  1.5310      14
qplot(p[,2], geom="histogram") 
## Warning: Removed 14 rows containing non-finite values (stat_bin).

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00000 0.06813 0.09988 0.10620 0.14320 0.33190      14
qplot(p[p[,3]<80,3], geom="histogram") 
## Warning: Removed 14 rows containing non-finite values (stat_bin).

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## 8.000e+00 2.300e+01 3.200e+01 3.849e+13 1.080e+02 6.001e+14        14
length(p[p[,3]<80,3])
## [1] 76

now with different parameters

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k= 0
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.2, ct = 5)
  if(s$newick$Nnode>5){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    i = i-1
    k=k+1}
}
k
## [1] 18
qplot(p[,1], geom="histogram") 
## Warning: Removed 18 rows containing non-finite values (stat_bin).

summary(p[,1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.5818  0.9055  1.1290  1.2170  1.4420  3.2360      18
qplot(p[,2], geom="histogram") 
## Warning: Removed 18 rows containing non-finite values (stat_bin).

summary(p[,2])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.1588  0.2163  0.2268  0.2629  0.6343      18
qplot(p[p[,3]<80,3], geom="histogram") 
## Warning: Removed 18 rows containing non-finite values (stat_bin).

summary(p[,3])
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## 3.000e+00 1.100e+01 2.300e+01 6.084e+13 4.200e+01 1.170e+15        18

this one has many issues normally (mu=0.4) and with small crown age is too difficult

n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
while(i < n_sim){
  s = sim_phyl(mu0 = 0.4, ct = 5,seed=i)
   if(s$newick$Nnode>3){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
    i=i+1
  }
  else{
    k=k+1}
}
k
n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
g=0
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.4, ct = 5,seed=i)
   if(s$newick$Nnode>5){
    g = g+1
  }
  else{
    k=k+1}
}
k
g
qplot(p[,1], geom="histogram") 
summary(p[,1])
qplot(p[,2], geom="histogram") 
summary(p[,2])
qplot(p[p[,3]<80,3], geom="histogram") 
summary(p[,3])
n_sim = 100
p=matrix(ncol=3,nrow=n_sim)
k=0
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0,ct = 5)
   if(s$newick$Nnode>5){
    pars = c(8,1,70)
    p[i,] = subplex(par = pars, fn = llik, n = s$n, E = s$E, t = s$wt)$par
  }
  else{
    i = i-1
    k=k+1}
}
k
qplot(p[,1], geom="histogram") 
summary(p[,1])
qplot(p[,2], geom="histogram") 
summary(p[,2])
qplot(p[p[,3]<80,3], geom="histogram") 
summary(p[,3])