Distribución de probabilidad de programación en R

Question

Estoy tratando de simular la respuesta a este problema en R.

Del noroeste Minnesota es propensa a las inundaciones en la primavera. Supongamos que el 20% de todos los agricultores están asegurado contra daños por inundaciones. Cuatro granjeros son seleccionados al azar. ¿Cuál es la probabilidad de que al menos dos agricultores tienen seguro contra inundaciones?

Gracias por cualquier ayuda de cualquier persona.

Answer 1

Esta respuesta le guiará a través de los pasos que son comunes a la mayoría de las simulaciones estocásticas, que muestra cómo crear el código en pequeño, sencillo, de fácil probado en trozos en cualquier plataforma de software. El proceso se ilustra con R código. Puede ejecutar los fragmentos de código a medida que avanza a ver lo que producen.

FWIW, he aquí un compacto, rápido y sucio solución que comprime los cálculos en una línea. Imprime una estimación basada en una simulación con n de iteraciones. (Figura alrededor de cuatro millones de iteraciones por segundo de tiempo de CPU.)

n <- 1e6; mean(rowSums(matrix(runif(4*n) <= 0.20, nrow=n)) >= 2)

---

Paso 1: Simular un agricultor.

Crear un cuadro con las entradas, una por cada agricultor. En las entradas se escriben "el asegurado" y "no asegurado." La proporción de "el asegurado" en la caja debe ser de 20%. (Para una descripción detallada de las entradas-en-un-modelo de cuadro de variables aleatorias, ver Lo que se entiende por una variable aleatoria

En casi todos los sistema de software que apoya la generación de números aleatorios, hay una función para dibujar un billete de salida de una caja muy grande (y reemplazarlo después de que el billete se observa). Estas entradas tienen valores de punto flotante de $0$ a sólo un poco menos de $1$ escrito en ellos. Para cualquier intervalo de $0 \le a \le b \lt 1$, la proporción de billetes con valores de entre $a$$b$$b-a$. Usted puede aprovechar esta función, mediante la escritura de "el asegurado" en el 20% de todos los números entre el$0$$1$. Por ejemplo, podría escribir "asegurado" en todos los números entre el$0$$20/100$.

Este es el Uniforme de distribución (entre $0$$1$). En R esta función se llama runif. Un único sorteo de una entrada de este cuadro podría ser programado como

label <- ifelse(runif(1) <= 20/100, "insured", "not insured")

Es más rápido y más conveniente prescindir de la etiqueta, sin embargo, y reemplazarlo con el número de $1$ (reservando $0$ para todas las otras entradas). Esto le da la más simple R código

label.indicator <- runif(1) <= 20/100

porque en R (como en muchos sistemas) de un verdadero valor es igual a $1$ y un valor falso con $0$.

Paso 2: Simular cuatro agricultores.

Sólo las entradas para el sorteo de la caja de forma independiente. Eso se hace con un bucle. En R el bucle se realiza (de manera muy eficiente, detrás de las escenas) cuando la solicitud de varios valores de runif. Lo que es importante saber es que estos valores son (pseudo) independiente: parecen no dependen una de la otra. Por lo tanto,

label.indicators <- runif(4) <= 20/100

se simula el dibujo de las entradas de cada cuatro agricultores de este cuadro. Se produce una serie de cuatro números del conjunto {FALSE, TRUE} o equivalentemente $\{0,1\}$.

Paso 3: Calcular la estadística (el agricultor cuenta).

El número de agricultores que en cualquier grupo es encontrado por la suma de sus indicadores.

farmer.count <- sum(label.indicators)

Esto es debido a que cada asegurado agricultor contribuye a $1$ a la suma y cada uno de los asegurados agricultor contribuye a $0$. La suma meramente cuenta del asegurado los agricultores. La eficiencia de la utilización de indicadores, en lugar de las etiquetas, es evidente aquí.

Paso 4: Repita muchas veces.

Esto es un bucle en la mayoría de las plataformas. En algunos, como R, es más rápido para dibujar un montón de billetes y agruparlos en grupos de cuatro. (Esto es debido a que a menudo es casi tan rápido a generar muchos valores aleatorios como para generar una sola: hay menos gastos generales involucrados.) Cada grupo representa una iteración de la simulación. El código siguiente pone de los cuatro boletos para cada iteración en las filas de una matriz y, a continuación, suma de cada fila (como en el Paso 3).

n.sim <- 1e6   # Number of iterations
n.farmers <- 4 # Number of farmers per iteration
simulation <- rowSums(matrix(runif(n.sim * n.farmers) <= 20/100, nrow=n.sim))

Paso 5: Post-proceso.

La oportunidad de observar dos o más asegurados de los agricultores puede ser estimado como la proporción de iteraciones en las que dos o más asegurados los agricultores fueron encontradas en la muestra. Como antes, esto se encuentra eficazmente por las pruebas y sumando (o prueba y el promedio):

estimate <- mean(simulation >= 2)

Paso 6: Evaluar.

Sólo han realizado un equipo en el experimento. Como cualquier otro experimento, los resultados son variables, por lo que debe proporcionar un error estándar para el resultado. En este experimento Binomial, el error estándar de la estimación $\hat p$ ($n$ iteraciones) es

$$\text{se}(\hat p) = \sqrt{\hat p\left(1 - \hat p\right) / n}.$$

Compute this and print out the results:

se.estimate <- sqrt(estimate * (1-estimate) / n.sim)
print(c(Estimate=estimate, SE=se.estimate), digits=2)

When the random seed is set to $17$ (ver el código completo más abajo) y un millón de iteraciones, el resultado es

Estimate       SE 
 0.18127  0.00039 

---

Aquí está la solución completa. Está estructurado para permitir la variación de los parámetros de manera que, mediante la repetición de la misma (el cálculo se toma menos de un segundo), se puede estudiar la forma en que los resultados dependen de los parámetros para obtener más sentido intuitivo de lo que está sucediendo.

#
# Specify the problem.
#
p <- 20/100    # Chance of insurance
k <- 2         # Minimum number of insured farmers
n.farmers <- 4 # Number of farmers per iteration
n.sim <- 1e6   # Number of iterations
#
# Simulate.
#
set.seed(17)   # Optional: provides a reproducible result
simulation <- rowSums(matrix(runif(n.sim * n.farmers) <= p, nrow=n.sim))
#
# Post-process and report.
#
estimate <- mean(simulation >= k)
se.estimate <- sqrt(estimate * (1-estimate) / n.sim)
print(c(Estimate=estimate, SE=se.estimate), digits=2)

Answer 2

Aquí está una posible solución a tu pregunta

prop=0
T=1e6
for (t in 1:T){
  insrd=sample(0:1,4,rep=TRUE,prob=c(8,2))
  prop=prop+(sum(insrd)>1)}
print(prop/T)

con respuesta 0.1809. Pero usted puede calcular el valor exacto [0.1808] de esta probabilidad por darse cuenta de que el sorteo de cuatro agricultores es una variable de aleatoria de B(4,0.2) Binomial cuando el número de agricultores asegurados fuera de los cuatro.