Recall the multiplier method assumes we have some known endpoints along with estimates of \(p,q\), and uses this information to create a variance weighted estimate of \(N\). Repeated sampling produces a range of estimates which can be used to plot a distribution as above. Let’s start with uniform \(p, q\). To reduce computation time, we return to smaller values for \(A, B, a,b\).

To use MM code, we rewrite the data to be read as a tree by the code which produces weighted MM estimates:

treedata <- read.table("MM_subtree.txt", sep = ",", header = TRUE)
treedata %>%
  knitr::kable() %>%
  kable_styling()
From To Est Total TrueEst Pop Desc Prob ProbLabel
N A 0 0 50 TRUE A 0.5 a
N AA 0 0 NA FALSE AA 0.5 aa
AA B 0 0 150 TRUE B 0.5 b
AA BB 0 0 NA FALSE BB 0.5 bb

Now we set the labels required:

treedata$Estimate <- treedata$Est
treedata$Population <- treedata$Pop
# display tree with true numbers
tree <- FromDataFrameNetwork(treedata)
SetEdgeStyle(tree, penwidth=1, 
             label = function(node) paste(node$ProbLabel, node$Prob, sep = ":\n"))
SetNodeStyle(tree, fontsize=12, penwidth=2, width=1,
             label = function(node) paste(node$Desc, node$TrueEst, sep = ":\n"))
plot(tree)

We now add the weights equation to the calculation of the root estimate, to obtain one weighted estimate given the estimates and variances of leaf paths.

## [1] "using weighted mean estimates"
##    levelName        uncertainty
## 1 N           3222 (-Inf - Inf)
## 2  ¦--A        1090 (51 - 2470)
## 3  °--AA          NaN (NA - NA)
## 4      ¦--B  7140 (188 - 43700)
## 5      °--BB      NaN (NA - NA)
## [1] "A : 0.646887954667278" "B : 0.353112045332722"

The estimates oF \(N\) given by each leaf is printed within the leaf, with a corresponding range. The weighted estimates of \(N\) are contained in \(N\), with weights determined as printed above. Negative lower bounds on the variance are possible since variance is calculated using \(\pm 2*\sqrt{Var}\), where \(Var = \frac{1}{e^T\Sigma^{-1} e}\), where \(\Sigma\) is the covariance matrix for the leaves used in estimation. We cannot plot a distribution of estimates of \(N\) with one pass, since mean values are used with weights to generate one estimate of \(N\) on each application of the multiplier method and replicates are in the values of estimates using each path seperately. However, we can repeat the method many times, generating multiple estimates of \(N\) to produce a similar plot. This takes some time, so I ran the code, saved the output, and then plotted, but the commented code is included for completeness:

# initialize vector to store estimates of N
# MMests <- vector("numeric",length = 5000)
# for (i in 1:length(MMests)){
#   weightedPopTree(tree,sample_length=1000,method='weightedEstimate')
#   tree$Do(function(node) {
#     if(node$isRoot){
#       node$uncertainty <-pop.conf.int(tree)
#       }
#     else{
#       node$uncertainty <- prettyConfidenceIntervals(node$overdose_samples)
#       }}
#     )
#   MMests[i] <- tree$root$estimate
# }

This method does not bound possible estimates of \(N\). Though unlikely, some values estimated can be very high. To help visualize, we cut the top 5% of values to make the distribution more comparable when plotting. Also, since weights can be negative, we also ignore any negative estimates of N that could theoretically exist as these values are non-sensical.

#MMests <- round(MMests)
#MMests <- as.data.frame(MMests)
#top_95 <- sort(MMests$MMests)[round(length(MMests$MMests) * 0.95)]
#MMests95 <- MMests[MMests$MMests <= top_95, ,drop = FALSE]
#MMests95 <- MMests95[MMests95$MMests >=0, , drop = FALSE]
#write.csv(MMests95, file = "uniformMMests95.csv", row.names = FALSE)
MMests95 <- read.csv2("uniformMMests95.csv")

p <- ggplot(MMests95, aes(x = MMests)) + geom_density(alpha=.2, fill="#FF6666") + ggtitle("Distribution of weighted MM estimates for uniform f(p), f(q)", subtitle = "5000 replicates")
p

The mode of the distribtuion of values estimated using the multiplier method for uniform \(f(p),f(q)\) does not coincide well to that suggested by the distribution \(f(N|A,B)\). To help compare these distributions, we can restrict our plot to the range \(N = (10,750)\), which was the upper and lower endpoints of the discrete uniform distribution used in \(f(N|A,B)\).

#MMests_restrict <- MMests[MMests$MMests >=10, , drop = FALSE]
#MMests_restrict <- MMests_restrict[MMests_restrict$MMests <=750, , drop = FALSE]
#write.csv(MMests_restrict, file = "uniformMMests.csv", row.names = FALSE)
MMests_restrict <- read.csv2("uniformMMests.csv")

p <- ggplot(MMests_restrict, aes(x = MMests)) + geom_density(alpha=.2, fill="#FF6666") + ggtitle("Distribution of weighted MM estimates for uniform f(p), f(q)", subtitle = "5000 replicates")
p

#ggsave("uniformMMests.pdf")

We can see this distribution is quite different from that of \(f(N|A,B)\) when \(f(p), f(q)\) are uniform. The lower bound given by the data (200) is not respected using this method, and the mode falls closer to 500 (and appears nearly bimodal, with a second peak at 600).

It should be noted, the data itself does not make sense with uniform binary splitting - while a population at node \(B\) suggests \(N\) should be around 600 if the population split equality at each level, estimating from \(A\) suggests \(N\) should be around 100. We are looking at how each method resolves this discrepancy. The estimate from \(B\) is reflected in this distribution with much more weight than an estimate from \(A\), suggesting \(B\) is influencing the estimate of \(N\) through the weighted average, despite the larger variance of estimating two unknown branches.

We now move on to the variations of plots, the first being where \(f(q)\sim Beta(150,1)\). To do this in each case in the MM code, we overwrite the data in the \(Estimate\) and \(Total\) columns, so that \(Beta(Estimate+1, Total-Estimate+1)\) gives the desired distribution.

peakedq1 <- treedata
peakedq1[3,3] <- 149
peakedq1[3,4] <- 149
peakedq1

Now we plot the distribution of \(N\) estimates for peaked \(f(q)\), with point mass on 1:

This plot shows more certainty towards the lower mode than the uniform plot, but still does not deviate entirely to the lower end of our range, as seen in the distribution of \(f(N|A,B)\) given this choice of \(f(q)\).

Now we plot the distribution of \(N\) estimates for peaked \(f(q)\) with point mass on 0:

We see here that for this choice of \(f(q)\) the distribution becomes much more uncertain over the entire range, but skews much of the weight towards the upper end of the range, as it did for this choice of \(f(q)\) in \(f(N|A,B)\).

New Method

We now try something a little different. For each pass of the multiplier method, 1000 values of each of \(p,q\) are generated, resulting in a covariance matrix which generates the weights. Given \(A,B\), an estimate of \(N\) is generated using a weighted sum of the average \(p,q\). Repeating this process 5000 times may not be the best comparison to distribution in notebook #1 or the Bayesian model in notebook #3. Therefore, we take another approach to compare to these methods. We generate 5000 values of \(p,q\) and the corresponding weights as in a singular pass of the weighted MM. However, instead of using average \(p,q\) values to compute one estimate, \(\hat{N}\), we use the paired \(p,q\) values from each sampling in the weighted average, resulting in 5000 estimates of \(N\). These are plotted as an alternative to the method above, which should intuitively result in a narrower distribution around the mode with increased compute-power, suggesting it may not be a desirable comparison to the distribution of notebooks #1 and #3.

set.seed(1)

# import data
uniformdata <- treedata
uniformdata$Estimate <- uniformdata$Est
uniformdata$Population <- uniformdata$Pop

# create tree with uniform data
tree <- FromDataFrameNetwork(uniformdata)

# calculate estimates from A,B using 5000 simulated p,q values
weightedPopTree(tree,sample_length=5000,method='weightedEstimate')
## [1] "using weighted mean estimates"
tree$Do(function(node) {
  if(node$isRoot){
    node$uncertainty <-pop.conf.int(tree)
    }
  else{
    node$uncertainty <- prettyConfidenceIntervals(node$overdose_samples)
    }}
  )

# create dataframe of N estimates from each leaf
x <- tree$Get('overdose_samples', filterFun = isLeaf, traversal = 'post-order')
# Ignore NAs and then also choose leaves with population=TRUE values only (these are the data)
if(is.list(x)){    
  mat.x <- NULL
  for (i in 1:length(x)) {
    if (length(x[[i]] > 0)) {
      mat.x <- cbind(mat.x, x[[i]])
      colnames(mat.x)[dim(mat.x)[2]] <- rownames(x)[[i]]
    }
    }
  x <- mat.x
  }
x <- as.data.frame(x)
  
# rename column to be of the correct leaves (with population=TRUE values)
g <- which(tree$Get('Population', filterFun = isLeaf, traversal = 'post-order'))
pop.x <- x %>% 
  select(g)

pop.x <- round(pop.x)
uniform.pop.x <- pop.x

# generate values \hat{N} to plot
weights <- as.matrix(pop.ko.weights(tree))
uniformweights <- weights
uniformNs <- as.data.frame(as.matrix(pop.x) %*% t(weights))
uniformNs <- round(uniformNs)
uniformNs.restrict <- uniformNs[uniformNs$V1>=10 & uniformNs$V1<=750, , drop = FALSE]
p <- ggplot(uniformNs.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for uniform distribution for q, uniform for p", subtitle = "One MM pass and set of weights, 5000 p,q estimates")
p

#ggsave("altuniformMMests.pdf")
peakedq1 <- treedata
peakedq1[3,3] <- 149
peakedq1[3,4] <- 149
peakedq1$Estimate <- peakedq1$Est
peakedq1$Population <- peakedq1$Pop
peakedq1
# create tree with peaked data
tree <- FromDataFrameNetwork(peakedq1)

# calculate estimates from A,B using 5000 simulated p,q values
weightedPopTree(tree,sample_length=5000,method='weightedEstimate')
## [1] "using weighted mean estimates"
tree$Do(function(node) {
  if(node$isRoot){
    node$uncertainty <-pop.conf.int(tree)
    }
  else{
    node$uncertainty <- prettyConfidenceIntervals(node$overdose_samples)
    }}
  )

# create dataframe of N estimates from each leaf
x <- tree$Get('overdose_samples', filterFun = isLeaf, traversal = 'post-order')

# Ignore NAs and then also choose leaves with population=TRUE values only (these are the data)
if(is.list(x)){    
  mat.x <- NULL
  for (i in 1:length(x)) {
    if (length(x[[i]] > 0)) {
      mat.x <- cbind(mat.x, x[[i]])
      colnames(mat.x)[dim(mat.x)[2]] <- rownames(x)[[i]]
    }
    }
  x <- mat.x
  }
x <- as.data.frame(x)
  
# rename column to be of the correct leaves (with population=TRUE values)
g <- which(tree$Get('Population', filterFun = isLeaf, traversal = 'post-order'))
pop.x <- x %>% 
  select(g)

pop.x <- round(pop.x)
peakedq1.pop.x <- pop.x


# generate values \hat{N} to plot
weights <- as.matrix(pop.ko.weights(tree))
peakedq1weights <- weights
peakedq1Ns <- as.data.frame(as.matrix(pop.x) %*% t(weights))
peakedq1Ns <- round(peakedq1Ns)
peakedq1Ns.restrict <- peakedq1Ns[peakedq1Ns$V1>=10 & peakedq1Ns$V1<=750, , drop = FALSE] # restrict to the range from a=10 to b=750
p <- ggplot(peakedq1Ns.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for peaked distribution for q, uniform for p", subtitle = "Point mass on 1. One MM pass and set of weights, 5000 p,q estimates.")
p

#ggsave("altpeakedq1MMests.pdf")

We see that this doesn’t seem to have much of an effect on the mode, even though the estimates from node \(B\) are significantly higher in the uniform case. There is however, more confidence that the estimate of \(N\) is smaller. This is because in the peaked case, the majority of the mass of \(f(q)\) sits on one, so we have higher confidence that the population at the parent of \(B\) is not likely to be significantly higher than that at \(B\):

print("Uniform estimates from A,B:")
## [1] "Uniform estimates from A,B:"
summary(uniform.pop.x)
##        A                  B          
##  Min.   :    50.0   Min.   :    153  
##  1st Qu.:    66.0   1st Qu.:    391  
##  Median :    98.0   Median :    778  
##  Mean   :   393.3   Mean   :  11320  
##  3rd Qu.:   198.0   3rd Qu.:   2316  
##  Max.   :114781.0   Max.   :9319242
print("Peaked q estimates from A,B:")
## [1] "Peaked q estimates from A,B:"
summary(peakedq1.pop.x)
##        A                  B            
##  Min.   :    50.0   Min.   :    150.0  
##  1st Qu.:    66.0   1st Qu.:    200.0  
##  Median :   100.0   Median :    303.0  
##  Mean   :   415.4   Mean   :   1605.6  
##  3rd Qu.:   197.2   3rd Qu.:    595.2  
##  Max.   :104249.0   Max.   :2193652.0

Once we examine the weighted estimates of \(N\), however, we can see the estimates from \(B\) have been downweighted in the uniform case, resulting in a similar distribution of estimates \(\hat{N}\) because of increased confidence in the estimate of \(B\) in the peaked case:

print("Uniform weighted estimates of N:")
## [1] "Uniform weighted estimates of N:"
summary(uniformNs)
##        V1           
##  Min.   :     88.0  
##  1st Qu.:    221.8  
##  Median :    416.0  
##  Mean   :   3975.7  
##  3rd Qu.:   1029.2  
##  Max.   :3055623.0
print("Weights in uniform case:")
## [1] "Weights in uniform case:"
uniformweights
##                         
## [1,] 0.6721284 0.3278716
print("Peaked q weighted estimates of N:")
## [1] "Peaked q weighted estimates of N:"
summary(peakedq1Ns)
##        V1           
##  Min.   :    101.0  
##  1st Qu.:    162.0  
##  Median :    245.0  
##  Mean   :   1011.4  
##  3rd Qu.:    474.2  
##  Max.   :1098401.0
print("Weights in peaked f(q) case:")
## [1] "Weights in peaked f(q) case:"
peakedq1weights
##                         
## [1,] 0.4992985 0.5007015

We move the point mass in the peak distribution to observe the effects:

peakedq0 <- treedata
peakedq0[3,3] <- 0  #this may need to be changed?
peakedq0[3,4] <- 149
peakedq0$Estimate <- peakedq0$Est
peakedq0$Population <- peakedq0$Pop
peakedq0
# create tree with peaked data
tree <- FromDataFrameNetwork(peakedq0)

# calculate estimates from A,B using 5000 simulated p,q values
weightedPopTree(tree,sample_length=5000,method='weightedEstimate')
## [1] "using weighted mean estimates"
tree$Do(function(node) {
  if(node$isRoot){
    node$uncertainty <-pop.conf.int(tree)
    }
  else{
    node$uncertainty <- prettyConfidenceIntervals(node$overdose_samples)
    }}
  )

# create dataframe of N estimates from each leaf
x <- tree$Get('overdose_samples', filterFun = isLeaf, traversal = 'post-order')

# Ignore NAs and then also choose leaves with population=TRUE values only (these are the data)
if(is.list(x)){    
  mat.x <- NULL
  for (i in 1:length(x)) {
    if (length(x[[i]] > 0)) {
      mat.x <- cbind(mat.x, x[[i]])
      colnames(mat.x)[dim(mat.x)[2]] <- rownames(x)[[i]]
    }
    }
  x <- mat.x
  }
x <- as.data.frame(x)
  
# rename column to be of the correct leaves (with population=TRUE values)
g <- which(tree$Get('Population', filterFun = isLeaf, traversal = 'post-order'))
pop.x <- x %>% 
  select(g)

pop.x <- round(pop.x)
peakedq0.pop.x <- pop.x


# generate values \hat{N} to plot
weights <- as.matrix(pop.ko.weights(tree))
peakedq0weights <- weights
peakedq0Ns <- as.data.frame(as.matrix(pop.x) %*% t(weights))
peakedq0Ns <- round(peakedq0Ns)
peakedq0Ns.restrict <- peakedq0Ns[peakedq0Ns$V1>=10 & peakedq0Ns$V1<=5000, , drop = FALSE] # restrict to the range from a=10 to b=750
p <- ggplot(peakedq0Ns.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for peaked distribution for q, uniform for p", subtitle = "Point mass on 0. One MM pass and set of weights, 5000 p,q estimates.")
p

#ggsave("altpeakedq1MMests.pdf")

Interestingly, we would have an essentially blank plot if we were to restrict to the range \((10,750)\) as much higher values of \(N\) are being estimated. In the summary for this case, we do expect to see larger estimates \(\hat{N}\) from \(B\), since the peaked distribution \(f(q)\) with mass on 0 suggests the population at \(B\) only represents a small fraction of the population at the parent of \(B\):

print("Uniform estimates from A,B:")
## [1] "Uniform estimates from A,B:"
summary(uniform.pop.x)
##        A                  B          
##  Min.   :    50.0   Min.   :    153  
##  1st Qu.:    66.0   1st Qu.:    391  
##  Median :    98.0   Median :    778  
##  Mean   :   393.3   Mean   :  11320  
##  3rd Qu.:   198.0   3rd Qu.:   2316  
##  Max.   :114781.0   Max.   :9319242
print("Peaked q estimates from A,B (mass on 0):")
## [1] "Peaked q estimates from A,B (mass on 0):"
summary(peakedq0.pop.x)
##        A                  B            
##  Min.   :    50.0   Min.   :3.561e+03  
##  1st Qu.:    67.0   1st Qu.:3.324e+04  
##  Median :   101.0   Median :8.565e+04  
##  Mean   :   489.1   Mean   :1.901e+06  
##  3rd Qu.:   212.0   3rd Qu.:2.643e+05  
##  Max.   :174975.0   Max.   :3.496e+09

This is indeed supported by the estimates of \(N\) from leaf \(B\). Although the overall weighted estimates are heavily effected by these larger values, the weighting of paths is not as intuitive as I would have hoped. To clarify, let’s consider the variance of estimates from \(A,B\). Weights are based on the sample precision matrix generated using the estimates \(\hat{N}\) from each leaf. Let \(q \sim Beta (150,1)\), the case with mass on 1, without loss of generality. Since \(p \sim Uniform(0,1)\), we expect \(Var(p) = 1/12 \approx 0.083\), while \(Var(q) = \frac{150}{(151)^2(152)} \approx 4.33*10^{-5}\). The effect of the longer path to leaf \(B\) does affect it’s overall weighting, placing more weight on the data from \(A\), but when comparing only estimates from \(B\), we should see smaller variation in the peaked case than the uniform case. I am not sure why this is occurring.

Considering we have strong data from leaf \(B\), I would have expected to see a relative increase in that weight from the uniform case, however this does not appear to be happening, and the coefficient of variation is actually higher for estimates from that leaf:

print("Uniform weighted estimates of N:")
## [1] "Uniform weighted estimates of N:"
summary(uniformNs)
##        V1           
##  Min.   :     88.0  
##  1st Qu.:    221.8  
##  Median :    416.0  
##  Mean   :   3975.7  
##  3rd Qu.:   1029.2  
##  Max.   :3055623.0
print("Weights in uniform case:")
## [1] "Weights in uniform case:"
uniformweights
##                         
## [1,] 0.6721284 0.3278716
print("Peaked q weighted estimates of N (mass on 0):")
## [1] "Peaked q weighted estimates of N (mass on 0):"
summary(peakedq0Ns)
##        V1           
##  Min.   :     1068  
##  1st Qu.:     9455  
##  Median :    24066  
##  Mean   :   527613  
##  3rd Qu.:    73903  
##  Max.   :969638574
print("Weights in peaked f(q) case (mass on 0):")
## [1] "Weights in peaked f(q) case (mass on 0):"
peakedq0weights
##                         
## [1,] 0.7226716 0.2773284
print("Ratio of coefficients of variation:")
## [1] "Ratio of coefficients of variation:"
cv.uniform <- sd(uniform.pop.x$B)/mean(uniform.pop.x$B)
cv.peaked <- sd(peakedq0.pop.x$B)/mean(peakedq0.pop.x$B)
cv.peaked/cv.uniform
## [1] 1.839072

We may get a better sense of the differences in these three plots if we do not cut the upper end of the range so low. We re-plot here with the upper end set to \(b=5000\), instead:

uniformNs.restrict <- uniformNs[uniformNs$V1>=10 & uniformNs$V1<=5000, , drop = FALSE]
p <- ggplot(uniformNs.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for uniform distribution for q, uniform for p", subtitle = "One MM pass and set of weights, 5000 p,q estimates")
peakedq1Ns.restrict <- peakedq1Ns[peakedq1Ns$V1>=10 & peakedq1Ns$V1<=5000, , drop = FALSE]
q <- ggplot(peakedq1Ns.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for peaked distribution for q, uniform for p", subtitle = "Mass on 1. One MM pass and set of weights, 5000 p,q estimates")
peakedq0Ns.restrict <- peakedq0Ns[peakedq0Ns$V1>=10 & peakedq0Ns$V1<=5000, , drop = FALSE]
r <- ggplot(peakedq0Ns.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for peaked distribution for q, uniform for p", subtitle = "Mass on 0. One MM pass and set of weights, 5000 p,q estimates")
p

q

r

Now, we add a “middle case” - that is, let \(q \sim Beta(15,15)\).

middleq <- treedata
middleq[3,3] <- 14
middleq[3,4] <- 28
middleq$Estimate <- middleq$Est
middleq$Population <- middleq$Pop
middleq
# create tree with peaked data
tree <- FromDataFrameNetwork(middleq)

# calculate estimates from A,B using 5000 simulated p,q values
weightedPopTree(tree,sample_length=5000,method='weightedEstimate')
## [1] "using weighted mean estimates"
tree$Do(function(node) {
  if(node$isRoot){
    node$uncertainty <-pop.conf.int(tree)
    }
  else{
    node$uncertainty <- prettyConfidenceIntervals(node$overdose_samples)
    }}
  )

# create dataframe of N estimates from each leaf
x <- tree$Get('overdose_samples', filterFun = isLeaf, traversal = 'post-order')

# Ignore NAs and then also choose leaves with population=TRUE values only (these are the data)
if(is.list(x)){    
  mat.x <- NULL
  for (i in 1:length(x)) {
    if (length(x[[i]] > 0)) {
      mat.x <- cbind(mat.x, x[[i]])
      colnames(mat.x)[dim(mat.x)[2]] <- rownames(x)[[i]]
    }
    }
  x <- mat.x
  }
x <- as.data.frame(x)
  
# rename column to be of the correct leaves (with population=TRUE values)
g <- which(tree$Get('Population', filterFun = isLeaf, traversal = 'post-order'))
pop.x <- x %>% 
  select(g)

pop.x <- round(pop.x)
middleq.pop.x <- pop.x


# generate values \hat{N} to plot
weights <- as.matrix(pop.ko.weights(tree))
middleqweights <- weights
middleqNs <- as.data.frame(as.matrix(pop.x) %*% t(weights))
middleqNs <- round(middleqNs)
middleqNs.restrict <- middleqNs[middleqNs$V1>=10 & middleqNs$V1<=750, , drop = FALSE] # restrict to the range from a=10 to b=750
p <- ggplot(middleqNs.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for central distribution for q, uniform for p", subtitle = "Beta(15,15). One MM pass and set of weights, 5000 p,q estimates.")
middleqNs.restrict <- middleqNs[middleqNs$V1>=10 & middleqNs$V1<=5000, , drop = FALSE] # restrict to the range from a=10 to b=750
q <- ggplot(middleqNs.restrict, aes(x = V1)) + geom_density(alpha=.2, fill="#FF6666") + 
   ggtitle("MM estimates for central distribution for q, uniform for p", subtitle = "Beta(15,15). One MM pass and set of weights, 5000 p,q estimates.")
p

q

The mode here looks like it sits just below 250, while using the distribution \(f(N|A,B)\), the mode occurs closer to 350. The MM is also much more uncertain about the prospect of \(N\) having a higher vale.

Now we take a look at the summary of \(\hat{N}\) estimates from \(B\) in the middle case, and the coefficient of variation:

print("Uniform estimates from A,B:")
## [1] "Uniform estimates from A,B:"
summary(uniform.pop.x)
##        A                  B          
##  Min.   :    50.0   Min.   :    153  
##  1st Qu.:    66.0   1st Qu.:    391  
##  Median :    98.0   Median :    778  
##  Mean   :   393.3   Mean   :  11320  
##  3rd Qu.:   198.0   3rd Qu.:   2316  
##  Max.   :114781.0   Max.   :9319242
print("Middle case q estimates from A,B:")
## [1] "Middle case q estimates from A,B:"
summary(middleq.pop.x)
##        A                  B         
##  Min.   :    50.0   Min.   :   210  
##  1st Qu.:    67.0   1st Qu.:   402  
##  Median :   101.0   Median :   613  
##  Mean   :   409.4   Mean   :  2325  
##  3rd Qu.:   202.0   3rd Qu.:  1231  
##  Max.   :235380.0   Max.   :625175
print("Uniform weighted estimates of N:")
## [1] "Uniform weighted estimates of N:"
summary(uniformNs)
##        V1           
##  Min.   :     88.0  
##  1st Qu.:    221.8  
##  Median :    416.0  
##  Mean   :   3975.7  
##  3rd Qu.:   1029.2  
##  Max.   :3055623.0
print("Weights in uniform case:")
## [1] "Weights in uniform case:"
uniformweights
##                         
## [1,] 0.6721284 0.3278716
print("Middle case q weighted estimates of N:")
## [1] "Middle case q weighted estimates of N:"
summary(middleqNs)
##        V1          
##  Min.   :   131.0  
##  1st Qu.:   270.0  
##  Median :   411.0  
##  Mean   :  1346.8  
##  3rd Qu.:   803.5  
##  Max.   :306021.0
print("Weights in middle case:")
## [1] "Weights in middle case:"
middleqweights
##                         
## [1,] 0.5105496 0.4894504
print("Ratio of coefficients of variation:")
## [1] "Ratio of coefficients of variation:"
cv.uniform <- sd(uniform.pop.x$B)/mean(uniform.pop.x$B)
cv.middle <- sd(middleq.pop.x$B)/mean(middleq.pop.x$B)
cv.middle/cv.uniform
## [1] 0.3958381

We see a corresponding increase in the weight on path \(B\), a sign that the method is recognizing there is more information contained in this branch than in the uniform case. However, this was not echoed in the peaked case with mass on 0, so we should not be so quick to celebrate this.

cv.peakedq1 <- sd(peakedq1.pop.x$B)/mean(peakedq1.pop.x$B)
cv.peakedq1/cv.uniform
## [1] 1.327889

In fact, the COV is higher in all cases that are not uniform. Why is this so?