Foro de debate

matrices por bloques

matrices por bloques

by Irene Castro Conde -
Number of replies: 6
Hola, necesito introducir una matriz por bloques de gran dimensión y no se cómo hacerlo de forma sencilla. Hay alguna función o método que pueda utilizar. También quisiera encontrar una forma de generar una variable normal multidimensional de forma más eficiente que con mvrnorm. Gracias.
In reply to Irene Castro Conde

Re: matrices por bloques

by Manuel Muñoz Márquez -
Hola Irene:

Puedes crear la matriz y luego ir asignándola por bloques, te pongo un ejemplo:
> m<-matrix(0, nrow=3, ncol=3)
> m
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
[3,] 0 0 0
> m[1:2, 1:2] <- diag(1)
> m[3,3]<-2
> m
[,1] [,2] [,3]
[1,] 1 1 0
[2,] 1 1 0
[3,] 0 0 2

Con respecto a lo de mvrnorm, ¿qué problema tienes? ¿qué te hace pensar que existe un método más eficiente que el que usa el programador de dicha función?

Saludos.
In reply to Manuel Muñoz Márquez

Re: matrices por bloques

by Irene Castro Conde -
Gracias por la ayuda pero necesito otra forma más directa de introducir la matriz.Necesito un código que me sirva para todos los casos, es decir, si necesito 2 ,3,6,.. bloques. Porque de la forma que me dices solo sirve para un caso particular. Por ejemplo para 2 bloques lo he programado así:

ro<-0.3
s<-1000
k<-500#dimensión de los bloques
sigma2<-matrix(0,nrow=s,ncol=s)
sigma2[1:k,1:k]<-diag(1-ro,k)+matrix(ro,k,k)
sigma2[k+1:k,k+1:k]<-diag(1-ro,k)+matrix(ro,k,k)

Pero necesito hacer simulaciones cambiando el número de bloques luego tendría que hacer un código para cada caso y lo que busco es algo más eficiente.

En cuanto a la normal, el generador mvrnorm tarda demasiado tiempo en ejecutarse en mis simulaciones. creo que utilizando la descomposición de Cholesky puede ser más eficiente.

In reply to Irene Castro Conde

Re: matrices por bloques

by Manuel Muñoz Márquez -
Hola:

No veo el problema de crear una función que haga lo que quieres.

mymatrix <- function(s, k, ro)
{
k<-500#dimensión de los bloques
sigma2<-matrix(0,nrow=s,ncol=s)
sigma2[1:k,1:k]<-diag(1-ro,k)+matrix(ro,k,k)
sigma2[k+1:k,k+1:k]<-diag(1-ro,k)+matrix(ro,k,k)
return sigma2
}
mymatrix(1000, 500, .3)

Si quieres más bloques pues otro parámetro,...

Tal vez si nos cuentas un poco más de tu problema podamos ayudarte un poco mejor.

¿Qué paquete está usando? La función
rmvnorm(n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)),
method=c("eigen", "svd", "chol"))
del paquete mvtnorm admite como método el de Cholesky.

Saludos.
In reply to Manuel Muñoz Márquez

Re: matrices por bloques

by Irene Castro Conde -
Hola,

estaba utilizando mvrnorm, con el que me has dicho ya me funciona mucho mejor. Gracias.

Te cuento un poco mi historia. Formo parte de una investigación de estadística sobre contrastes múltiples y me han pasado un código donde se realiza una simulación y mi tarea es optimizarlo. Entonces, un punto que interesa es ir variando la forma de la matriz sigma en cuanto al número de bloques para ver distintos resultados. Tal como me lo han dado, cada matriz estña definida por separado y se utiliza en cada simulación la que se quiera. Lo que me interesa es hacer una función del estilo que me has escrito tú pero en la que se pueda variar el número de bloques y no veo la forma de hacerlo. Porque en tu función si en vez de k=500 pones k=200, que se corresponde con 5 bloques, la función no contruye 5 bloques sino 2. No se si me explico.He pensado que igual se puede hacer añandiendo un bucle for dentro de la función. Decir que tengo muy poca experiencia programando en R y he encontrado este foro buscando un poco de ayuda.

Saludos.

In reply to Irene Castro Conde

Re: matrices por bloques

by Manuel Muñoz Márquez -
Hola:

Los bucles en R son muy poco eficientes y cuando es posible es conveniente evitarlos usando vectorización, estoy convencido de que lo que pretendes se puede hacer sin bucles, pero bueno, podemos empezar con bucles.

Si no he entendido mal tu problema se resuelve con la función:
mymatrix <- function(s, k, ro)
{
sigma2<-matrix(0,nrow=s,ncol=s)
for (i in seq(1, s, k))
{
j <- i +k -1
sigma2[i:j,i:j]<-diag(1-ro,k)+matrix(ro,k,k)
}
sigma2
}
> mymatrix(6, 2, .3)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.0 0.3 0.0 0.0 0.0 0.0
[2,] 0.3 1.0 0.0 0.0 0.0 0.0
[3,] 0.0 0.0 1.0 0.3 0.0 0.0
[4,] 0.0 0.0 0.3 1.0 0.0 0.0
[5,] 0.0 0.0 0.0 0.0 1.0 0.3
[6,] 0.0 0.0 0.0 0.0 0.3 1.0
> mymatrix(6, 3, .3)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.0 0.3 0.3 0.0 0.0 0.0
[2,] 0.3 1.0 0.3 0.0 0.0 0.0
[3,] 0.3 0.3 1.0 0.0 0.0 0.0
[4,] 0.0 0.0 0.0 1.0 0.3 0.3
[5,] 0.0 0.0 0.0 0.3 1.0 0.3
[6,] 0.0 0.0 0.0 0.3 0.3 1.0

Por cierto, tu matriz resultante tiene valor 1 en la diagonal ¿es eso correcto?

Saludos.
In reply to Manuel Muñoz Márquez

Re: matrices por bloques

by Irene Castro Conde -

Hola, muchas gracias por la ayuda, esto es exactamente lo que andaba buscando, lo que me pasaba es que no caía en usar una secuencia. En cuanto a la vectorización, ya me han comentado que es más eficiente pero no se muy bien como usarla en este caso. En el siguiente código, por ejemplo, se podría eliminar el bucle como indico a continuación?:

for (j in 1:M) {
for (k in 1:s){
t12[j,k]=t.test(xx1[,k],xx2[,k],var.equal=T)$statistic}
}

for (j in 1:M) {
t12[j,1:s]=t.test(xx1[,1:s],xx2[,1:s],var.equal=T)$statistic
}

Saludos.