#### Function for simulating data
generate.data = function(alphax, betax, alphay, betay, eta, N) {
  z = rnorm(N,0,1)
  x = y = numeric(N)

  ## px is an N x length(alphax) matrix.
  ## Each row has the TRUE cummulative probabilities for each subject.
  px = (1 + exp(- outer(alphax, betax*z, "+"))) ^ (-1)
  aa = runif(N)
  for(i in 1:N)
    x[i] = sum(aa[i] > px[,i])
  x = as.numeric(as.factor(x))
  ## x = x+1 may have category gaps if there are small probability categories.

  py = (1 + exp(- outer(alphay, betay*z+eta[x], "+"))) ^ (-1)
  aa = runif(N)
  for(i in 1:N)
    y[i] = sum(aa[i] > py[,i])
  y = as.numeric(as.factor(y))
  ## y = y+1 may have category gaps if there are small probability categories.

  return(list(x=x, y=y, z=z))
}
  
N = 500
Nemp = 1000
NREPL = 1000
#
alphay = c(-1, 0, 1)
betay = -.5
alphax = c(-1, 0, 1, 2)
betax = 1
#
eta = rep(0,5)

data = generate.data(alphax, betax, alphay, betay, eta, N)
#COBOT.scores(data)