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]]
Hola, yo creo que deberías hacerlo vía cópulas. Mira en http://cran.r-project.org/web/packages/copula/index.html Es la manera habitual de generar distribuciones multivariantes, a partir de unas distribuciones marginales fijadas asignando una medida de correlación entre ellas.> 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 >jm~ _______________________________ J. Miguel Marin http://www.est.uc3m.es/jmmarin Dep. of Statistics University Carlos III of Madrid Spain (E.U.)
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------------ próxima parte ------------ Se ha borrado un adjunto en formato HTML... URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20101216/987debce/attachment-0002.html> ------------ próxima parte ------------ A non-text attachment was scrubbed... Name: Biomial2.pdf Type: application/pdf Size: 152811 bytes Desc: no disponible URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20101216/987debce/attachment-0001.pdf> ------------ próxima parte ------------ Se ha borrado un adjunto en formato HTML... URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20101216/987debce/attachment-0003.html>
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
>
>
>
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]]