Skip to main content
edited body
Source Link
Zen
  • 25.1k
  • 4
  • 87
  • 126

Essentially the same idea of Bill Huberwhuber, nowthis time using Dirichlets: the lines of the matrix are independent Dirichlets, with each Dirichlet putting more mass on the corresponding diagonal element.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if x does not satisfy the constraints. If the rejection rate is too high, increase the value of A0. Here is a sample x.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490

P.S. VectorizePlease, vectorize that pesky for.

Essentially the same idea of Bill Huber, now using Dirichlets.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if x does not satisfy the constraints. If the rejection rate is too high, increase the value of A0. Here is a sample x.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490

P.S. Vectorize that pesky for.

Essentially the same idea of whuber, this time using Dirichlets: the lines of the matrix are independent Dirichlets, with each Dirichlet putting more mass on the corresponding diagonal element.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- A0
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if x does not satisfy the constraints. If the rejection rate is too high, increase the value of A0. Here is a sample x.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490

P.S. Please, vectorize that pesky for.

added 42 characters in body
Source Link
Zen
  • 25.1k
  • 4
  • 87
  • 126

Essentially the same idea of Bill Huber, now using Dirichlets.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if xx does not satisfy the constraints. If the rejection rate is too high, increase the value of A0A0. Here is a sample xx.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490

P.S. Vectorize that pesky for.

Essentially the same idea of Bill Huber, now using Dirichlets.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if x does not satisfy the constraints. If the rejection rate is too high, increase the value of A0. Here is a sample x.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490

Essentially the same idea of Bill Huber, now using Dirichlets.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if x does not satisfy the constraints. If the rejection rate is too high, increase the value of A0. Here is a sample x.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490

P.S. Vectorize that pesky for.

Source Link
Zen
  • 25.1k
  • 4
  • 87
  • 126

Essentially the same idea of Bill Huber, now using Dirichlets.

rdirichlet <- function(a) {
    x <- sapply(a, function(a) rgamma(1, a, 1))
    return(x / sum(x))
}

n <- 5
A0 <- 10
a <- matrix(data = 1, nrow = n, ncol = n)
diag(a) <- rep(A0, n)
x <- matrix(nrow = n, ncol = n)
for (i in 1:n) x[i,] <- rdirichlet(a[i,])

Reject if x does not satisfy the constraints. If the rejection rate is too high, increase the value of A0. Here is a sample x.

> x
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,] 0.484437196 0.03511667 0.04919437 0.12718905 0.304062711
[2,] 0.008626386 0.76520244 0.03231747 0.02454215 0.169311552
[3,] 0.176631389 0.01971251 0.55780424 0.07952712 0.166324747
[4,] 0.003109732 0.12056624 0.09335086 0.77330892 0.009664246
[5,] 0.037015097 0.02485376 0.04536731 0.05083834 0.841925490