n <- 2500

x.use <- rnorm(n)
w.use <- runif(n,-1,1)
marginals.use <- c("ZINB", "ZIGA")

# simulate data
y.use <- scdeco.sim.cop(marginals=marginals.use, x=x.use,
                    eta1.true=c(-2, 0.8), eta2.true=c(-2, 0.8),
                    beta1.true=c(1, 0.5), beta2.true=c(1, 1),
                    alpha1.true=7, alpha2.true=3,
                    tau.true=c(0, 1), w=w.use)
                    
                    
######## plotting the data                     

library(gridExtra)
library(ggplot2)
grid.arrange(p5, p6, ncol = 3)


group<-as.numeric(cut(x.use, breaks=quantile(x.use, prob=seq(0, 1, by=0.16))))

dat<-data.frame(y.use, x.use, group)
colnames(dat)<-c("Gene A", "Gene B", "X", "group")


p1<-ggplot(data = subset(dat, group==1),  
  aes(x =`Gene A`, y = `Gene B`)) +  
  geom_point()
p1

p2<-ggplot(data = subset(dat, group==2),  
  aes(x =`Gene A`, y = `Gene B`)) +  
  geom_point()

p3<-ggplot(data = subset(dat, group==3),  
  aes(x =`Gene A`, y = `Gene B`)) +  
  geom_point()
p4<-ggplot(data = subset(dat, group==4),  
  aes(x =`Gene A`, y = `Gene B`)) +  
  geom_point()
p5<-ggplot(data = subset(dat, group==5),  
  aes(x =`Gene A`, y = `Gene B`)) +  
  geom_point()
p6<-ggplot(data = subset(dat, group==6),  
  aes(x =`Gene A`, y = `Gene B`)) +  
  geom_point()

  
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3)

####### Estimation 

mcmc.out <- scdeco.cop(y=y.use, x=x.use, marginals=marginals.use, w=w.use,
                     n.mcmc=5000, burn=1000, thin=10)
                     
                     
                     lowerupper <- t(apply(mcmc.out, 2, quantile, c(0.025, 0.5, 0.975)))
estmat <- cbind(lowerupper[,1],
                c(c(-2, 0.8), c(-2, 0.8), c(1, 0.5), c(1, 1), 7, 3, c(-0.2, .3)),
                lowerupper[,c(2,3)])
colnames(estmat) <- c("lower", "trueval", "estval", "upper")
estmat