non-conformable arguments

Hi, I'm having problems with this error: "non-conformable arguments". I know it means that I am multiplying matrices with different dimensions, but I'm not sure how I can stop it from happening given the code works for a few iterations and then fails. Any help much appreciated! Thanks!

library("MASS")     # mvrnorm() for multivariate normal
set.seed(15)        # random seed set for reproducibility
n <- 60              # sampling 60 each from 4 separate distributions
m1 <- c(1.5, 1.5)   # upper right x, y means
S1 <- matrix(c(0.3, 0.05, 0.05, 0.3), ncol = 2) # variance of c1
# sampling n from each cluster as per its mean mu and variance Sigma
clus1 <- mvrnorm(n = n, mu = m1, Sigma = S1)
m2 <- c(1.5, -1.5)  # lower right
S2 <- matrix(c(0.5, -0.08, -0.08, 0.2), ncol = 2)
clus2 <- mvrnorm(n = n, mu = m2, Sigma = S2)
m3 <- c(-1.5, 1.5)  # upper left
S3 <- matrix(c(0.1, 0.03, 0.03, 0.1), ncol = 2)
clus3 <- mvrnorm(n = n, mu = m3, Sigma = S3)
m4 <- c(-1.5, -1.5) # lower left
S4 <- matrix(c(0.8, 0.50, 0.50, 0.8), ncol = 2)
clus4 <- mvrnorm(n = n, mu = m4, Sigma = S4)
datc <- rbind(clus1, clus2, clus3, clus4) # 240 observations altogether
# run the CRP Gibbs function in Appendix B.
alpha <- 0.01
mu0 <- matrix(rep(0, 2), ncol = 2, byrow = TRUE)
sigma0 <- diag(2) * 3^2
sigma_y <- diag(2) * 1
c_init <- rep(1, nrow(datc))
results <- crp_gibbs(data = datc, alpha = alpha, mu0 = mu0,
                     sigma0 = sigma0, sigma_y = sigma_y, c_init = rep(1, nrow(datc)))

crp_gibbs <- function(data,alpha=0.01,mu0,sigma0,sigma_y,c_init,maxIters=1000)
{
  # data: an N by D matrix of data points
  # A small alpha encourages fewer clusters
  #    alpha <- 0.01
  # sigma_y: measurement error of data y, assumed known, same for all clusters.
  #    sigma_y <- diag(data_dim) * 1
  # mu0, sigma0: prior mean and variance around unknown mean mu
  #    mu0 <- matrix(rep(0, data_dim), ncol = data_dim, byrow = TRUE)
  #    sigma0 <- diag(data_dim) * 3^2
  # c_init: initial assignments of data points to clusters
  #
  require(mvtnorm)
  # dimension of the data points
  data_dim <- ncol(data)
  # number of data points
  N <- nrow(data)
  #####
  # Priors
  #####
  # prior precision on the unknown mean mu.
  tau0 <- solve(sigma0) # prior precision on $mu$, inverse of prior covariance
  # cluster-specific precision, assumed known, all clusters are assumed to
  # share identical measurement error of y ~ N(mu, sigma_y).
  tau_y <- solve(sigma_y)
  # initialize the CRP Gibbs sampler
  z <- c_init # initial cluster membership assignments
  n_k <- as.vector(table(z)) # initial data counts at each cluster, Eq (4).
  Nclust <- length(n_k) # initial number of clusters
  ##

  # Chinese Restaurant Process (CRP) Gibbs sampler begins
  ##
  res <- matrix(NA, nrow = N, ncol = maxIters) # cluster membership storage
  pb <- txtProgressBar(min = 0, max = maxIters, style = 3)
  for(iter in 1:maxIters) { # maxIters also prevents endless loops
    for(n in 1:N) { # one data point (customer) at a time
      # when nth customer enters the Chinese restaurant, we need to first
      # un-assign his/her initial cluster membership, then use Eq (12)
      # [already occupied table] and (13) [new table] to calculate the
      # updated probability of table assignment.
      c_i <- z[n] # what is the nth persons table assignment?
      n_k[c_i] <- n_k[c_i] - 1 # remove the nth person from table
      # if the table becomes empty when the nth person is removed,
      # then that table/cluster is removed.
      if( n_k[c_i] == 0 )
      {
        n_k[c_i] <- n_k[Nclust]  # last cluster to replace this empty cluster
        loc_z <- ( z == Nclust ) # who are in the last cluster?
        z[loc_z] <- c_i          # move them up to fill just emptied cluster
        n_k <- n_k[ -Nclust ]    # take out the last cluster, now empty
        Nclust <- Nclust - 1     # decrease total number of clusters by 1
      }
      z[n] <- -1 # ensures z[n] doesnt get counted as a cluster #####
      ##
      # Now we are ready to update table assignment by Eqs (12) and (13).
      ##
      # log probabilities for the clusters, add previously unoccupied table
      logp <- rep( NA, Nclust + 1 )
      # loop over already occupied tables 1:J and calculate pr as per Eq (13).
      for( c_i in 1:Nclust ) {
        tau_p <- tau0 + n_k[c_i] * tau_y # cluster precision as per Eq (4)
        sig_p <- solve(tau_p) # cluster variance, inverse of tau_c
        # find all of the points in this cluster
        loc_z <- which(z == c_i)
        # sum all the points in this cluster
        if(length(loc_z) > 1) {
          sum_data <- colSums(data[z == c_i, ]) }
        else {
          sum_data <- data[z == c_i, ]
        }
        #
        # We need to use the predictive distribution of each already
        # occupied table to predict the next customer sitting there.
        #
        # Use Eq (4) to find the conditional posterior distribution for
        # the cluster means, (y * n_k * tau_j + mu0 * s0) /tau_p,and
        # then use the predictive distribution of y_j in Eq (12) to
        # predict new data value c_i from c-i.
        #
        mean_p <- sig_p %*% (tau_y %*% sum_data + tau0 %*% t(mu0))
        logp[c_i] <- log(n_k[c_i]) +
          dmvnorm(data[n,], mean = mean_p, sigma = sig_p + sigma_y, log = TRUE) #Eq (12) }
      #
      # We are done looping over already occupied tables. Next, we use
      # Eq (13) to calcualte the log probability of a previously
      # unoccupied, 'new' table. Essentially,it is the prior predicitive
      # distribution of the DP.
      logp[ Nclust+1 ] <- log(alpha) +
        dmvnorm(data[n,], mean = mu0, sigma = sigma0 + sigma_y, log = TRUE) #Eq (13)
      # transform unnormalized log probabilities into probabilities
      max_logp <- max(logp)
      logp <- logp - max_logp
      loc_probs <- exp(logp)
      loc_probs <- loc_probs / sum(loc_probs)
      
      
      if(is.na(loc_probs)) {
        print(c_i)
        print(n)
        print(iter)
        print(loc_probs)
      }
      
      if(is.na(loc_probs)){
        loc_probs[is.na(loc_probs)] <-0.01
      }
      # draw a sample of which cluster this new customer should belong to
      newz <- sample(1:(Nclust+1), 1, replace = TRUE, prob = loc_probs)
      # spawn a new cluster if necessary
      if(newz == Nclust + 1) {
        n_k <- c(n_k, 0)
        Nclust <- Nclust + 1
      }
      z[n] <- newz
      n_k[newz] <- n_k[newz] + 1 # update the cluster n_k
    }
    setTxtProgressBar(pb, iter) # update text progress bar after each iter
    res[, iter] <- z     # cluster membership of N observations
  }
  close(pb)            # close text progress bar
  invisible(res)       # return results, N by maxIters matrix
}

}

Hi, when I run this I get an

Error in crp_gibbs(data = datc, alpha = alpha, mu0 = mu0, sigma0 = sigma0, : could not find function "crp_gibbs"

so I'm guessing that crp_gibbs is a custom function (I see a reference to 'Appendix B' so I'm guessing it's not part of a package). Since this seems to be the function that's causing your error (I don't get an error from any but the last line), are you able to share that function?

Also, is it your function or someone else's? If it's someone else's (like the author of a paper) you might be better asking them directly.

Hi Gabriel,

The function is included in the code above. I run the function first and then the code at the top.

Thanks

Seems to be from CRPClustering:

@Dr_Loz I am very very sorry - I didn't scroll down and because the call to your crp_gibbs is before the definition it threw an error. I see it now.

@williaml I looked at the source for CRPClustering::crp_gibbs and it's not the same; this looks like a function with the same name maybe implementing the same algorithm. With that said, the CRPClustering package does give a sample dataset which looks very similar to the example data here but drawn from 3 distributions. That also gives an error if I try to run it through this thread's crp_gibbs function.

I don't have an answer at this point but I can observe that the error is coming from tau_y %*% sum_data and when I ran the function on that other dataset there was a requires numeric/complex matrix/vector arguments error rather than the conformable arguments error with datc.

tau_y is a fixed 2x2 (in this case) matrix so that's not going to be the problem. sum_data might be running into a condition where it's logical:

solve(diag(2)) %*% TRUE
# Error in solve(diag(2)) %*% TRUE : non-conformable arguments

(the same happens for FALSE). Perhaps there's a place where a 0 inadvertently gets turned into a FALSE or 1 into TRUE.

NULL gives that other error:

Error in solve(diag(2)) %*% NULL
# Error in solve(diag(2)) %*% NULL : 
#   requires numeric/complex matrix/vector arguments

I hope this gives you something to go on - let me know if that helps at all and I can have another look through it, but can't make any promises.

1 Like