Generación de variables aleatorias continuas en R

Blog_técnicas

Blog_simulación_vaCont

######GENERACIÓN DE VARIABLES ALEATORIAS CONTINUAS#########
### EN PARTICULAR, GENERAMOS UNA DISTRIBUCIÓN NORMAL Y A PARTIR DE ESTA, LA DITRIBUCIÓN CHI-CUADRADO
#cargar el package "randtoolbox"
 
#########Generar una va Normal a través del método de aceptación/rechazo#######
#queremos generar una va normal Z~N(0,1), el valor absoluto de Z tiene densidad de probabilidad f(x)=[2/(sqrt(2*phi))]*exp(-x^2/2) 0<x<infinito
#tenemos a g como la función de densidad exponencial con media 1 (lambda=1): g(x)=exp(-x) 0<x<infinito
#por tanto, f(x)/g(x)=sqrt(2/phi)*exp(x-x^2/2), y el valor máximo de f(x)/g(x) tiene lugar cuando x maximiza (x-x^2/2).
#calculando c, se llega a la expresión: c=max(f(g)/g(x))=f(1)/g(1)=sqrt(2exp(1)/phi) (con exp(1)=e)
#tenemos que: f(x)/(c*g(x))=exp(x-(x^2/2)-(1/2))=exp(-(x-1)^2/2)
#Sea Y la variable aleatoria exponencial generada a partir de g, seleccionamos X=Y si U<=exp((-(Y-1)^2/2)), donde U es un número aleatorio
 
#Asimismo, podemos modificar este último paso considerando su equivalente -log(U)>=((Y-1)^2)/2,
#y siendo que -log(U) también es una exponencial con parámetro 1, obtendríamos Y2>=(Y1-1)^2/2.
#Operando llegamos a la expresión: Y2-((Y1-1)^2)/2>0, y en caso que se cumpla la desigualdad, seleccionamos Y=Y2-((Y1-1)^2)/2.
#Finalmente, generamos un número aleatorio U y escogemos Z=Y1 si U<=1/2 o Z=-Y1 si U>1/2.
 
########algoritmo para generar una va exponencial con parámetro 1 y una va normal con media 0 y varianza 1, independiente de la va anterior.
#1)generar y1, va exponencial Y1 con lambda=1
#2)generar y2, va exponencial Y2 con lambda=1
#3)si y2-((y1-1)^2/2)>0, seleccionar y=y2-((y1-1)^2/2); en caso contrario volver al paso 1
#4)generar un número aleatorio u y seleccionar Z=y1 si u<=1/2 o Z=y2 si u>1/2
#4)repetir los pasos hasta alcanzar en número deseado de valores independientes de la variable Z.
#5)finalmente obtenemos Z -va con distrib. N(0,1)- e Y -va con distrib~Exp(1)-
 
#N representa la cantidad de valores que deseo generar
#se debe tomar un n suficientemente grande: n=N*2
normal<-function(dim=1,mod=2^31-1,mult=7^5,incr=0,echo=FALSE,N=100,lambda=0.5)
{
Y<-rep(NA,N)
for(i in 1:N)  #genero Y
  {
   while(is.na(Y[i])==TRUE)
   {
      y1<-exponencial(n=1)
      y2<-exponencial(n=1)
   if(is.na(Y[i])==TRUE & y2[1]-((y1[1]-1)^2/2)>0) {Y[i]<-(y2[1]-((y1[1]-1)^2/2))} else {Y[i]<-Y[i]}
   }
  }
Z<-rep(0,length(Y))
u<-congruRand(n=length(Y),dim,mod,mult,incr,echo)
for(k in 1:length(Y)) #genero Z
  {
  if(u[k]<=(1/2)){Z[k]<-Y[k]} else {Z[k]<--Y[k]}
  }
Z
}
 
#La distribución chi2 (chi-cuadrado) es la distribución de la suma de variables aleatorias normales al cuadrado. Sea Z_i números aleatorios i.i.d. normal(0,1), y seleccionando:
#chi2=sum(Z^2) suma de i=1 hasta n #Tenemos chi2 con distribución chi-cuadrado y n grados de libertad.
#run representa la cantidad de valores que deseo generar
#N representa el parámetro de chi2. Ej se denota: X^2(10)
chi_squared<-function(dim=1,mod=2^31-1,mult=7^5,incr=0,echo=FALSE,N=10,run=100)
{
chi2<-rep(0,run)
  for(i in 1:(run))
  {
    z<-normal(dim,mod,mult,incr,echo,N)
    chi2[i]<-sum(z^2)
  }
chi2
}
 
#####generación de variables aleatorias exponenciales
# Antes elegir setSeed(seed) sino se obtiene por default (y da mejores resultados)
exponencial<-function(n=10,dim=1,mod=2^31-1,mult=7^5,incr=0,echo=FALSE,lambda=1)   #generar va exponenciales con lambda=1
{
U<-congruRand(n,dim,mod,mult,incr,echo)
X<-NULL
X<-as.vector(rep(0,length(U)))
    for(i in 1:length(X))
    {
      X[i]<-(-log(U[i])*(1/lambda))
    }
X
}
 
############GENERACIÓN UNA LA VARIABLE ALEATORIA GAMMA(N,lambda)#############
# cargar library("randtoolbox")
 
#El objetivo es generar una variable aleatoria (va) Z con distribución gamma(n,lambda).
#Para la función de distribución de esta va, no es posible realizar una adecuada transformación inversa de esta función,
#por ello utilizaremos el hecho de que la va Z~gamma(n,lambda) puede plantearse como una suma de n va exponenciales independientes de parámetro lambda, X~Exp(lambda).
#Por tanto podemos utilizar estas va exponenciales para generar Z: Z=sum(Xi[1:n])
#En particular, podemos generar a Z~gamma(n,lambda) a partir de la generación de n números aleatorios U=u1…un,
#y seleccionamos Z=-(1/lambda)*log(u1)-…-(1/lambda)*log(un)=-(1/lambda)*log(u1…un)
 
#Algoritmo:
#1) generamos un vector de n números pseudoaleatorios uniformes U=u1…un
#2) elegimos Z=-(1/lambda)*log(prod(U[1:n]))
 
gamma<-function(n=10,dim=1,mod=2^31-1,mult=7^5,incr=0,echo=FALSE,lambda=1,N=10)
{
Z<-rep(0,N)
  for(i in 1:N)
  {
    U<-congruRand(n,dim,mod,mult,incr,echo)
    X<-NULL
    X<-as.vector(rep(0,length(U)))
    X<--(log(prod(U))*(1/lambda))
    Z[i]<-X
  }
Z
}
 
#########TÉCNICA DE COMPOSICIÓN DISCRETA##########
#El objetivo es simular el valor de una variable aleatoria (va) cuya función masa de probabilidad (fmp) está dada por:
#P[X=j]=(alpha*pj)+((1-alpha)*qj) donde 0<alpha<1 , j>=0
#siendo que se define X como X1 (con probabilidad alpha) y X2 (con probabilidad (1-alpha)), X1 y X2 son va
#Para estas variables, cuya distribución es una mixtura o composición de otras, la transformada inversa y técnica de aceptación-rechazo, no suelen ser apropiadas.
#Por ello, utilizaremos la técnica de composición: generamos el valor de cada va generando un Nºaleat U y luego obteniendo el valor de x1 si U<alpha y el de X2 si U>alpha
 
#EJ: pj=P[X=j]=   i)0.05  si j=1:5 o ii) 0.15  si j=6:10
#podemos observar que pj=0.5*p1j y 0.5*p2j, donde
# p1j=0.1 para j=1:10 y
# p2j=  i) 0 si j=1:5 o ii) 0.2 si j=6:10
#####algoritmo####
#1) generar un número aleatorio uniforme u1
#2) generar un número aleatorio uniforme u2
#3) si u1<=0.05 elegimos X=as.integer(10*u2)+1; en caso contrario elegimos X=as.integer(5*u2)+6
#se elige antes la semilla con setSeed(seed)
composition.aprox<-function(n=10,dim=1,mod=64,mult=13,incr=0,echo=FALSE,alfa=0.5)
{
out<-NULL
  X<-NULL
  X<-as.vector(rep(0,n))
        for(i in 1:n)
        {
          U1<-congruRand(n=1,dim,mod,mult,incr,echo)
          U2<-congruRand(n=1,dim,mod,mult,incr,echo)
          if(U1<=alfa) {X[i]<-as.integer(10*U2)+1} else {X[i]<-as.integer(5*U2)+6}
        }
X
}
 
#########VARIABLES ALEATORIAS DE POISSON CON PARÁMETRO LAMBDA#############
# cargar library("randtoolbox")
 
############GENERACIÓN DE UNA VA CON DISTRIBUCIÓN POISSON (~P(lambda))############
#####generación de una distribución de Poisson(lambda) a través del método de transformación inversa
#Algoritmo:
#1) generar un número aleatorio U=U1 .#si n=1 genero 1 número, si n>1 genero un vector
#2) i=0, p=exp(-lambda), F=p
#3) si U<F, elegimos X=i
#4) p=p*lambda/(i+1),    F=F+p,        i=i+1
#5) volvemos al paso 3
 
poisson.distr<-function(n=1,dim=1,mod=2^31-1,mult=7^5,incr=0,echo=FALSE,lambda=5)
{
U<-congruRand(n,dim,mod,mult,incr,echo)
X<-rep(NA,n)
for(j in 1:n)
  {
  i=0
  p=exp(-lambda)
  F=p
  if(U[j]<F){X[j]=i}
  else {
  while(U[j]>=F)
      {
      p=lambda*p/(i+1)
      F=F+p
      i=i+1
      }
      X[j]=i}
  }
X
}
 
 
########GENERACIÓN DE LAS PRIMERAS T UNIDADES DE UN PROCESO POISSON DE PARÁMETRO lambda#########
#generaremos sucesivamente los intervalos de tiempo, y finalizaremos cuando su suma sea mayor a T.
#para éste tipo de proceso, los tiempos entre eventos sucesivos son va exponenciales independientes con parámetro lambda
#los valores generados en los primeros n eventos temporales son sum(Xi[1:n]), con j=1…n.
 
#Algoritmo: para generar todos los eventos temporales (I) que ocurren en (0,T) de un proceso de Poisson con parámetro lambda.  S(I) representa el evento más reciente.
#1) Comenzamos con t=0 e I=0,
#2) generamos el número aleatorio u
#3) Fijamos t=t-(1/lambda)*log(u) y si t>T detenemos la simulación
#4) seleccionamos I=I+1, S(I)=t y volvemos al paso 2
#El valor obtenido de I representará el número de eventos que ocurren en el tiempo T, y los valores S(1)…S(I) serán los eventos temporales en orden creciente.
 
#se debe fijar setSeed(seed) antes de comenzar
poisson.proc<-function(n=1,dim=1,mod=2^31-1,mult=7^5,incr=0,echo=FALSE,lambda=.1,Tiempo=100)   #Tiempo=T
{
ti=0  #comenzamos en t=0
In=0  #representará el número de eventos que ocurren en el tiempo T (se renueva en cada iteración)
S<-NULL
S_In<-NULL   #será el vector de los eventos temporales en orden creciente (se renueva en cada iteración)
    while(ti<=Tiempo)    #detenemos la simulación si t>T
    {
        U<-congruRand(n,dim,mod,mult,incr,echo)  #generamos el número aleatorio u
        ti=ti-(1/lambda)*log(U)
        In=In+1
        S_In=ti
        S<-c(S,S_In)
    }
out<-NULL
out<-list("número de eventos en T"=In,"Tiempos de cada evento"=S)
out
}
#observación: el tiempo del último evento es >T, entonces debería eliminarlo,no? S<-S[-length(S)] e In<--1