Skip to content

Instantly share code, notes, and snippets.

@brando90
Last active December 23, 2024 03:22
Show Gist options
  • Save brando90/72a1e57de2774c6a9942e67d246f54f7 to your computer and use it in GitHub Desktop.
Save brando90/72a1e57de2774c6a9942e67d246f54f7 to your computer and use it in GitHub Desktop.
method1_replace_0s_1s_with_rounding_err.py
import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# Generate synthetic data with latent Beta distribution
np.random.seed(42)
alpha_true, beta_true = 2, 5 # True Beta distribution parameters
n_samples = 1000
resolution = 1e-4
# Latent values from the Beta distribution
latent_values = beta.rvs(alpha_true, beta_true, size=n_samples)
# Observed values after rounding to the nearest resolution
observed_values = np.round(latent_values / resolution) * resolution
# Replace exact zeros with a small offset
offset = 1e-5
adjusted_values = np.where(observed_values == 0, offset, observed_values)
# Plot latent vs observed values
plt.hist(latent_values, bins=50, alpha=0.5, label="Latent Beta values", color="blue")
plt.hist(adjusted_values, bins=50, alpha=0.5, label="Adjusted Observed values", color="orange")
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.legend()
plt.title("Latent vs Adjusted Observed Values")
plt.show()
# Fit a Beta distribution to adjusted values using maximum likelihood estimation
def negative_log_likelihood(params, data):
alpha, beta_param = params
if alpha <= 0 or beta_param <= 0:
return np.inf # Ensure positivity of parameters
return -np.sum(beta.logpdf(data, alpha, beta_param))
# Initial guesses for alpha and beta
initial_guess = [1.0, 1.0]
# Minimize the negative log-likelihood
result = minimize(negative_log_likelihood, initial_guess, args=(adjusted_values,))
fitted_alpha, fitted_beta = result.x
# Print results
fitted_alpha, fitted_beta
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment