transdimensional-EM

A trandimensional MCEM

set up observed tree

s=sim_phyl(seed=3)
s2 = s$tree.extant

initial parameter and 1st MC

pars=mle_tree(s$tree)
n_trees=1000
time =proc.time()
recs = sim_srt(wt=s2$wt,pars=pars,n_trees = n_trees)
proc.time() - time

##    user  system elapsed 
##       3       0       3

observe dimensions

size = vector(mode = 'numeric',length = n_trees)
for(i in 1:n_trees){
  rec = recs[[i]]
  size[i] = length(rec$wt)
}
qplot(size,geom='histogram')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

table(size)

## size
##  68  74  80  82  84  86  88  90  92  94  96  98 100 102 104 106 108 110 
##   1   2   3   3   8  13  19  22  30  41  53  57  75  65  79 102  82  72 
## 112 114 116 118 120 122 124 126 128 130 132 138 142 
##  60  59  37  33  21  22  13  10   6   3   5   3   1

length(s$tree$wt)

## [1] 102

generation the meta set of trees

dim <- function(list){
  length(list$wt)
}
checkdim <- function(list,i){
  dim(list) == i
}
meta_srt <- function(recs){
  size = vector(mode = 'numeric',length = length(recs))
  for(i in 1:n_trees){
    rec = recs[[i]]
    size[i] = length(rec$wt)
  }
  count = as.data.frame(table(size))
  Mrecs = list()
  for(i in count$size){
    Mrecs = c(Mrecs,list(tree=recs[lapply(recs,dim) == i]))
  }
  return(list(Mrecs=Mrecs,count=count))
}

Mrecs = meta_srt(recs)
#count = Mrecs$count
#Mrecs=Mrecs$Mrecs
#length(Mrecs[[1]])
#length(Mrecs[[18]])
#ahora, estimar parametros para cada grupo y ponerlos en su df con su dimension y peso...

estimations_matrix <- function(Mrecs){
  count = Mrecs$count
  Mrecs=Mrecs$Mrecs
  m=length(Mrecs)
  M = matrix(ncol=5,nrow=m)
  for(i in 1:m){
      p = subplex(par=c(2,1,80),fn=llik_st,setoftrees = Mrecs[[i]])
      M[i,1] = as.numeric(as.character(count$size[i]))
      M[i,2] = count$Freq[i]
      M[i,3:5] = p$par
  }
  return(M)
}

M = estimations_matrix(Mrecs)
sum(M[,1]*M[,3])/sum(M[,1])

## [1] 0.8297868

sum(M[,1]*M[,4])/sum(M[,1])

## [1] 0.08665098

sum(M[,1]*M[,5])/sum(M[,1])

## [1] 39.9701

subplex(par=c(2,1,80),fn=llik_st,setoftrees = recs)$par

## [1]  0.64903121  0.08382141 41.80198210

pars

## [1]  0.88886163  0.08410373 40.74476402

# Ok, we have the estimations and matrix for next iteration... 
m = nrow(M)
MM = list()
for(i in 1:m){
  n_trees = M[i,2]
  pars = M[i,3:5]
  MM = c(MM,sim_srt(wt=s2$wt,pars=pars,n_trees = n_trees))
}
n_trees = length(MM)
size = vector(mode = 'numeric',length = n_trees)
for(i in 1:n_trees){
  rec = MM[[i]]
  size[i] = length(rec$wt)
}
qplot(size,geom='histogram')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Mrecs = meta_srt(recs)
#count = Mrecs$count
#Mrecs=Mrecs$Mrecs
#m=length(Mrecs)
M = estimations_matrix(Mrecs)
M

##       [,1] [,2]      [,3]       [,4]     [,5]
##  [1,]   68    1 0.6522617 0.04720154 37.01687
##  [2,]   74    2 0.7349055 0.05024834 39.45220
##  [3,]   80    3 0.7654955 0.06102321 36.63922
##  [4,]   82    3 0.6698479 0.06104050 39.81561
##  [5,]   84    8 0.7961511 0.06364183 38.40640
##  [6,]   86   13 0.7564933 0.06524713 39.70007
##  [7,]   88   19 0.7234102 0.06668620 40.03981
##  [8,]   90   22 0.6647717 0.06933196 40.47570
##  [9,]   92   30 0.6842927 0.07123600 40.45750
## [10,]   94   41 0.6308719 0.07299010 41.59209
## [11,]   96   53 0.7046054 0.07498968 40.40724
## [12,]   98   57 0.7255015 0.07723244 40.42927
## [13,]  100   75 0.7389280 0.08013433 40.17690
## [14,]  102   65 0.7198494 0.08019416 40.96863
## [15,]  104   79 0.7343407 0.08264250 40.85585
## [16,]  106  102 0.6589164 0.08384837 41.90242
## [17,]  108   82 0.7601050 0.08524965 40.86969
## [18,]  110   72 0.8178037 0.08794716 40.21846
## [19,]  112   60 0.7075206 0.08982214 41.05115
## [20,]  114   59 0.8485831 0.09205627 40.20944
## [21,]  116   37 0.8868974 0.09214597 39.94341
## [22,]  118   33 0.8867810 0.09273478 40.36408
## [23,]  120   21 0.9189229 0.09473646 40.14669
## [24,]  122   22 0.9173494 0.09693036 40.10554
## [25,]  124   13 0.9469329 0.10232616 39.41625
## [26,]  126   10 1.0050369 0.10110139 40.10378
## [27,]  128    6 0.9773943 0.10198337 39.64559
## [28,]  130    3 0.9343398 0.10368236 40.07864
## [29,]  132    5 1.0228318 0.10427636 40.14036
## [30,]  138    3 1.1147411 0.11080832 39.32738
## [31,]  142    1 1.0066756 0.12356054 37.92075

sum(M[,1]*M[,3])/sum(M[,1])

## [1] 0.8297868

sum(M[,1]*M[,4])/sum(M[,1])

## [1] 0.08665098

sum(M[,1]*M[,5])/sum(M[,1])

## [1] 39.9701

Applied transdimensional MCEM (TMCEM)

s=sim_phyl(seed=3)
s2 = s$tree.extant
init_pars = c(2,1,80)
n_trees=1000
time = proc.time()
recs = sim_srt(wt=s2$wt,pars=init_pars,n_trees = n_trees)
proc.time()-time

##    user  system elapsed 
## 150.604   0.012 150.610

#for 

Mrecs = meta_srt(recs)
#count = Mrecs$count
#Mrecs=Mrecs$Mrecs
time=proc.time()
M = estimations_matrix(Mrecs)
proc.time() - time

##    user  system elapsed 
##  82.360   0.020  82.379

superlista = vector(mode='list',length = 40)
time2=proc.time()
for(l in 1:40){
time = proc.time()
  m = nrow(M)
  recs = list()
  for(i in 1:m){
    n_trees = M[i,2]
    pars = M[i,3:5]
    recs = c(recs,sim_srt(wt=s2$wt,pars=pars,n_trees = n_trees))
  }
  n_trees = length(recs)
  size = vector(mode = 'numeric',length = n_trees)
  for(i in 1:n_trees){
    rec = recs[[i]]
    size[i] = length(rec$wt)
  }
  print(qplot(size,geom='histogram'))
  Mrecs = meta_srt(recs)
 # count = Mrecs$count
 # Mrecs=Mrecs$Mrecs
  M = estimations_matrix(Mrecs)
  superlista[[l]] = M
  print(proc.time() - time)
  sum(M[,1]*M[,3])/sum(M[,1])
  sum(M[,1]*M[,4])/sum(M[,1])
  sum(M[,1]*M[,5])/sum(M[,1])
}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
## 197.420   0.032 197.459

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
## 172.616   0.016 172.637

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
## 150.016   0.004 150.050

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
## 127.464   0.000 127.465

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
## 114.644   0.000 114.655

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  92.200   0.000  92.207

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  80.512   0.004  80.545

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  68.428   0.000  68.452

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  57.516   0.000  57.516

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  51.688   0.000  51.683

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  45.436   0.000  45.438

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  40.984   0.000  40.985

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  38.164   0.000  38.190

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  37.108   0.000  37.146

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  32.556   0.100  32.670

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  29.592   0.004  29.596

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  28.576   0.152  28.753

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  26.248   0.000  26.251

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  24.776   0.000  24.778

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  24.308   0.000  24.328

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  25.252   0.004  25.278

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  23.400   0.000  23.408

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  22.064   0.000  22.078

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  21.984   0.000  22.019

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  21.300   0.004  21.327

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  20.004   0.000  20.003

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.580   0.000  19.585

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  18.956   0.000  18.957

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.304   0.000  19.307

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.872   0.000  19.887

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.316   0.004  19.334

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.192   0.228  19.419

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  18.052   0.000  18.053

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  18.556   0.000  18.562

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  17.988   0.004  17.994

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  18.604   0.004  18.609

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  18.244   0.000  18.247

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.252   0.000  19.268

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.744   0.012  19.763

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##    user  system elapsed 
##  19.788   0.000  19.790

proc.time()-time2

##     user   system  elapsed 
## 1840.744    0.572 1841.741

lambda=vector(mode='numeric',length=length(superlista))
mu = vector(mode='numeric',length=length(superlista))
for(i in 1:length(superlista)){
  M = superlista[[i]]
  lambda[i] = sum(M[,1]*M[,3])/sum(M[,1])
  mu[i] = sum(M[,1]*M[,4])/sum(M[,1])
  print(sum(M[,1]*M[,5])/sum(M[,1]))
}

## [1] 100.8514
## [1] 105.8426
## [1] 107.6146
## [1] 106.8492
## [1] 104.1293
## [1] 98.9889
## [1] 93.6636
## [1] 87.4043
## [1] 82.12475
## [1] 77.43614
## [1] 72.75213
## [1] 68.67908
## [1] 64.44687
## [1] 60.9607
## [1] 57.79565
## [1] 55.3054
## [1] 53.01928
## [1] 50.88001
## [1] 48.81282
## [1] 47.27042
## [1] 45.67818
## [1] 44.75304
## [1] 43.14697
## [1] 42.37498
## [1] 41.2965
## [1] 40.67438
## [1] 40.25462
## [1] 39.99555
## [1] 39.51122
## [1] 38.90683
## [1] 38.78959
## [1] 38.47266
## [1] 38.39976
## [1] 38.06462
## [1] 38.14314
## [1] 37.97935
## [1] 38.11994
## [1] 38.02256
## [1] 37.92704
## [1] 37.82183

qplot(lambda,mu)

TMCEM

rd = length(s$tree$wt)
rd

## [1] 102

sub_meta <- function(recs,rdim){
  size = vector(mode = 'numeric',length = length(recs))
  for(i in 1:n_trees){
    rec = recs[[i]]
    size[i] = length(rec$wt)
  }
  count = as.data.frame(table(size))
  Mrecs = list()
#  for(i in count$size){
    Mrecs = recs[lapply(recs,dim) == rdim]
#  }
  return(list(Mrecs=Mrecs,count=count))
}

pars=c(0.9,0.1,40)
n_it = 50
MMM = matrix(nrow=n_it,ncol=3)
n_trees=50000
for(j in 1:n_it){
  #print(pars)
  recs = sim_srt(wt=s2$wt,pars=pars,n_trees = n_trees)
  n_trees = length(recs)
  size = vector(mode = 'numeric',length = n_trees)
  for(i in 1:n_trees){
    rec = recs[[i]]
    size[i] = length(rec$wt)
  }
  print(qplot(size,geom='histogram'))
  sMrecs = sub_meta(recs,rd)
  sMrecs = sMrecs$Mrecs
  sMrecs
  p = subplex(par=c(2,1,80),fn=llik_st,setoftrees = sMrecs)
  pars = p$par
  MMM[j,] = pars
}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

MMM

##            [,1]       [,2]     [,3]
##  [1,] 0.6511724 0.07956196 42.13626
##  [2,] 0.5938594 0.08212105 42.22457
##  [3,] 0.6235101 0.08239771 41.65961
##  [4,] 0.6203222 0.08233729 41.63032
##  [5,] 0.6218696 0.08231023 41.64987
##  [6,] 0.5914099 0.08233795 42.17856
##  [7,] 0.5931867 0.08252984 42.18321
##  [8,] 0.5927653 0.08248465 42.18324
##  [9,] 0.6225664 0.08261316 41.63275
## [10,] 0.6217351 0.08229371 41.64714
## [11,] 0.6581938 0.08216057 41.13895
## [12,] 0.6207232 0.08213455 41.64950
## [13,] 0.6211107 0.08233681 41.64003
## [14,] 0.6223863 0.08234757 41.65016
## [15,] 0.6208616 0.08228887 41.63913
## [16,] 0.5917039 0.08232017 42.18341
## [17,] 0.6240415 0.08253844 41.65414
## [18,] 0.5912057 0.08226863 42.18292
## [19,] 0.5923098 0.08257776 42.16975
## [20,] 0.5924918 0.08254560 42.17445
## [21,] 0.6238504 0.08253695 41.65348
## [22,] 0.6208320 0.08228639 41.63897
## [23,] 0.5899934 0.08246440 42.15170
## [24,] 0.5679145 0.08250687 42.73767
## [25,] 0.6251565 0.08270215 41.65133
## [26,] 0.6229688 0.08224236 41.66511
## [27,] 0.6236454 0.08226130 41.66960
## [28,] 0.6207684 0.08235244 41.63272
## [29,] 0.6211857 0.08227319 41.64364
## [30,] 0.6221897 0.08223154 41.65714
## [31,] 0.6209567 0.08226390 41.64150
## [32,] 0.6213979 0.08240907 41.63530
## [33,] 0.6219205 0.08238689 41.64185
## [34,] 0.6221263 0.08236275 41.64582
## [35,] 0.5934054 0.08233496 42.20175
## [36,] 0.6242931 0.08240507 41.66814
## [37,] 0.5668064 0.08228021 42.74228
## [38,] 0.6269755 0.08266461 41.67394
## [39,] 0.6219755 0.08241330 41.64185
## [40,] 0.6210751 0.08238042 41.63471
## [41,] 0.5923540 0.08226017 42.19588
## [42,] 0.5929803 0.08271330 42.16569
## [43,] 0.6240340 0.08252229 41.65475
## [44,] 0.6203573 0.08241830 41.62441
## [45,] 0.5909301 0.08237742 42.16976
## [46,] 0.5921833 0.08258790 42.16745
## [47,] 0.6224892 0.08249566 41.64159
## [48,] 0.6567531 0.08217423 41.12306
## [49,] 0.6185624 0.08210098 41.62945
## [50,] 0.6202857 0.08230280 41.63162

mle_tree(s$tree)

## [1]  0.88886163  0.08410373 40.74476402

s=sim_phyl(mu0=0.4,seed=3)
rd = length(s$tree$wt)
rd

## [1] 287

s2 = s$tree.extant
pars=c(1.2,0.3,40)
n_it = 50
MMM = matrix(nrow=n_it,ncol=3)
n_trees=10000
for(j in 1:n_it){
  #print(pars)
  recs = sim_srt(wt=s2$wt,pars=pars,n_trees = n_trees)
  n_trees = length(recs)
  size = vector(mode = 'numeric',length = n_trees)
  for(i in 1:n_trees){
    rec = recs[[i]]
    size[i] = length(rec$wt)
  }
  print(qplot(size,geom='histogram'))
  sMrecs = sub_meta(recs,rd)
  sMrecs = sMrecs$Mrecs
  sMrecs
  p = subplex(par=c(2,1,80),fn=llik_st,setoftrees = sMrecs)
  pars = p$par
  MMM[j,] = pars
}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

MMM

##           [,1]      [,2]     [,3]
##  [1,] 1.087239 0.3143710 37.48299
##  [2,] 1.071940 0.3220699 36.07030
##  [3,] 1.090918 0.3246159 35.32217
##  [4,] 1.131392 0.3240413 34.78653
##  [5,] 1.176092 0.3226267 34.33702
##  [6,] 1.205662 0.3220935 34.02351
##  [7,] 1.211292 0.3243962 33.74556
##  [8,] 1.236106 0.3247029 33.54922
##  [9,] 1.232833 0.3268427 33.33283
## [10,] 1.251501 0.3258224 33.34709
## [11,] 1.247896 0.3263675 33.24202
## [12,] 1.271033 0.3268295 33.14423
## [13,] 1.280133 0.3253345 33.12600
## [14,] 1.300602 0.3246287 33.08191
## [15,] 1.296330 0.3265938 32.91103
## [16,] 1.319807 0.3279351 32.68053
## [17,] 1.313371 0.3286779 32.58483
## [18,] 1.334889 0.3265363 32.62886
## [19,] 1.349140 0.3255544 32.62208
## [20,] 1.362774 0.3241111 32.65958
## [21,] 1.376007 0.3231800 32.73204
## [22,] 1.392845 0.3244325 32.61224
## [23,] 1.414476 0.3245159 32.48426
## [24,] 1.425561 0.3236990 32.49314
## [25,] 1.414747 0.3235606 32.54351
## [26,] 1.415369 0.3222123 32.57895
## [27,] 1.408988 0.3232344 32.45904
## [28,] 1.398293 0.3245735 32.45220
## [29,] 1.381155 0.3279197 32.25754
## [30,] 1.384600 0.3276393 32.35612
## [31,] 1.391675 0.3270588 32.35386
## [32,] 1.397541 0.3236624 32.48764
## [33,] 1.396265 0.3243382 32.51724
## [34,] 1.384132 0.3251238 32.48880
## [35,] 1.379865 0.3272897 32.27485
## [36,] 1.395200 0.3256359 32.32352
## [37,] 1.376958 0.3250698 32.36121
## [38,] 1.397178 0.3230896 32.48139
## [39,] 1.400389 0.3241574 32.47425
## [40,] 1.390997 0.3261717 32.38811
## [41,] 1.385926 0.3266155 32.17641
## [42,] 1.381946 0.3270737 32.17309
## [43,] 1.389940 0.3257508 32.29418
## [44,] 1.255971 0.3258267 32.59900
## [45,] 1.287059 0.3273331 32.71031
## [46,] 1.262473 0.3278451 32.62121
## [47,] 1.291168 0.3287884 32.66239
## [48,] 1.306374 0.3272277 32.78325
## [49,] 1.314460 0.3259234 32.80293
## [50,] 1.332139 0.3235436 32.88893

mle_tree(s$tree)

## [1]  0.9148556  0.4672226 30.9270569