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)
My binomial: 0.375
# Función de distribución acumulada para los números menores e igual a 2
print(7/8)
dist.cdf(2)
0.875
# 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)
Distribución Simulada: [0.25 0.4 0.3 0.05]
Distribución Teórica: [0.125, 0.3750000000000001, 0.3750000000000001, 0.125]
plot_hist(200)
Distribución Simulada: [0.14 0.425 0.29 0.145]
Distribución Teórica: [0.125, 0.3750000000000001, 0.3750000000000001, 0.125]
plot_hist(20000)
Distribución Simulada: [0.12485 0.3692 0.3824 0.12355]
Distribución Teórica: [0.125, 0.3750000000000001, 0.3750000000000001, 0.125]
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()
/tmp/ipykernel_84/2716258221.py:2: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
ax = fig.gca(projection='3d')
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)
3.989422804014327
x = -10
y = distribution.pdf(x)
print(y)
0.0
x = 10
y = distribution.pdf(x)
print(y)
0.0
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))
P(X <= 4) = 0.5
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))
P(X >= 4) = 0.5
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)
[36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55]
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())
45.5 45.5
3.9 3.9
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()