##' @title Dirichlet distribution: random generator
##' @param n Number of draws
##' @param alpha Dirichlet parameter: a vector of positive number
##' @return A matrix with \code{n} rows and \code{length(alpha)} columns
##' @export
rdirichlet <- function (n=1, alpha)
{
if(is.vector(alpha))
l <- length(alpha)
else
l <- nrow(alpha)
x <- matrix(rgamma(l * n, alpha), ncol = l, byrow = TRUE)
sm <- apply(x, 1, sum)
return(x/ sm )
}
##' @rdname dpairbeta
##' @param n The number of points on the simplex to be generated.
##' @param dimData the dimension of the sample space, which is \eqn{1 + } the dimension of the simplex.
##' @keywords datagen distribution
##' @export
rpairbeta <-
function(n=1, dimData=3, par= c(1, rep(1,3)) )
{
p <- dimData
if((choose(p,2)+1) != length(par) )
{ stop("wrong parameter length")}
alpha <- par[1]
beta <- par[-1]
res=matrix(0,ncol=p,nrow=n) #to store the result
for(ll in 1:n)
{
##choose uniformly a pair i,j:
i_0=floor(runif(1,1,(p+1)))
jStar=floor(runif(1,1,p))
if(jStar<i_0){j_0=jStar}else{j_0=jStar+1}
i=min(i_0,j_0);j=max(i_0,j_0)
place=(i-1)*(p-(i)/2)+j-i ##place=\sum_{k=0}^{i_1}(p-k) + (j-i)
##simulate r=wi+wj and theta=wi/r
r_ij=rbeta(1,2*alpha+1,(p-2)*alpha) #2*alpha+1
theta=rbeta(1,beta[place],beta[place])
wi=theta*r_ij
wj=(1-theta)*r_ij
#simulate the remaining components (p-2), uniformly on the simplex : "sum(w_i)=1-r"
remains=(1-r_ij)*rdirichlet(1,rep(1,(p-2)))
#put together r,theta and the remainings
W0=append(as.vector(remains),wi,after=(i-1))
W=append(W0,wj,after=(j-1))
#keep the result...
res[ll,]= W
}
return(res)
}
## @importFrom MCMCpack rdirichlet
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.