Muchas gracias a Olivier y Juan Miguel por su valiosa ayuda.

Feliz tarde,
Jorge Ivan Velez


2010/12/16 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@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 <>
>> 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:
>> 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@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>>
>>
>>
>>
>>
>

	[[alternative HTML version deleted]]

