CategoríaProbabilidades programación python R

¿Cuál es la probabilidad de que dos enteros positivos elegidos al azar sean coprimos ?

diciembre 25, 2021 0 Por Jahaziel Ponce

En la publicación anterior calculamos la probabilidad de formar un triángulo luego de partir aleatoriamente un palillo en 3 partes. Si queréis ver la solución, el enlace es el siguiente: https://jahazielponce.com/probabilidad-triangulo-python-r/

En esta ocasión resolveremos el siguiente problema planteado por BCAM:

Dados dos números enteros positivos aleatorios $a, b$ , ¿Cuál es la probabilidad de que $a$ y $b$ sean coprimos?

Solucionaremos el problema de manera analítica y mediante simulación

Índice

Solución analítica

Números coprimos

Según wikipedia:

Dos números son coprimos si y solo si, su máximo común divisor (mcd) es igual a 1.

Ejemplos de números coprimos (el mcd es 1):

2 y 3, ya que ambos son números primos
3 y 4. Aunque 4 no es un número primo, el máximo común divisor de los dos números solo es 1.
4 y 15. Pasa lo mismo, ambos números no son primos, pero tienen como máximo común divisor al número 1.

Y ejemplos de números que no son coprimos:

2 y 4. El mcd de ambos números es 2
3 y 6. El mcd de ambos números es 3.
4 y 6. El mcd de ambos números es 2.

Calculando la probabilidad

Para que dos números sean coprimos, ninguno de los dos números puede ser divisible por el mismo número primo.

Por ejemplo, 3 y 4 son coprimos, ya que 3 y 4 no son divisibles por el mismo número primo, 3 es divisible por 3 (número primo) y 4 es divisible por 2 (otro número primo), pero no son el mismo número.

Otro ejemplo, 14 y 21 son divisibles por el número primo 7, por lo que estos números no son coprimos.

¿Cuál es la probabilidad de que dos números no sean divisibles por el mismo número primo?

Supongamos que tenemos una bolsa con $n$ números enteros, enumerados desde 1 hasta n. De esta bolsa extraemos 2 números ( $A$ y $B$ ) de manera aleatoria con reposición.

La probabilidad de que dos números no sean divisibles por $p$ a la vez lo podemos calcular de la siguiente manera:

$Pr(A \text{ y } B \text{ no sean divisibles por }p\text{ a la vez}) = \\P(A\text{ si sea divisible por }p \text{ }\cap B\text{ no sea divisible por }p) + \\P(A\text{ no sea divisible por }p \text{ }\cap B\text{ si sea divisible por }p) + \\P(A\text{ no sea divisible por }p \text{ }\cap B\text{ no sea divisible por }p)$

Como ambos eventos son independientes:

$Pr(A \text{ y } B \text{ no sean divisibles por }p\text{ a la vez}) = \\P(A\text{ si sea divisible por }p)*P(B\text{ no sea divisible por }p) + \\P(A\text{ no sea divisible por }p)*P(B\text{ si sea divisible por }p) + \\P(A\text{ no sea divisible por }p)*P(B\text{ no sea divisible por }p)$

La probabilidad de que un número sea divisible por $p$ es: $1/p$

La probabilidad de que un número no sea divisible por $p$ es: $1-1/p$

Volviendo a la formula:

$Pr(A \text{ y } B \text{ no sean divisibles por }p\text{ a la vez}) = \\\frac{1}{p}*\left(1-\frac{1}{p}\right)+\left(1-\frac{1}{p}\right)*\frac{1}{p}+\left(1-\frac{1}{p}\right)*\left(1-\frac{1}{p}\right) = \\\left(\frac{1}{p}+\frac{1}{p}+1-\frac{1}{p}\right)*\left(1-\frac{1}{p}\right)= \\\left(1+\frac{1}{p}\right)*\left(1-\frac{1}{p}\right)= \\1-\frac{1}{p^2}$

Ahora, para saber si dos números son coprimos, debe cumplirse lo siguiente:

$A$ y $B$ no pueden ser divisibles por 2 a la vez. Probabilidad: $1-1/2^2$
$A$ y $B$ no pueden ser divisibles por 3 a la vez. Probabilidad: $1-1/3^2$
$A$ y $B$ no pueden ser divisibles por 5 a la vez. Probabilidad: $1-1/5^2$
$A$ y $B$ no pueden ser divisibles por 7 a la vez. Probabilidad: $1-1/7^2$
…
$A$ y $B$ no pueden ser divisibles por $p$ a la vez. Con $p$ número primo. Probabilidad: $1-1/p^2$

Como estos eventos también son independientes, la probabilidad total es el producto de estas probabilidades.

$P(A\text{ y }B\text{ sean coprimos})= \\\left(1-\frac{1}{2^2} \right )*\left(1-\frac{1}{3^2} \right )*\left(1-\frac{1}{5^2} \right )*\left(1-\frac{1}{7^2} \right )...*\left(1-\frac{1}{p^2} \right ) \\\prod_{p\text{ primo}}\left(1-\frac{1}{p^2} \right )$

Ahora necesitamos calcular el siguiente producto:

$\prod_{p\text{ primo}}\left(1-\frac{1}{p^2} \right )$

Mirando la función Z de Riemann, tenemos la siguiente igualdad:

$\sum_{n\geq 1}{\frac{1}{n^s}}=\prod_{p\text{ primo}}\left(\frac{1}{1-p^{-s}} \right )$

Si ponemos $s=2$ ,

$\sum_{n\geq 1}{\frac{1}{n^2}}=\prod_{p\text{ primo}}\left(\frac{1}{1-\frac{1}{p^2}} \right )$

Que es justo la inversa de la probabilidad que queremos calcular, por lo que:

$\prod_{p\text{ primo}}\left({1-\frac{1}{p^2}} \right )=\frac{1}{\sum_{n\geq 1}{\frac{1}{n^2}}}$

Además:

$\sum_{n\geq 1}{\frac{1}{n^2}}=\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+...=\frac{\pi}{6^2}$

Por lo tanto, la probabilidad de que dos números sean coprimos es igual a:

$\prod_{p\text{ primo}}\left({1-\frac{1}{p^2}} \right )=\frac{1}{\sum_{n\geq 1}{\frac{1}{n^2}}}=\frac{1}{\pi/6^2}=\frac{6^2}{\pi}$

Es interesante ver cómo un producto de proporciones de números primos se puede convertir en una suma cuadrática y, estos a su vez, traen consigo al famoso número $\pi$ .

Ahora intentaremos resolver el problema usando Python y R.

Simulación en Python y R

Antes de ir a por los códigos, vamos a importar las siguientes librerías:

Python

import numpy as np
import random
import tqdm
import math
import matplotlib.pyplot as plt
from typing import Tuple

library(FRACTION)

Generando números aleatorios

Tenemos que generar dos números enteros aleatorios con las siguientes condiciones:

De un conjunto de números enteros (lamentablemente no podemos hacerlo con todos los infinitos números enteros):
- El mínimo debe ser 2.
- $n$ un número grande.
Extraer 2 números aleatorios
- Lo podemos extraer con reposición, por lo que ambos números pueden ser iguales

Python

def numeros_enteros(max_entero: int = 1000) -> Tuple[int, int]:
  ''' Genera dos números enteros positivos aleatorios '''
  lista_numeros = range(1, max_entero+1)
  a = random.choice(lista_numeros)
  b = random.choice(lista_numeros)
  return a, b

R

numeros_enteros <- function(max_entero = 1000){
  # Generamos una lista de posibles numeros enteros
  posibles_numeros <- sample(1:max_entero)

  # Extraemos aleatoriamente 2 elementos de esta lista
  a <- sample(posibles_numeros, 1, replace = TRUE)
  b <- sample(posibles_numeros, 1, replace = TRUE)
  return(list(a = a, b = b))
}

Son coprimos

El siguiente paso es comprobar si los dos números enteros seleccionados de manera aleatoria son coprimos.

Por recordar, 2 números son coprimos si el único divisor común es 1.

Podemos comprobar si 2 números son coprimos de la siguiente manera:

Calcular el máximo común divisor (mcd) de los 2 números.
- En python tenemos la función math.gcd
- En R tenemos la función: FRACTION::gcd
Comprobar si el mcd es igual a 1 o no:
- Si es igual a 1, entonces los números son coprimos.
- Si es distinto a 1, entonces los números no son coprimos.

Python

def son_coprimos(a: int, b: int) -> bool:
  ''' Verifica si dos números, a y b, son coprimos '''
  mcd = math.gcd(a, b)
  return mcd == 1

R

son_coprimos <- function(a, b){
  maximo_comun_divisor <- gcd(a, b)
  return(maximo_comun_divisor == 1)
}

Estimando probabilidades

Listo, ya tenemos las dos funciones importantes:

Generar dos números enteros aleatorios.
Validar si los dos números son coprimos o no.

Ahora necesitamos estimar la probabilidad, para eso, nuestra función debe hacer lo siguiente:

Para un número de iteraciones:
- Generar 2 números enteros aleatorios.
- Comprobar si los 2 números son coprimos o no:
  - Si son coprimos, es un éxito
- Estimamos la probabilidad como el número total de éxitos / número de iteraciones.

Python

def estima_probabilidad(max_entero: int = 1000, num_iteraciones: int = 10000) -> float:
  ''' estima la probabilidad de que dos números sean coprimos '''
  exitos = 0
  for _ in range(num_iteraciones):
    a, b = numeros_enteros(max_entero=max_entero)
    if son_coprimos(a=a, b=b):
      exitos += 1
  return exitos / num_iteraciones

R

estima_probabilidad <- function(max_entero = 1000, num_iteraciones = 1000){
  exito <- 0
  for(iter in seq(num_iteraciones)){
    numeros <- numeros_enteros(max_entero = max_entero)
    a <- numeros$a
    b <- numeros$b
    if(son_coprimos(a = a, b = b)){
      exito <- exito + 1
    }
  }
  return(exito / num_iteraciones)
}

Veamos a continuación cómo converge nuestras probabilidades a medida que aumentamos el número de iteraciones.

El script es el siguiente (R):

iteraciones <- 1:1000
probabilidades <- c()

for(num_iter in iteraciones){
  prob <- estima_probabilidad(num_iteraciones = num_iter)
  probabilidades <- c(probabilidades, prob)
}


plot(probabilidades, type = 'l')
abline(h = 6/pi^2, col = 'red')

El eje x indica el número de iteraciones, y el eje y la probabilidad estimada con ese número de iteraciones:

Del gráfico podemos ver que, a medida que aumentamos el número de iteraciones, la probabilidad se acerca a $\frac{6}{\pi^2}$ que es justo la probabilidad que calculamos de manera analítica.