1. En una guardería hay N niños y sólo una persona encargada de cuidarlos. De vez en cuando un bebé empieza a llorar exigiendo la atención de la niñera. Si la niñera está ocupada atendiendo a otro bebé, el nuevo bebé que llora debe esperar su turno. Si en el tiempo t un bebé está tranquilo entonces la probabilidad de que él empiece a llorar y exija ser atendido en el intervalo de tiempo (t,t+h) es igual a λh+o(h)λh+o(h). Si en el tiempo t un bebé está siendo atendido por la niñera entonces la probabilidad de que él se calme en el intervalo de tiempo (t,t+h] es igual a μh+o(h). Supóngase que Xt:=st número de bebés que están exigiendo ser atendidos en el tiempo tn.

(Ayuda: En este caso tenemos un proceso de nacimiento y muerte para el cual las tasas de nacimiento y muerte son, respectivamente, iguales a λi=(N−i)λ si 0≤i≤N y μi=μ para todo i≥1, es decir se trata de una cola de espera del tipo M/M/1/N.

# Función para calcular la distribución estacionaria en una cola M/M/1/N
mm1n_prob <- function(N, lambda, mu) {
  # Vector para almacenar probabilidades
  pi <- numeric(N + 1)
  
  # Estado base
  pi[1] <- 1  # Normalizamos con pi_0 = 1
  
  # Usamos la relación de balance de flujo
  for (i in 1:N) {
    pi[i + 1] <- pi[i] * ((N - i) * lambda) / mu
  }
  
  # Normalizamos dividiendo por la suma total de probabilidades
  pi <- pi / sum(pi)
  
  return(pi[1])  # Retorna la probabilidad de que X = 0
}

# Parámetros del primer caso
N1 <- 15
lambda1 <- 0.7
mu1 <- 0.2
p0_case1 <- mm1n_prob(N1, lambda1, mu1)
cat("Probabilidad de que haya 0 bebés en espera (Caso 1):", p0_case1, "\n")
## Probabilidad de que haya 0 bebés en espera (Caso 1): 2.082356e-19
# Parámetros del segundo caso
N2 <- 5
lambda2 <- 0.3
mu2 <- 0.5
p0_case2 <- mm1n_prob(N2, lambda2, mu2)
cat("Probabilidad de que haya 0 bebés en espera (Caso 2):", p0_case2, "\n")
## Probabilidad de que haya 0 bebés en espera (Caso 2): 0.0624438

La probabilidad para el caso 1 es de aproximadamente 0.02 La probabilidad para el caso 2 es de aproximadamente 0.06

  1. Escriba un programa de R o Python u otras para simular un proceso de Poisson no homogéneo con intensidad λ1(t)=sin⁡(t)+2 y λ2(t)=1/3*(sqrt(t/3)) ​en el intervalo[0,5000].
# Función para simular un proceso de Poisson no homogéneo
simulate_nhpp <- function(lambda, t_max) {
  t <- 0
  events <- numeric(0)
  
  while (t < t_max) {
    M <- max(lambda(t))
    t <- t - log(runif(1)) / M
    if (t > t_max) break
    if (runif(1) < lambda(t) / M) {
      events <- c(events, t)
    }
  }
  
  return(events)
}

# Definir las funciones de intensidad
lambda1 <- function(t) sin(t) + 2
lambda2 <- function(t) (1/3) * sqrt(t / 3)

# Simular los procesos
set.seed(123)
events1 <- simulate_nhpp(lambda1, 5000)
events2 <- simulate_nhpp(lambda2, 5000)

# Mostrar los primeros eventos
head(events1)
## [1] 0.6231314 0.9692001 0.9909314 1.2160154 1.4186320 1.4333976
head(events2)
## numeric(0)
  1. Haga un programa en R que simule al proceso de nacimiento puro {Xt,t≥0} con X(0)=1. Si caso 1: λ=iλ, para λ>0, y caso 2: λ=λi^(2), para λ>0, comparar estos dos casos.
# Función para simular un proceso de nacimiento puro
simulate_birth_process <- function(lambda_func, t_max, X0 = 1) {
  t <- 0
  X <- X0
  times <- c(t)
  states <- c(X)
  
  while (t < t_max) {
    rate <- lambda_func(X)
    if (rate == 0) break
    t <- t + rexp(1, rate)
    X <- X + 1
    times <- c(times, t)
    states <- c(states, X)
  }
  
  return(list(times = times, states = states))
}

# Definir las funciones de tasa
lambda1_func <- function(i) i * 0.1
lambda2_func <- function(i) i^(2) * 0.01

# Simular los procesos
set.seed(123)
process1 <- simulate_birth_process(lambda1_func, 100)
process2 <- simulate_birth_process(lambda2_func, 100)

# Graficar los resultados
plot(process1$times, process1$states, type = "s", col = "blue", xlab = "Tiempo", ylab = "Estado", main = "Proceso de Nacimiento Puro")
lines(process2$times, process2$states, type = "s", col = "red")
legend("topleft", legend = c("Caso 1", "Caso 2"), col = c("blue", "red"), lty = 1)

  1. Sea X=(Xt)t≥0​ una cadena de Markov de tiempo continuo con conjunto de estados S={0,1,2,3,4,5} y el generador Q. Hallar (si existe) la distribución estacionaria de X.
# Definir la matriz generadora Q según la imagen
Q <- matrix(c(
  -6, 6, 0, 0, 0, 0,
   4, -10, 6, 0, 0, 0,
   0, 8, -14, 6, 0, 0,
   0, 0, 8, -14, 6, 0,
   0, 0, 0, 8, -14, 6,
   0, 0, 0, 0, 8, -8
), nrow = 6, byrow = TRUE)

# Resolver el sistema para la distribución estacionaria
A <- t(Q)  # Transponer Q
A <- rbind(A, rep(1, 6))  # Agregar la ecuación de restricción sum(pi) = 1
b <- c(rep(0, 6), 1)  # Vector del lado derecho del sistema

pi <- qr.solve(A, b)  # Resolver el sistema

# Mostrar la distribución estacionaria
pi
## [1] 0.17933450 0.26900175 0.20175131 0.15131349 0.11348511 0.08511384