The complete phylogenetic tree

n_sim = 10000
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.3408  0.7474  0.8378  0.8470  0.9379  1.6640

qplot(p[,2], geom="histogram") 

summary(p[,2])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.03678 0.08971 0.10100 0.10170 0.11260 0.29200

qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 5.000e+00 3.900e+01 4.000e+01 7.454e+11 4.100e+01 7.538e+14

length(p[p[,3]<60,3])

## [1] 9936

if we do it for 100 sumulations only? (like the paper)

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.3621  0.7390  0.8554  0.8415  0.9471  1.2130

qplot(p[,2], geom="histogram") 

summary(p[,2])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.06845 0.09033 0.10160 0.10280 0.11290 0.14280

qplot(p[p[,3]<80,3], geom="histogram") 

summary(p[,3])

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 3.500e+01 3.900e+01 4.000e+01 6.379e+11 4.100e+01 5.103e+13

length(p[p[,3]<60,3])

## [1] 98

and if I start closer to the solution?

n_sim = 10000
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.4110  0.7565  0.8483  0.8584  0.9500  2.1570

qplot(p[,2], geom="histogram") 

summary(p[,2])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1128  0.1847  0.2009  0.2022  0.2177  0.4522

qplot(p[p[,3]<60,3], geom="histogram") 

summary(p[p[,3]<60,3])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.539  37.870  39.860  39.710  41.630  59.850

now with different parameters

n_sim = 10000
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.2)
  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.4334  0.7986  0.8926  0.9048  0.9938  2.3010

qplot(p[,2], geom="histogram") 

summary(p[,2])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2178  0.3751  0.4007  0.4045  0.4287  0.6999

qplot(p[p[,3]<60,3], geom="histogram") 

summary(p[p[,3]<60,3])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.827  31.910  37.600  36.420  42.590  59.990

n_sim = 10000
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0.4)
  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.2415  0.6972  0.7867  0.8011  0.8906  1.6400

qplot(p[,2], geom="histogram") 

summary(p[,2])

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.000e+00 9.000e-18 1.920e-17 1.056e-16 3.440e-17 8.169e-13

qplot(p[p[,3]<60,3], geom="histogram") 

summary(p[p[,3]<60,3])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   37.09   39.68   40.24   40.88   41.31   59.91

n_sim = 10000
p=matrix(ncol=3,nrow=n_sim)
for(i in 1:n_sim){
  s = sim_phyl(mu0 = 0)
  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.2415  0.6972  0.7867  0.8011  0.8906  1.6400

qplot(p[,2], geom="histogram") 

summary(p[,2])

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.000e+00 9.000e-18 1.920e-17 1.056e-16 3.440e-17 8.169e-13

qplot(p[p[,3]<60,3], geom="histogram") 

summary(p[p[,3]<60,3])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   37.09   39.68   40.24   40.88   41.31   59.91