[R-es] Generacion de binomiales correlacionadas

[Olivier Nuñez]

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
>
>
>