#############################################################
  # data BiRepeated was simulated using the following codes
  #############################################################
  set.seed(10)
  n <- 200
  x <- seq(0.1, 0.6, 0.1)
  K <- length(x)
  sum1 <- exp(2.5)
  sum2 <- exp(3)
  
  library(mvtnorm)
  set.seed(3)
  uu <- matrix(0,n,2)
  part1 = NULL
  part2 = NULL
  
  y1 <- rep(0, n*K)
  y2 <- rep(0, n*K)
  
  for(i in 1:n){
    S <- matrix(c(0.2,0,0,0.2),2,2)
    u <- rmvnorm(1, c(0,0), S)
    uu[i,] <-u
    for(j in 1:K){ 
      # the regression coeff can not be ~0, o.w, the prediction is bad. 
      lp1 <- u[1]-1 + (u[2]+1)*x[j] 
      a11 <- sum1*exp(lp1)/(1+exp(lp1))
      a12 <- sum1/(1+exp(lp1))
      y1[(i-1)*K+j] <- rbeta(1,a11, a12)
      
      lp2 <- u[1]-2 + (u[2]+2)*x[j]
      a21 <- sum2*exp(lp2)/(1+exp(lp2))
      a22 <- sum2/(1+exp(lp2))
      y2[(i-1)*K+j] <- rbeta(1,a21, a22)     
    }
  }
  
  y1 <- round(y1,3)
  y2 <- round(y2,3)
  id <- rep(1:n, each=K)
  x <- rep(x, n)
  BiRepeated <- data.frame(id, y1,y2,x)