## R CODE FOR SAMPLING 3-dimensional Beta(nu1,nu2) 
## with given pairwise correlation coefficients
## and Beta(nu1,nu2) marginals
## 
## Written by Vanja Dukic and Nevena Maric
## Please cite: 
## Dukic, V. and Maric, N., "Minimum correlation in construction of multivariate distributions”, Physical Review E, 87(3):032114
## 
## NOV 20 2012
##
## README.txt
####################################################################

# EXAMPLE OF USE:
source('cphi.r')
source('beta3D.r')

# specify the desired sample size, dimension, and parameters
N=10000; d = 3; nu1 = 4; nu2 = 7; 

# specify the desired pairwise correlation parameters
r12=0.4
r13=0.3
r23=0.2

# sampling
y=beta3D(nsamples=N,nu1=nu1,nu2=nu2,r12=r12, r23=r23, r13=r13,d=3)

# check the sample pairwise correlation coefficients 
cor(y[,1],y[,2])
cor(y[,1],y[,3])
cor(y[,2],y[,3])


## VISUALIZATION 
install.packages('scatterplot3d')
require(scatterplot3d)
library(scatterplot3d)

par(mfrow=c(1,1),cex.axis=2,cex=1,cex.lab=2)
scatterplot3d(y[,1],y[,2],y[,3], angle=60,color='purple',pch=16,  highlight.3d=TRUE, #main="3D Beta, with correlation 0.25",
xlim = c(0,1), ylim = c(0,1), zlim = c(0,1), xlab = 'First component', ylab = '', zlab = 'Third component')

x11()
par(mfrow=c(1,4))
hist(rbeta(10000,nu1,nu2))
hist(y[,1])
hist(y[,2])
hist(y[,3])