sum to zero transform

[bob-carpenter]

@spinkney proposed a new sum-to-zero transform:

stan-dev/stanc3#1111 (comment)

Given time, we should evaluate this as well as the one we currently use (which sets the last value to the negation of the sum of the first values).

[spinkney]

There's another one

d <- 4

B <- matrix(NA, nrow = d, ncol = d)
B[d, ] <- rep((1/sqrt(d)), d)

for (j in 1:(d-1)){
  for (i in 1:d){
    if (i < j)
      B[j, i] <- 0
    if (i == j)
      B[j, i] <- ((d - j)/(d - j + 1))^(1/2)
    if (i > j)
      B[j, i] <- - 1/(((d-j+1)*(d-j))^(1/2))
  }
}

# this is just taking 1000 rows of size d and showing
# that average deviation from 0 is really small
mean(cbind(matrix(runif((d-1) * 1000), 1000, d - 1), 0) %*% B)
-6.526446e-16

In Stan you'd construct the B matrix and then a d - 1 vector. In tp, append 0 to the vector and then multiply the vector by the matrix B to get the sum-to-zero constraint.

Cite: user277126 (https://stats.stackexchange.com/users/277126/user277126), Sampling Normal variables with linear constraints and given variances - Fraser (1951), URL (version: 2021-12-28): https://stats.stackexchange.com/q/558510

Which in turn comes from
Fraser, D. (1951). Normal Samples With Linear Constraints and Given Variances. Canadian Journal of Mathematics, 3, 363-366. doi:10.4153/CJM-1951-041-9

[mjhajharia]

thanks @spinkney I'll check this out today!!

[spinkney]

I think "best" for a random sum-to-zero uniform vector is using this program below which I ported from Matlab at https://www.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum. The method of generating a d - 1 vector and then subtracting the sum does not preserve uniformity.

In the comments of that Matlab code there is mention of another program when you don't have symmetrical bounds (which I don't think we care about). You can read about that at https://www.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum?tab=discussions#discussions_1687518.

randfixedsum <- function(n, a, b) {
 # m <- samples
  m <- 1  # number of samples, just set to 1
  s <- 0  # constraint on sum
  
  if ( n < 1 )
   return('n must be a whole number and m a non-negative integer.')
  
  if (a >= b )
   return('a must be < b')

  # Rescale to a unit cube: 0 <= x(i) <= 1
  s <- (s - n * a) / (b - a)
  
  # Construct the transition probability table, t.
  # t[i,j] will be utilized only in the region where j <= i + 1.
  k <- max(min(floor(s), n - 1), 0)  # Must have 0 <= k <= n-1
  s <- max(min(s, k + 1), k) # Must have k <= s <= k+1
  s1 <- s - k:(k-n + 1) # s1 & s2 will never be negative
  s2 = (k + n):(k + 1) - s
  w = matrix(0, n, n + 1)
  realmax <- .Machine$double.xmax
  w[1, 2] <- realmax
 
  t <- matrix(0, n - 1, n) 
  
  tiny <- .Machine$double.xmin
  for (i in 2:n) {
    tmp1 <- w[i - 1, 2:(i+1)] * s1[1:i] / i
    tmp2 <- w[i - 1, 1:i] * s2[(n-i+1):n] / i
    w[i, 2:(i+1)] <- tmp1 + tmp2
    tmp3 <- w[i, 2:(i+1)] + tiny # In case tmp1 & tmp2 are both 0,
    tmp4 <- (s2[(n - i + 1):n] > s1[1:i]) # then t is 0 on left & 1 on right
    t[i-1, 1:i] <- (tmp2 / tmp3) * tmp4 + (1 - tmp1 / tmp3 ) * !tmp4
  }
  
  # Derive the polytope volume v from the appropriate
  # element in the bottom row of w.
  v <- n^(3/2) * (w[n, k+2] / realmax) * (b - a)^(n - 1)
  # Now compute the matrix x.
  x <- matrix(0, n, m)
  if(m == 0)
    return(x) # If m is zero, quit with x = []
  
  rt = matrix(runif((n - 1) * m), n - 1, m) # For random selection of simplex type
  rs = matrix(runif((n - 1) * m), n - 1, m) # For random location within a simplex
  s = matrix(rep(s, m), 1, m)
  j = matrix(rep(k + 1), 1, m) # For indexing in the t table
  sm = matrix(0, 1, m)
  pr = matrix(1, 1, m) # Start with sum zero & product 1
  
  for (i in (n - 1):1) {  # Work backwards in the t table
    e = rt[n-i,] <= t[i,j] # Use rt to choose a transition
    sx = rs[n-i,] ^ (1/i) # Use rs to compute next simplex coord.
    sm = sm + (1-sx) *pr * (s/(i+1)) # Update sum
    pr = sx *pr # Update product
    x[n-i,] = sm + pr*e #Calculate x using simplex coords.
    s = s - e
    j = j - e #Transition adjustment
  }
  
  x[n,] = sm + pr * s # Compute the last x
  # Randomly permute the order in the columns of x and rescale.
  rp = matrix(runif(n * m), n, m) # Use rp to carry out a matrix 'randperm'
  p = sort(rp, index.return = T)$ix # The values placed in ig are ignored
  mat <- matrix(rep(seq(from = 0, to = n * (m - 1), by = n), n), n, m, byrow = T)
  x = (b-a) * x[p + mat] + a # % Permute & rescale x
  return( x )
}

sample_mat <- matrix(0, 10, 100)
for (i in 1:100) {
sample_mat[, i] <- randfixedsum(10, -5, 5)
}

hist(sample_mat)
colSums(sample_mat)
colSums(sample_mat)
  [1] -4.440892e-16  6.217249e-15 -1.065814e-14  4.440892e-15  4.440892e-15 -8.881784e-16  0.000000e+00
  [8] -4.440892e-15 -3.108624e-15  3.552714e-15 -4.440892e-16 -4.440892e-16  0.000000e+00  1.776357e-15
 [15]  1.332268e-15 -2.664535e-15  2.220446e-15  0.000000e+00  8.437695e-15 -4.884981e-15 -8.881784e-16
 [22]  0.000000e+00  8.881784e-16 -3.552714e-15 -4.440892e-15 -4.440892e-16  2.220446e-15  8.881784e-16
 [29]  0.000000e+00  3.108624e-15  7.105427e-15  9.769963e-15  0.000000e+00  1.776357e-15 -4.440892e-16
 [36] -4.884981e-15 -3.996803e-15 -1.332268e-15  2.664535e-15 -8.881784e-16 -1.776357e-15  3.552714e-15
 [43] -2.664535e-15 -3.108624e-15  0.000000e+00  8.881784e-16 -3.552714e-15  3.108624e-15 -4.884981e-15
 [50]  5.329071e-15  1.776357e-15 -5.773160e-15  1.776357e-15  2.664535e-15  2.664535e-15 -7.993606e-15
 [57]  8.881784e-16 -7.993606e-15 -1.776357e-15 -6.217249e-15  0.000000e+00  0.000000e+00 -8.881784e-16
 [64] -1.776357e-15 -4.440892e-15  1.776357e-15 -3.996803e-15  1.776357e-15  2.220446e-15  0.000000e+00
 [71] -2.664535e-15  1.776357e-15 -4.884981e-15 -1.776357e-15  2.220446e-15  1.776357e-15  3.552714e-15
 [78]  1.776357e-15 -2.664535e-15 -8.881784e-16 -8.881784e-16 -3.552714e-15 -1.332268e-15  6.661338e-15
 [85] -4.440892e-16  0.000000e+00 -8.881784e-16 -3.108624e-15  7.105427e-15  3.996803e-15  3.108624e-15
 [92]  2.664535e-15 -3.108624e-15  8.881784e-16  0.000000e+00 -1.332268e-15 -2.664535e-15 -8.881784e-16
 [99]  4.440892e-15 -8.881784e-16
> 

Otherwise, I can generate a sum-to-zero vector with a standard normal distribution via

d <- 5
v <- matrix(-1/(d - 1), d, d)
diag(v) <- 1
eigen(v)

x <- mvrnorm(1, mu=rep(0,d), Sigma=v)
sum(x)

But I think we want the uniform method.