# Modelo de difusión de Bernoulli-Laplace. Considere dos urnas, cada una de ellas conteniendo m bolitas; # b de ellas son negras y las restantes (2m − b) son rojas. Decimos que el sistema está en estado i si la # primera urna contiene exactamente i bolitas negras y m − i rojas, mientras que la segunda contiene # b − i negras y (m − b + i) bolitas rojas. Cada experimento consiste en elegir una bolita al azar de cada # urna e intercambiarlas. Sea X n el estado del sistema después de n intercambios # empezamos calculando la matriz de transicion m <- 50; b <- 40; transicion <- matrix(0,nrow=b+1,ncol=b+1); transicion[1,1]=(m-b)/m transicion[1,2]=b/m transicion[b+1,b+1]=(m-b)/m transicion[b+1,b]=b/m for (i in 1:(b-1)){ transicion[i+1,i] <- (i*(m-b+i))/(m^2); transicion[i+1,i+1] <- (i*(b-i))/(m^2) + (m-i)*(m-b+i)/(m^2); transicion[i+1,i+2] <- ((m-i)*(b-i))/(m^2); } #Escribo lo que sabemos que ya vale la invariante invariante <- rep(NA, b+1); for (i in 1:(b+1)){ invariante[i] <- choose(b,i-1) * choose(2*m-b, m-i+1) invariante[i] <- invariante[i]/choose(2*m,m) } invariante #busco t para que la t-potencia de la matriz diste menos de epsilon de la medida invariante epsilon = 0.00000001 t = 1 p = transicion while (max(abs(invariante-p[1,]))>epsilon){ t=t+1; p = p%*%transicion } print(t) #calculo la medida invariante resolviendo las ecuaciones de balance Id=diag(41) A=t(transicion - Id) A[41,]=rep(1,41) b=c(rep(0,40),1) x=t(solve(A,b)) #Verifico que todo lo que debería ser igual difiera menos que epsilon print(max(abs(x-p[1,]))epsilon){ Xn=Xn%*%transicion; n=n+1 } print(n)