[R-es] Generacion de binomiales correlacionadas

Jorge,

redacté a partir del articulo de Biswas y Hwang (2002) y de sus  
notaciones,
el código (mejorable) de una función para generar N realizaciones de  
dos variables binomiales Y_i  ~ Binomial(n_i, p_i)  ( i =1,2) y con  
correlación rho:


r2Binom <- function(N=100,n1=10,n2=10,p1=.5,p2=.5,rho=0){
m=min(n1,n2)
M=max(n1,n2)
min.alpha=max(-(1-p2)/(1+p1-p2),-p2/(1+p2-p1))
if(p1 > p2) max.alpha=p2/(p1-p2) else max.alpha=(1-p2)/(p2-p1+1e-30)
c=sqrt(M/m)*sqrt(p2*(1-p2)/(p1*(1-p1)))
rho.min=max(round(min.alpha/(1+min.alpha)/c,2),-1)
rho.max=min(round(max.alpha/(1+max.alpha)/c,2),1)
if(rho>rho.max |rho<rho.min) stop("la correlación rho ha de ser  
incluida entre ",rho.min," y ",rho.max,".")
z=rho*sqrt(M/m)*sqrt(p2*(1-p2)/p1/(1-p1))
alpha=z/(1-z)
Y1j=matrix(rbinom(M*N,1,p1),M,N)
p2y=(p2+alpha*(p2-p1)+alpha*Y1j)/(1+alpha)
Y2j=apply(p2y,2,FUN=function(x) rbinom(M,1,x))
Y1=apply(Y1j[1:n1,],2,sum)
Y2=apply(Y2j[1:n2,],2,sum)
return(cbind(Y1,Y2))
}



 > Y=r2Binom(2000,n1=10,n2=8,p1=.3,p2=.5,rho=-.4)
 > summary(Y)
        Y1              Y2
  Min.   :0.000   Min.   :0.000
  1st Qu.:2.000   1st Qu.:3.000
  Median :3.000   Median :4.000
  Mean   :3.010   Mean   :4.029
  3rd Qu.:4.000   3rd Qu.:5.000
  Max.   :7.000   Max.   :8.000
 > cor(Y)
            Y1         Y2
Y1  1.0000000 -0.3881115
Y2 -0.3881115  1.0000000
Un saludo. Olivier


PD: Hay que tomar en cuenta que hay restricciones sobre el rango de  
valores de rho dados p1 y p2.
Lo cual es bastante lógico. Cuentame si te funciona.

____________________________________

Olivier G. Nuñez
Email: onunez en iberstat.es
Tel : +34 663 03 69 09
Web: http://www.iberstat.es

____________________________________




El 16/12/2010, a las 3:51, Jorge Ivan Velez escribió:

Gracias Olivier. Ese fue uno de los primeros papers que encontre y  
sobre el que comence a trabajar  :)  Al parecer no existe (bueno,  
al menos no la he encontrado aun) una manera simple de generar, via  
R, directamente lo que necesito.

Feliz noche,
Jorge


2010/12/15 Olivier Nuñez <onunez en iberstat.es>
Es un problema delicado, y como alternativa a lo que JuanMi propone,
puedes también utilizar la metodología expuesta en el articulo  
adjunto.
Consiste básicamente en sumar combinaciones de Bernouilli.
Un saludo. Olivier

--  

____________________________________

Olivier G. Nuñez
Email: onunez en iberstat.es
Tel : +34 663 03 69 09
Web: http://www.iberstat.es

____________________________________




El 15/12/2010, a las 23:08, Jorge Ivan Velez escribió:

Buenas tardes,

Estoy interesado en generar observaciones de una distribucion  
binomial
bivariada en la que hay _cierto_ grado de correlacion (denotemoslo  
rho).
Podria por favor alguien indicarme como hacerlo en R?

Este es el contexto. Supongamos que se tienen dos experimentos en  
los que la
variable respuesta _sigue_ una distribucion binomial, i.e., X_i
~Binomial(n_i, p_i), i=1,2 y que, por ahora, n_i y p_i son  
conocidos. El
interes principal es calcular T = P(X_1=1, X_2=1).  Si las X_i's  
fueran
independientes (caso 1), seria suficiente generar binomiales con los
parametros correspondientes a cada tipo de experimento (via rbinom) y
calcular lo que se necesita usando table(). Ahora, si X_i's *no* son
independientes (caso 2), el calculo de T no puede hacerse de la  
misma forma.
Observe que independiente del modelo (caso) bajo el cual estemos  
trabajando,
T _podria_ calcularse como una funcion de rho, n_i y p_i.  
Desafortunadamente
no tengo claro como _modelar_ el segundo caso utilizando R.

Alguna sugerencia?

Muchas gracias,
      Jorge Ivan Velez

[[alternative HTML version deleted]]

_______________________________________________
R-help-es mailing list
R-help-es en r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es