Where does ward.d2 and cubic clustering start misclassifying for small sample sizes?

Question

Given
Two univariate Gaussian distributions with unit variance and means and separated by distance 'dx', from which samples of size N are drawn from each distribution.

Find
the minimum value of dx such that on ~300 repeated tests using Ward's clustering with Euclidean distance, and the cubic cluster criterion there is a greater than 99% chance of not having a false classification as a function of "N". Assume that one of the clusters has its center at 0, and therefore the other is at 'dx', that the maximum value of dx is less than 10, and the range of N is between 5 and 35.

I can simulate this, and get an approximation, but I am wondering if there is an analytic form derived from first principles.

Answer 1

So here is how I initially made the data:

## libraries and purpose
library(NbClust)
library(pracma)

## main
set.seed(1)

#how many trials per sample-size
N_tests <- 3000

#list of sample sizes
N_samp <- c(seq(from=5, to=25, by=1),
            seq(from=30, to=100, by=5))

#range of separation (dx)
separation <- seq(from=2.5,to=6.5,length.out=100)

#initialization
err_sum <- matrix(0,nrow = length(separation),ncol=length(N_samp))

#main loop
for(k in 1:length(N_samp)){

     for(i in 1:length(separation)){

          target <- c(rep.int(1,N_samp[k]),
                      rep.int(2,N_samp[k]))

          err <- numeric(length(target)/2)

          for (j in 1:N_tests){

               #make data
               x<-c(rnorm(N_samp[k],mean = 0,sd=1),
                    rnorm(N_samp[k],mean=separation[i],sd=1) )

               #cluster it
               res<-NbClust(x,
                            distance = "euclidean",
                            min.nc=1, max.nc=4,
                            method = "ward.D2",
                            index = "ccc")

               #count misclassification
               idx <- which(res$Best.partition - target != 0, arr.ind = T)

               if (length(idx)==0){
                    err[j] <- 0
               } else {
                    err[j] <- length(idx)/length(target)
               }
          }

          err_sum[i,k] <- mean(err)
     }

}

data1 <- data.frame(separation,err_sum)

for (i in 1:(length(N_samp))){
     names(data1)[1+i] <- paste("n_", as.character(N_samp[i]),sep="")

     zidx <- which(data1[,i+1]>= 1/N_samp[i])

     zval <- (0*data1[,i+1])
     zval[zidx] <- 1

     data1 <- cbind(data1,zval)
     names(data1)[ncol(data1)] <- paste("z_", as.character(N_samp[i]),sep="")

}

for (i in 1:(length(N_samp))){

     dzval <- c(0, diff(data1[,37+i]))

     data1 <- cbind(data1,dzval)
     names(data1)[ncol(data1)] <- paste("dz_", as.character(N_samp[i]),sep="")

}

write.csv(x = data1,file = "myout.csv")

This gave me a csv file, which I imported into JMP, stacked the DZ columns, looked for non-zero values, and then averaged them to get a table of separation that has 50/50 likelihood of not containing a mis-classification, and then fit it in the bivariate platform as shown here:

I think the log-log domain, third order polynomial fit is the best one. I don't believe it is much more than a transformed version of a taylor series approximation to the truth, but it is what I have.

$ dx \left( N,P \right) |_{P = 50 \%} = e^{1.3989 - 0.1221 \cdot ln \left(N\right) + 0.07444 \cdot ln \left(N\right)^{2} - 0.007329 \cdot ln \left(N\right)^{3}}$

And yes, I used the level of local oscillations in the spline to determine a decent fit parameter. Someday I will apply a Fourier method to it to find the best parameter, or I will find someone else who has done it, but until then I will use it as an "eyeball norm" to augment fit parameters.

The sampling delta for x was about 0.04 units. I would expect error in the mean to be on that scale. I would also expect the trendline error to be less than that.