Matemáticas - Data Science: Probabilidad

import numpy as np # Librería de manejo de matrices y arreglos de números from numpy.random import binomial # Generador aleatorio de números basado en la distribución binomial from scipy.stats import binom # Permite implementar la función binomial from math import factorial # Permite implementar la función factorial import matplotlib.pyplot as plt # Permite realizar graficos de nuestros datos

# Definición de la deistribución binomial def my_binomial(k, n, p): return factorial(n)/(factorial(k) * factorial(n-k)) * pow(p, k) * pow(1-p, n-k)

print(f'My binomial: {my_binomial(2, 3, 0.5)}') dist = binom(3, 0.5) dist.pmf(2)

# Función de distribución acumulada para los números menores e igual a 2 print(7/8) dist.cdf(2)

# Simulación con 100 lanzamientos de moneda equilibrada # (ejecuta esta celda varias veces para observar la variación en los resultados) p = 0.5 n = 3 binomial(n, p)

# Vamos a hacer un experimento generando una muestra de conjuntos de lanzamientos de a 3 monedas arr = [] for _ in range(100): arr.append(binomial(n, p))

def plot_hist(num_trials): values = [0, 1, 2, 3] arr = [] for _ in range(num_trials): arr.append(binomial(n, p)) simulated_distribution = np.unique(arr, return_counts=True)[1]/len(arr) teoric_distribution = [binom(3, 0.5).pmf(k) for k in values] plt.bar(values, simulated_distribution, color='springgreen') plt.bar(values, teoric_distribution, alpha=0.5, color='darkblue') plt.title(f'{num_trials} experimentos') plt.show() print(f'Distribución Simulada: {np.unique(arr, return_counts=True)[1]/len(arr)}') print(f'Distribución Teórica: {[binom(3, 0.5).pmf(k) for k in values]}')

plot_hist(20)

plot_hist(200)

plot_hist(20000)

import pandas as pd from scipy.stats import norm

def gaussian(x, mu, sigma): return 1/(sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * pow((x-mu)/sigma, 2))

x = np.arange(-4, 4, 0.1) y = gaussian(x, 1.0, 0.5) plt.plot(x, y)

distribution = norm(0, 1) x = np.arange(-4, 4, 0.1) y = [distribution.pdf(value) for value in x] # pdf (probability density function) plt.plot(x, y)

# Scipy nos ayuda a sacar la Distribución de probabilidad acumulada # ya que algunas funciones no son fáciles de integrar distribution = norm(0, 1) x = np.arange(-4, 4, 0.1) y = [distribution.cdf(value) for value in x] # cdf (cumulative distribution function) plt.plot(x, y)

# Para cargar archivos excel pega el siguiente comando: pip install xlrd # al final del archivo init.ipynb de la sección Environment df = pd.read_excel('s057.xls') df

# Trabajamos solo con la columna de tamaño de alas de mosca pero a partir de la posición 3 # ya que las primeras son solo etiquetas como poder en la tabla superior arr = df['Normally Distributed Housefly Wing Lengths'].values[3:] values, dist = np.unique(arr, return_counts=True) plt.bar(values, dist)

# Estimación de distribución mu = arr.mean() sigma = arr.std() x = np.arange(30, 60, 0.1) dist = norm(mu, sigma) y = [dist.pdf(value) for value in x] plt.plot(x, y) # Como la distribución está aplastada por ser probabilidades (van de 0 a 1) # tenemos que normalizar nuestor conteo para mostrar las 2 gráficas juntas values, dist = np.unique(arr, return_counts=True) plt.bar(values, dist/len(arr))

from numpy.random import normal # Generador de números aleatorios basado en la dist. normal

sample = normal(size=10000) # generador plt.hist(sample, bins=100) plt.show()

sample = normal(loc=50, scale=5, size=10000) # mu = 50, sigma = 5 mu = sample.mean() sigma = sample.std() dist = norm(mu, sigma) values = [value for value in range(30, 70)] probabilities = [dist.pdf(value) for value in values] plt.hist(sample, bins=100, density=True) plt.plot(values, probabilities) plt.show()

from numpy import hstack from sklearn.neighbors import KernelDensity # Contruimos una distribución bimodal sample1 = normal(loc=20, scale=5, size=300) sample2 = normal(loc=40, scale=5, size=700) sample = hstack((sample1, sample2)) model = KernelDensity(bandwidth=2, kernel='gaussian') sample = sample.reshape((len(sample), 1)) model.fit(sample)

values = np.asarray([value for value in range(1, 60)]) values = values.reshape((len(values), 1)) probabilities = model.score_samples(values) # probabilidad logarítmica probabilities = np.exp(probabilities) # inversión de probabilidad plt.hist(sample, bins=50, density=True) plt.plot(values, probabilities) plt.show()

from mpl_toolkits.mplot3d import Axes3D # Permite construir un gráfico 3D from matplotlib import cm # Librería para acceder a diferentes paletas de colores

def likelihood(y, yp): return yp*y + (1-yp)*(1-y)

fig = plt.figure() ax = fig.gca(projection='3d') Y = np.arange(0, 1, 0.01) YP = np.arange(0, 1, 0.01) Y, YP = np.meshgrid(Y, YP) # Generamos una malla de valores Z = likelihood(Y, YP) # Ya que la gráfica de una función de 2 variables es una superficie surf = ax.plot_surface(Y, YP, Z, cmap=cm.coolwarm) # Creamos una superficie fig.colorbar(surf, shrink=0.5, aspect=5) # Barra de colores para identificar los máximos plt.show()

from sklearn.datasets import load_iris from sklearn.linear_model import LogisticRegression

atrib_names = ['sepal length', 'sepal width', 'petal length', 'petal width'] X, y = load_iris(return_X_y=True)

X[:10]

y[:100] # Como queremos un problema binario solo escogemos las 2 primeras clases 0 y 1

clasif = LogisticRegression(random_state=10, solver='liblinear').fit(X[:100], y[:100]) # random_state: generador aleatorio que inicializa las variables de nuestro modelo # solver: elegimos el método de optimización de nuestro modelo, en este caso *liblinear*

clasif.coef_

model_coefs = pd.DataFrame(clasif.coef_, columns=atrib_names) model_coefs

clasif.predict(X[:100, :])

clasif.predict_proba(X[:10, :])

from math import factorial # Definición de la distribución binomial def my_binomial(k, n, p): return factorial(n)/(factorial(k)*(factorial(n-k)))*pow(p,k)*pow(1-p, n-k) def my_cumulative_binomial(k, n, p): total = 0 for i in range(k+1): total += my_binomial(i, n, p) return total

my_binomial(3, 12, 0.5)

my_cumulative_binomial(5, 10, 0.5)

my_cumulative_binomial(5, 10, 0.5)

my_binomial(3, 12, 0.3)

my_cumulative_binomial(5, 10, 0.3)

my_cumulative_binomial(5, 10, 0.3)

import numpy as np from scipy.stats import norm import matplotlib.pyplot as plt mu = 4 sigma = 0.1 distribution = norm(mu, sigma)

x = 4 y = distribution.pdf(x) print(y)

x = -10 y = distribution.pdf(x) print(y)

x = 10 y = distribution.pdf(x) print(y)

X = np.arange(3.5, 4, 0.01) y = [distribution.cdf(value) for value in X] plt.plot(X, y) print('P(X <= 4) =',distribution.cdf(4))

X = np.arange(4, 4.5, 0.01) y = [distribution.cdf(value) for value in X] plt.plot(X, y) print('P(X >= 4) =', round(distribution.cdf(5)- distribution.cdf(3.9999), 2))

import numpy as np from numpy.random import binomial import matplotlib.pyplot as plt def generate_binomial_trial(trials, coin_toss): ''' El resultado de esta función es generar un conjunto de experimentos binomiales (trials) y de cada uno obtener las cantidades de éxitos en cada secuencia (ej. lanzar monedas). * trial: es una secuencia de <coin_toss> lanzaientos de moneda * coin_toss: es el número de monedas lanzadas en cada trial ''' arr = [] for _ in range(trials): arr.append(binomial(coin_toss, 0.5)) values, dist = np.unique(arr, return_counts=True) return arr, values, dist

arr, values, dist = generate_binomial_trial(100000, 100) plt.bar(values, dist)

mu = np.mean(arr) sigma = np.std(arr) distribution = norm(mu, sigma) x = np.arange(20, 80, 0.1) y = [distribution.pdf(value) for value in x] plt.plot(x, y, color='orange') # Datos plt.bar(values, dist/len(arr))

import pandas as pd import numpy as np import matplotlib.pyplot as plt # Recuerda que este archivo lo puedes bajar de: # https://seattlecentral.edu/qelp/sets/057/057.html df = pd.read_excel('s057.xls') arr = df['Normally Distributed Housefly Wing Lengths'].values[3:] values, dist = np.unique(arr, return_counts=True) print(values) plt.bar(values, dist)

def optimal_mu(arr=arr): mu = 0 for i in arr: mu += i return mu/len(arr) def optimal_sigma(arr=arr): mu = optimal_mu(arr) sigma = 0 for i in arr: sigma += pow(i - mu, 2) return np.sqrt(sigma/len(arr))

print(optimal_mu(), arr.mean()) print(optimal_sigma(), arr.std())

from scipy.stats import norm values, dist = np.unique(arr, return_counts=True) plt.bar(values, dist/len(arr)) dist = norm(optimal_mu(), optimal_sigma()) x = np.arange(30, 60, 0.1) y = [dist.pdf(value) for value in x] plt.plot(x, y, color='orange')

from sklearn.datasets import make_blobs import matplotlib.pyplot as plt import numpy as np from scipy.stats import norm X, y = make_blobs(n_samples=10000, centers=2, n_features=2, random_state=1) # Esta función ajusta una gaussiana # a un conjunto de data def fit_distribution(data): mu = data.mean() sigma = data.std() dist = norm(mu, sigma) return dist plt.scatter(X[y==1][:, 0], X[y==1][:,1], label='1', color='red') plt.scatter(X[y==0][:, 0], X[y==0][:,1], label='0', color='blue') plt.legend()

# Calculamos priors def prior(c): return len(X[y==c])/len(X) # Tenemos 4 posibles distribuciones a ajustar (verosimilitud) def dist_X0(c): if c == 0: return fit_distribution(X[y==0][:, 0]) elif c == 1: return fit_distribution(X[y==1][:, 0]) def dist_X1(c): if c == 0: return fit_distribution(X[y==0][:, 1]) elif c == 1: return fit_distribution(X[y==1][:, 1]) def likelihood(X, c): return dist_X0(c).pdf(X[0]) * dist_X1(c).pdf(X[1]) # Posterior def probability(c, X): return prior(c) * likelihood(X, c) predictions = [np.argmax([probability(0, vector), probability(1, vector)]) for vector in X]

from sklearn.metrics import confusion_matrix confusion_matrix(y, predictions)

def class_distribution(x, y): print() return np.argmax([probability(0, [x, y]), probability(1, [x, y])]) class_distribution(-6, 0)

class_distribution(-4, 0)

plt.scatter(X[y==1][:, 0], X[y==1][:, 1], label='1', color='red', marker='*') plt.scatter(X[y==0][:, 0], X[y==0][:, 1], label='0', color='blue', marker='*') plt.scatter(-6, 0, color='red', marker='s', s=53) plt.scatter(-4, 0, color='blue', marker='s', s=53) plt.legend()