Generación de matrices aleatorias con restricciones de suma y maximalidad

Question

Me gustaría generar una matriz cuadrada aleatoria tal que las filas estén normalizadas a uno y los elementos diagonales sean el máximo de su columna. ¿Existe una forma eficiente de muestrear estas matrices de manera uniforme?

$2 \times 2$ Las matrices son fáciles de crear. La primera columna se genera mediante el muestreo de dos uniformes independientes de $[0,1]$ y moviendo el máximo a la primera entrada. La segunda columna es entonces el complemento de la primera. Mis intentos de seguir un procedimiento análogo en dimensiones superiores parecen muy ad hoc y, lo que es más problemático, dan una distribución muy diferente para la columna complementada final que para las demás. ¿Será mi mejor opción algún tipo de muestreo de rechazo?

La motivación es generar matrices de probabilidad, siendo las filas estados del mundo y las columnas señales observadas. Cada señal es más probable en su estado correspondiente, pero no es necesariamente la señal más probable en ese estado.

Editar: La distribución exacta no importa, salvo que parezca natural. Sin embargo, me gustaría evitar que la diagonal también sea máxima para las filas.

Aquí está mi propio intento ah-hoc anterior.

gen.likelihood.matrix <- function(T){
    # T: Number of states

    mat <- matrix(runif(T^2), nrow=T) # Intialize with uniform variates
    diag(mat) <- 1 # To guarantee diag is max in column

    candidate.mat <- matrix(nrow=T, ncol=T)
    # Loop to check for constraint satisfaction
    for(k in 1:100){
        weights <- runif(T-1)
        weights <- weights / sum(weights)

        for(i in 1:T){
            # Rescale columns
            candidate.mat[i,-T] <- mat[i,-T] * weights
            # Replace last column with complement to meet row constraint
            candidate.mat[i,T] <- 1 - sum(candidate.mat[i,-T])
        }

        # Does the last column meet the maximality constraint?
        if(candidate.mat[T,T] > max(candidate.mat[-T,T])){
            # Shuffle to minimize the distortion from last column
            reorder <- sample(1:T,T)
            return(candidate.mat[reorder,reorder])
        }
    }

    # Recurse to get new initial values if necessary
    return(genLikeliMatrix(T))
}

Esto funciona en términos de darme matrices eficientemente, pero obtengo una extraña distribución bimodal en la diagonal debido a cómo manejo la última columna.

Answer 1

Esencialmente la misma idea de whuber, esta vez usando Dirichlets: las líneas de la matriz son Dirichlets independientes, con cada Dirichlet poniendo más masa en el elemento diagonal correspondiente.

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,])

Rechazar si x no satisface las restricciones. Si el porcentaje de rechazo es demasiado alto, aumente el valor de A0 . He aquí una muestra 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.D. Por favor, vectorice ese molesto for .

Answer 2

Bien, para avanzar en estos esfuerzos (pero con cierta desconfianza) ofrezco este enfoque: generar primero los elementos diagonales. Hacerlos grandes constantes . Generar todos los elementos no diagonales iid según cualquier distribución (no negativa) que se desee. Normalizar las filas. Compruebe la condición de columna-máxima. Repetir si se viola.

Haciendo que las constantes iniciales sean suficientemente grandes, el número esperado de repeticiones puede hacerse pequeño.

Está claro que los elementos diagonales son iid, los elementos no diagonales son iid, pero (por supuesto) las dos distribuciones difieren.

Aquí hay algo de código para jugar.

n <- 8                                     # Matrix dimension
y <- rep(1 + 3/sqrt(n),n)                  # A large constant compared to entries in x
x <- matrix(runif(n^2), ncol=n)            # Here, uniform distributions off diagonal
x[cbind(1:n,1:n)] <- y                     # (Paste in the diagonal)
z <- t(apply(x, 1, function(u) u/sum(u)))  # Normalize the rows
which(1:n != apply(z, 2, which.max))       # Find all columns violating the conditions

(Uno espera que integer(0) como salida; de lo contrario, saldrán los índices de las columnas cuyos máximos no son diagonales). He experimentado con n que van de 3 a 300.

Es instructivo trazar las columnas:

plot(z[1,], type="n")
apply(z, 2, function(u) lines(u, col=(256*runif(1))))

Answer 3

Primero, generaría una matriz $M$ tal que $M_{ij}$ se distribuye exponencialmente con media $\mu$ para todos $i \ne j$ o $0$ de lo contrario.

A continuación, calcula tu matriz \begin{align} X_{ij} &= \frac{e^{-M_{ij}}}{\sum_je^{-M_{ij}}} \end{align}

Para muestrear eficientemente de una distribución categórica (una lista de probabilidades con un total de 1), basta con construir un árbol de Huffman utilizando las probabilidades de cada resultado, llevando la cuenta de la probabilidad total en los descendientes a la izquierda de cada nodo. Tome una muestra de una distribución uniforme y encuentre el nodo de la hoja cuya CDF sea la menor de las mayores de su valor muestreado. (Esto tiene un coste amortizado de la entropía de su distribución categórica).

Answer 4

Primera generar al azar filas, a continuación, dividir cada fila por la suma de sus partes, de modo de llegar a la suma de uno de ellos. A continuación, el orden de las filas, de modo particular, en la primera fila que desea la fila que contiene el elemento en la primera posición, que es la máxima de la columna 1. Así, en la n-ésima fila que desea la fila que tiene el max de la n-ésima columna.

Answer 5

Sólo hay que generar cada fila a partir de un Distribución de Dirichlet entonces reordenar si es necesario para poner el valor más grande en la diagonal (o elegir los parámetros de tal manera que el más grande es probable que sea en la diagonal).