wehub-resource-sync
749 lines
25 KiB
Python
749 lines
25 KiB
Python
import math
|
|
import random
|
|
|
|
random.seed(42)
|
|
|
|
|
|
def sample_uniform(a, b):
|
|
return a + (b - a) * random.random()
|
|
|
|
|
|
def sample_exponential_inverse_cdf(lam):
|
|
u = random.random()
|
|
return -math.log(u) / lam
|
|
|
|
|
|
def verify_inverse_cdf(lam, n=10000):
|
|
samples = [sample_exponential_inverse_cdf(lam) for _ in range(n)]
|
|
empirical_mean = sum(samples) / len(samples)
|
|
theoretical_mean = 1.0 / lam
|
|
print(f" Exponential(lambda={lam}): empirical mean={empirical_mean:.4f}, "
|
|
f"theoretical={theoretical_mean:.4f}")
|
|
return samples
|
|
|
|
|
|
def normal_pdf(x, mu, sigma):
|
|
coeff = 1.0 / (sigma * math.sqrt(2 * math.pi))
|
|
exponent = -0.5 * ((x - mu) / sigma) ** 2
|
|
return coeff * math.exp(exponent)
|
|
|
|
|
|
def sample_normal_box_muller(mu=0.0, sigma=1.0):
|
|
u1 = random.random()
|
|
u2 = random.random()
|
|
z = math.sqrt(-2 * math.log(u1)) * math.cos(2 * math.pi * u2)
|
|
return mu + sigma * z
|
|
|
|
|
|
def rejection_sample(target_pdf, proposal_sample, proposal_pdf, M):
|
|
attempts = 0
|
|
while True:
|
|
x = proposal_sample()
|
|
u = random.random()
|
|
attempts += 1
|
|
if u < target_pdf(x) / (M * proposal_pdf(x)):
|
|
return x, attempts
|
|
|
|
|
|
def rejection_sample_batch(target_pdf, proposal_sample, proposal_pdf, M, n):
|
|
samples = []
|
|
total_attempts = 0
|
|
for _ in range(n):
|
|
x, attempts = rejection_sample(target_pdf, proposal_sample, proposal_pdf, M)
|
|
samples.append(x)
|
|
total_attempts += attempts
|
|
acceptance_rate = n / total_attempts
|
|
return samples, acceptance_rate
|
|
|
|
|
|
def truncated_normal_demo(mu, sigma, a, b, n=5000):
|
|
norm_const = sum(
|
|
normal_pdf(a + (b - a) * i / 1000, mu, sigma) * (b - a) / 1000
|
|
for i in range(1001)
|
|
)
|
|
M_val = max(
|
|
normal_pdf(x, mu, sigma) / (1.0 / (b - a))
|
|
for x in [a + (b - a) * i / 200 for i in range(201)]
|
|
)
|
|
|
|
def target(x):
|
|
if a <= x <= b:
|
|
return normal_pdf(x, mu, sigma)
|
|
return 0.0
|
|
|
|
def proposal():
|
|
return sample_uniform(a, b)
|
|
|
|
def proposal_pdf(x):
|
|
if a <= x <= b:
|
|
return 1.0 / (b - a)
|
|
return 0.0
|
|
|
|
samples, acc_rate = rejection_sample_batch(target, proposal, proposal_pdf, M_val, n)
|
|
return samples, acc_rate
|
|
|
|
|
|
def importance_sampling_estimate(f, target_pdf, proposal_pdf, proposal_sample, n):
|
|
weighted_sum = 0.0
|
|
weight_sum = 0.0
|
|
for _ in range(n):
|
|
x = proposal_sample()
|
|
w = target_pdf(x) / proposal_pdf(x)
|
|
weighted_sum += f(x) * w
|
|
weight_sum += w
|
|
unnormalized = weighted_sum / n
|
|
self_normalized = weighted_sum / weight_sum
|
|
return unnormalized, self_normalized
|
|
|
|
|
|
def importance_sampling_demo():
|
|
mu, sigma = 2.0, 1.0
|
|
a, b = -3.0, 7.0
|
|
|
|
def f(x):
|
|
return x ** 2
|
|
|
|
def target(x):
|
|
return normal_pdf(x, mu, sigma)
|
|
|
|
def proposal():
|
|
return sample_uniform(a, b)
|
|
|
|
def proposal_pdf(x):
|
|
if a <= x <= b:
|
|
return 1.0 / (b - a)
|
|
return 0.0
|
|
|
|
est_unnorm, est_selfnorm = importance_sampling_estimate(
|
|
f, target, proposal_pdf, proposal, n=50000
|
|
)
|
|
theoretical = mu ** 2 + sigma ** 2
|
|
print(f" E[X^2] under N({mu},{sigma}):")
|
|
print(f" Unnormalized IS: {est_unnorm:.4f}")
|
|
print(f" Self-normalized IS: {est_selfnorm:.4f}")
|
|
print(f" Theoretical: {theoretical:.4f}")
|
|
|
|
|
|
def monte_carlo_pi(n):
|
|
inside = 0
|
|
for _ in range(n):
|
|
x = random.uniform(-1, 1)
|
|
y = random.uniform(-1, 1)
|
|
if x * x + y * y <= 1:
|
|
inside += 1
|
|
return 4 * inside / n
|
|
|
|
|
|
def monte_carlo_integral(f, a, b, n):
|
|
total = 0.0
|
|
for _ in range(n):
|
|
x = sample_uniform(a, b)
|
|
total += f(x)
|
|
return (b - a) * total / n
|
|
|
|
|
|
def metropolis_hastings(target_log_pdf, x0, n_samples, burn_in, proposal_std=1.0):
|
|
samples = []
|
|
x = x0
|
|
accepted = 0
|
|
total = n_samples + burn_in
|
|
|
|
for i in range(total):
|
|
x_new = x + random.gauss(0, proposal_std)
|
|
log_alpha = target_log_pdf(x_new) - target_log_pdf(x)
|
|
|
|
if math.log(random.random() + 1e-300) < log_alpha:
|
|
x = x_new
|
|
if i >= burn_in:
|
|
accepted += 1
|
|
|
|
if i >= burn_in:
|
|
samples.append(x)
|
|
|
|
acceptance_rate = accepted / n_samples
|
|
return samples, acceptance_rate
|
|
|
|
|
|
def metropolis_hastings_2d(target_log_pdf, x0, y0, n_samples, burn_in, proposal_std=0.5):
|
|
samples = []
|
|
x, y = x0, y0
|
|
accepted = 0
|
|
total = n_samples + burn_in
|
|
|
|
for i in range(total):
|
|
x_new = x + random.gauss(0, proposal_std)
|
|
y_new = y + random.gauss(0, proposal_std)
|
|
log_alpha = target_log_pdf(x_new, y_new) - target_log_pdf(x, y)
|
|
|
|
if math.log(random.random() + 1e-300) < log_alpha:
|
|
x, y = x_new, y_new
|
|
if i >= burn_in:
|
|
accepted += 1
|
|
|
|
if i >= burn_in:
|
|
samples.append((x, y))
|
|
|
|
acceptance_rate = accepted / n_samples
|
|
return samples, acceptance_rate
|
|
|
|
|
|
def bimodal_log_pdf(x):
|
|
p1 = 0.4 * normal_pdf(x, -3, 1)
|
|
p2 = 0.6 * normal_pdf(x, 3, 1)
|
|
return math.log(p1 + p2 + 1e-300)
|
|
|
|
|
|
def gibbs_sampling_2d(rho, n_samples, burn_in):
|
|
x, y = 0.0, 0.0
|
|
samples = []
|
|
|
|
for i in range(n_samples + burn_in):
|
|
x = random.gauss(rho * y, math.sqrt(1 - rho ** 2))
|
|
y = random.gauss(rho * x, math.sqrt(1 - rho ** 2))
|
|
if i >= burn_in:
|
|
samples.append((x, y))
|
|
|
|
return samples
|
|
|
|
|
|
def softmax(logits):
|
|
max_l = max(logits)
|
|
exps = [math.exp(z - max_l) for z in logits]
|
|
total = sum(exps)
|
|
return [e / total for e in exps]
|
|
|
|
|
|
def sample_from_probs(probs):
|
|
r = random.random()
|
|
cumsum = 0.0
|
|
for i, p in enumerate(probs):
|
|
cumsum += p
|
|
if r <= cumsum:
|
|
return i
|
|
return len(probs) - 1
|
|
|
|
|
|
def temperature_sample(logits, temperature):
|
|
if temperature <= 0:
|
|
return logits.index(max(logits))
|
|
scaled = [z / temperature for z in logits]
|
|
probs = softmax(scaled)
|
|
return sample_from_probs(probs)
|
|
|
|
|
|
def temperature_distribution(logits, temperature):
|
|
if temperature <= 0:
|
|
result = [0.0] * len(logits)
|
|
result[logits.index(max(logits))] = 1.0
|
|
return result
|
|
scaled = [z / temperature for z in logits]
|
|
return softmax(scaled)
|
|
|
|
|
|
def top_k_sample(logits, k):
|
|
indexed = sorted(enumerate(logits), key=lambda x: -x[1])
|
|
top = indexed[:k]
|
|
top_logits = [l for _, l in top]
|
|
probs = softmax(top_logits)
|
|
idx = sample_from_probs(probs)
|
|
return top[idx][0]
|
|
|
|
|
|
def top_k_distribution(logits, k):
|
|
probs = softmax(logits)
|
|
indexed = sorted(enumerate(probs), key=lambda x: -x[1])
|
|
result = [0.0] * len(logits)
|
|
top_indices = [idx for idx, _ in indexed[:k]]
|
|
top_probs = [probs[idx] for idx in top_indices]
|
|
total = sum(top_probs)
|
|
for idx in top_indices:
|
|
result[idx] = probs[idx] / total
|
|
return result
|
|
|
|
|
|
def top_p_sample(logits, p):
|
|
probs = softmax(logits)
|
|
indexed = sorted(enumerate(probs), key=lambda x: -x[1])
|
|
cumsum = 0.0
|
|
selected = []
|
|
for token_idx, prob in indexed:
|
|
cumsum += prob
|
|
selected.append((token_idx, prob))
|
|
if cumsum >= p:
|
|
break
|
|
sel_probs = [pr for _, pr in selected]
|
|
total = sum(sel_probs)
|
|
sel_probs = [pr / total for pr in sel_probs]
|
|
idx = sample_from_probs(sel_probs)
|
|
return selected[idx][0]
|
|
|
|
|
|
def top_p_distribution(logits, p):
|
|
probs = softmax(logits)
|
|
indexed = sorted(enumerate(probs), key=lambda x: -x[1])
|
|
cumsum = 0.0
|
|
selected_indices = []
|
|
for token_idx, prob in indexed:
|
|
cumsum += prob
|
|
selected_indices.append(token_idx)
|
|
if cumsum >= p:
|
|
break
|
|
result = [0.0] * len(logits)
|
|
total = sum(probs[i] for i in selected_indices)
|
|
for i in selected_indices:
|
|
result[i] = probs[i] / total
|
|
return result
|
|
|
|
|
|
def reparam_sample(mu, sigma):
|
|
epsilon = random.gauss(0, 1)
|
|
z = mu + sigma * epsilon
|
|
return z, epsilon
|
|
|
|
|
|
def reparam_gradient(epsilon):
|
|
dz_dmu = 1.0
|
|
dz_dsigma = epsilon
|
|
return dz_dmu, dz_dsigma
|
|
|
|
|
|
def vae_forward_demo(mu, log_var):
|
|
sigma = math.exp(0.5 * log_var)
|
|
z, epsilon = reparam_sample(mu, sigma)
|
|
|
|
dz_dmu = 1.0
|
|
dz_dsigma = epsilon
|
|
dsigma_dlogvar = 0.5 * sigma
|
|
|
|
dz_dmu_total = dz_dmu
|
|
dz_dlogvar = dz_dsigma * dsigma_dlogvar
|
|
|
|
return z, epsilon, dz_dmu_total, dz_dlogvar
|
|
|
|
|
|
def gumbel_sample():
|
|
u = random.random()
|
|
while u == 0:
|
|
u = random.random()
|
|
return -math.log(-math.log(u))
|
|
|
|
|
|
def gumbel_max_sample(log_probs):
|
|
gumbels = [lp + gumbel_sample() for lp in log_probs]
|
|
return gumbels.index(max(gumbels))
|
|
|
|
|
|
def gumbel_softmax_sample(log_probs, temperature):
|
|
gumbels = [lp + gumbel_sample() for lp in log_probs]
|
|
scaled = [g / temperature for g in gumbels]
|
|
return softmax(scaled)
|
|
|
|
|
|
def gumbel_softmax_straight_through(log_probs, temperature):
|
|
soft = gumbel_softmax_sample(log_probs, temperature)
|
|
hard_idx = soft.index(max(soft))
|
|
hard = [0.0] * len(log_probs)
|
|
hard[hard_idx] = 1.0
|
|
return hard, soft
|
|
|
|
|
|
def stratified_sample_1d(n):
|
|
samples = []
|
|
for i in range(n):
|
|
u = random.random()
|
|
samples.append((i + u) / n)
|
|
return samples
|
|
|
|
|
|
def compare_sampling_variance(f, n, n_trials=200):
|
|
standard_estimates = []
|
|
stratified_estimates = []
|
|
|
|
for _ in range(n_trials):
|
|
standard_samples = [random.random() for _ in range(n)]
|
|
standard_est = sum(f(x) for x in standard_samples) / n
|
|
standard_estimates.append(standard_est)
|
|
|
|
strat_samples = stratified_sample_1d(n)
|
|
strat_est = sum(f(x) for x in strat_samples) / n
|
|
stratified_estimates.append(strat_est)
|
|
|
|
std_mean = sum(standard_estimates) / len(standard_estimates)
|
|
std_var = sum((e - std_mean) ** 2 for e in standard_estimates) / len(standard_estimates)
|
|
|
|
strat_mean = sum(stratified_estimates) / len(stratified_estimates)
|
|
strat_var = sum((e - strat_mean) ** 2 for e in stratified_estimates) / len(stratified_estimates)
|
|
|
|
return std_var, strat_var
|
|
|
|
|
|
def text_generation_demo(vocab, logits, length, method, **kwargs):
|
|
tokens = []
|
|
for _ in range(length):
|
|
if method == "greedy":
|
|
idx = logits.index(max(logits))
|
|
elif method == "temperature":
|
|
idx = temperature_sample(logits, kwargs.get("temperature", 1.0))
|
|
elif method == "top_k":
|
|
idx = top_k_sample(logits, kwargs.get("k", 5))
|
|
elif method == "top_p":
|
|
idx = top_p_sample(logits, kwargs.get("p", 0.9))
|
|
else:
|
|
idx = sample_from_probs(softmax(logits))
|
|
tokens.append(vocab[idx])
|
|
return " ".join(tokens)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
print("=" * 65)
|
|
print("SAMPLING METHODS")
|
|
print("=" * 65)
|
|
|
|
print("\n--- 1. Inverse CDF Sampling (Exponential) ---")
|
|
for lam in [0.5, 1.0, 2.0]:
|
|
verify_inverse_cdf(lam)
|
|
|
|
print("\n--- 2. Rejection Sampling (Truncated Normal) ---")
|
|
trunc_samples, acc = truncated_normal_demo(0, 1, -1, 2, n=5000)
|
|
trunc_mean = sum(trunc_samples) / len(trunc_samples)
|
|
print(f" Truncated N(0,1) on [-1, 2]: mean={trunc_mean:.4f}, "
|
|
f"acceptance rate={acc:.4f}")
|
|
|
|
print("\n--- 3. Importance Sampling ---")
|
|
importance_sampling_demo()
|
|
|
|
print("\n--- 4. Monte Carlo Estimation ---")
|
|
print(" Estimating pi:")
|
|
for n in [1000, 10000, 100000]:
|
|
pi_est = monte_carlo_pi(n)
|
|
error = abs(pi_est - math.pi)
|
|
print(f" N={n:>7d}: pi ~ {pi_est:.6f}, error={error:.6f}")
|
|
|
|
print(" Estimating integral of sin(x) from 0 to pi (true = 2.0):")
|
|
for n in [1000, 10000, 100000]:
|
|
est = monte_carlo_integral(math.sin, 0, math.pi, n)
|
|
error = abs(est - 2.0)
|
|
print(f" N={n:>7d}: estimate={est:.6f}, error={error:.6f}")
|
|
|
|
print("\n--- 5. Metropolis-Hastings MCMC ---")
|
|
print(" Sampling from bimodal distribution (mixture of Gaussians):")
|
|
for std in [0.5, 1.0, 3.0]:
|
|
samples_mh, acc_rate = metropolis_hastings(
|
|
bimodal_log_pdf, x0=0.0, n_samples=10000, burn_in=2000,
|
|
proposal_std=std
|
|
)
|
|
mh_mean = sum(samples_mh) / len(samples_mh)
|
|
mh_std = (sum((x - mh_mean) ** 2 for x in samples_mh) / len(samples_mh)) ** 0.5
|
|
print(f" proposal_std={std}: mean={mh_mean:.4f}, std={mh_std:.4f}, "
|
|
f"acceptance={acc_rate:.4f}")
|
|
|
|
print("\n Sampling from 2D Gaussian:")
|
|
def gaussian_2d_log_pdf(x, y):
|
|
return -0.5 * (x ** 2 + y ** 2)
|
|
|
|
samples_2d, acc_2d = metropolis_hastings_2d(
|
|
gaussian_2d_log_pdf, x0=5.0, y0=5.0,
|
|
n_samples=10000, burn_in=2000, proposal_std=1.0
|
|
)
|
|
xs = [s[0] for s in samples_2d]
|
|
ys = [s[1] for s in samples_2d]
|
|
print(f" mean_x={sum(xs)/len(xs):.4f}, mean_y={sum(ys)/len(ys):.4f}, "
|
|
f"acceptance={acc_2d:.4f}")
|
|
|
|
print("\n--- 6. Gibbs Sampling (Correlated 2D Gaussian) ---")
|
|
for rho in [0.0, 0.5, 0.9]:
|
|
gibbs_samples = gibbs_sampling_2d(rho, n_samples=10000, burn_in=1000)
|
|
gx = [s[0] for s in gibbs_samples]
|
|
gy = [s[1] for s in gibbs_samples]
|
|
empirical_corr = (
|
|
sum(a * b for a, b in gibbs_samples) / len(gibbs_samples)
|
|
- (sum(gx) / len(gx)) * (sum(gy) / len(gy))
|
|
)
|
|
print(f" rho={rho}: empirical correlation={empirical_corr:.4f}")
|
|
|
|
print("\n--- 7. Temperature Sampling ---")
|
|
token_logits = [3.0, 2.0, 1.5, 0.5, -1.0, -2.0]
|
|
vocab = ["the", "a", "this", "one", "that", "some"]
|
|
print(" Logits:", token_logits)
|
|
print(" Vocab:", vocab)
|
|
for temp in [0.1, 0.5, 1.0, 1.5, 3.0]:
|
|
dist = temperature_distribution(token_logits, temp)
|
|
formatted = [f"{p:.4f}" for p in dist]
|
|
print(f" T={temp:.1f}: [{', '.join(formatted)}]")
|
|
|
|
print("\n Generation samples at different temperatures:")
|
|
for temp in [0.1, 0.7, 1.0, 2.0]:
|
|
tokens_out = []
|
|
for _ in range(10):
|
|
idx = temperature_sample(token_logits, temp)
|
|
tokens_out.append(vocab[idx])
|
|
print(f" T={temp}: {' '.join(tokens_out)}")
|
|
|
|
print("\n--- 8. Top-k Sampling ---")
|
|
print(" Logits:", token_logits)
|
|
for k in [1, 2, 3, 6]:
|
|
dist = top_k_distribution(token_logits, k)
|
|
formatted = [f"{p:.4f}" for p in dist]
|
|
print(f" k={k}: [{', '.join(formatted)}]")
|
|
|
|
print("\n--- 9. Top-p (Nucleus) Sampling ---")
|
|
print(" Logits:", token_logits)
|
|
for p in [0.5, 0.8, 0.9, 0.95, 1.0]:
|
|
dist = top_p_distribution(token_logits, p)
|
|
nonzero = sum(1 for d in dist if d > 0)
|
|
formatted = [f"{d:.4f}" for d in dist]
|
|
print(f" p={p:.2f}: [{', '.join(formatted)}] ({nonzero} tokens)")
|
|
|
|
print("\n--- 10. Reparameterization Trick ---")
|
|
mu_val, log_var_val = 2.0, 0.5
|
|
|
|
print(f" VAE encoder output: mu={mu_val}, log_var={log_var_val}")
|
|
z_samples = []
|
|
for trial in range(5):
|
|
z, eps, dz_dmu, dz_dlogvar = vae_forward_demo(mu_val, log_var_val)
|
|
z_samples.append(z)
|
|
print(f" Trial {trial+1}: z={z:.4f}, eps={eps:.4f}, "
|
|
f"dz/dmu={dz_dmu:.4f}, dz/dlog_var={dz_dlogvar:.4f}")
|
|
|
|
print(f" Mean of z samples: {sum(z_samples)/len(z_samples):.4f} "
|
|
f"(expected ~{mu_val})")
|
|
print(" Gradients exist because z = mu + sigma * epsilon is differentiable.")
|
|
|
|
print("\n Verifying reparameterization matches direct sampling:")
|
|
sigma_val = math.exp(0.5 * log_var_val)
|
|
direct_samples = [sample_normal_box_muller(mu_val, sigma_val) for _ in range(10000)]
|
|
reparam_samples = [reparam_sample(mu_val, sigma_val)[0] for _ in range(10000)]
|
|
d_mean = sum(direct_samples) / len(direct_samples)
|
|
r_mean = sum(reparam_samples) / len(reparam_samples)
|
|
d_std = (sum((x - d_mean)**2 for x in direct_samples) / len(direct_samples)) ** 0.5
|
|
r_std = (sum((x - r_mean)**2 for x in reparam_samples) / len(reparam_samples)) ** 0.5
|
|
print(f" Direct: mean={d_mean:.4f}, std={d_std:.4f}")
|
|
print(f" Reparam: mean={r_mean:.4f}, std={r_std:.4f}")
|
|
|
|
print("\n--- 11. Gumbel-Softmax ---")
|
|
probs = [0.5, 0.3, 0.15, 0.05]
|
|
log_probs = [math.log(p) for p in probs]
|
|
labels = ["cat", "dog", "bird", "fish"]
|
|
print(f" True probs: {probs}")
|
|
|
|
print("\n Gumbel-Max (exact categorical) verification:")
|
|
counts = [0] * len(probs)
|
|
n_gumbel = 10000
|
|
for _ in range(n_gumbel):
|
|
idx = gumbel_max_sample(log_probs)
|
|
counts[idx] += 1
|
|
empirical = [c / n_gumbel for c in counts]
|
|
print(f" Empirical: [{', '.join(f'{p:.4f}' for p in empirical)}]")
|
|
print(f" True: [{', '.join(f'{p:.4f}' for p in probs)}]")
|
|
|
|
print("\n Gumbel-Softmax at different temperatures:")
|
|
for tau in [0.1, 0.5, 1.0, 5.0]:
|
|
soft = gumbel_softmax_sample(log_probs, tau)
|
|
formatted = [f"{s:.4f}" for s in soft]
|
|
max_idx = soft.index(max(soft))
|
|
print(f" tau={tau:.1f}: [{', '.join(formatted)}] -> {labels[max_idx]}")
|
|
|
|
print("\n Straight-through estimator:")
|
|
hard, soft = gumbel_softmax_straight_through(log_probs, temperature=0.5)
|
|
print(f" Hard (forward): {hard}")
|
|
print(f" Soft (backward): [{', '.join(f'{s:.4f}' for s in soft)}]")
|
|
|
|
print("\n--- 12. Stratified Sampling ---")
|
|
def test_fn(x):
|
|
return math.sin(math.pi * x)
|
|
|
|
print(" Comparing standard vs stratified Monte Carlo:")
|
|
print(f" Function: sin(pi*x) on [0,1], true integral = 2/pi = {2/math.pi:.6f}")
|
|
for n in [10, 50, 100]:
|
|
std_var, strat_var = compare_sampling_variance(test_fn, n)
|
|
ratio = std_var / strat_var if strat_var > 0 else float('inf')
|
|
print(f" N={n:3d}: standard_var={std_var:.8f}, "
|
|
f"stratified_var={strat_var:.8f}, ratio={ratio:.2f}x")
|
|
|
|
print("\n--- 13. Text Generation Demo ---")
|
|
gen_vocab = ["the", "cat", "sat", "on", "mat", "a", "dog", "ran", "big", "red"]
|
|
gen_logits = [3.0, 2.5, 2.0, 1.8, 1.5, 1.0, 0.5, 0.0, -0.5, -1.0]
|
|
|
|
print(f" Vocab: {gen_vocab}")
|
|
print(f" Logits: {gen_logits}")
|
|
print()
|
|
|
|
methods = [
|
|
("greedy", {}),
|
|
("temperature", {"temperature": 0.5}),
|
|
("temperature", {"temperature": 1.0}),
|
|
("temperature", {"temperature": 2.0}),
|
|
("top_k", {"k": 3}),
|
|
("top_p", {"p": 0.8}),
|
|
]
|
|
|
|
for method_name, params in methods:
|
|
label = method_name
|
|
if params:
|
|
label += f"({', '.join(f'{k}={v}' for k, v in params.items())})"
|
|
sequences = []
|
|
for run in range(3):
|
|
seq = text_generation_demo(gen_vocab, gen_logits, length=8,
|
|
method=method_name, **params)
|
|
sequences.append(seq)
|
|
print(f" {label}:")
|
|
for i, seq in enumerate(sequences):
|
|
print(f" Run {i+1}: {seq}")
|
|
unique = len(set(sequences))
|
|
print(f" Unique sequences: {unique}/3")
|
|
print()
|
|
|
|
print("\n--- 14. Visualizations ---")
|
|
try:
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
|
|
fig, axes = plt.subplots(3, 3, figsize=(18, 16))
|
|
|
|
ax = axes[0][0]
|
|
ax.set_title("Inverse CDF: Exponential Samples")
|
|
exp_samples = [sample_exponential_inverse_cdf(1.0) for _ in range(10000)]
|
|
ax.hist(exp_samples, bins=60, density=True, alpha=0.7, color="#4a90d9",
|
|
label="Samples")
|
|
xs_exp = [i * 0.05 for i in range(160)]
|
|
ys_exp = [math.exp(-x) for x in xs_exp]
|
|
ax.plot(xs_exp, ys_exp, "r-", linewidth=2, label="True PDF")
|
|
ax.set_xlabel("x")
|
|
ax.set_ylabel("Density")
|
|
ax.legend()
|
|
|
|
ax = axes[0][1]
|
|
ax.set_title("Rejection Sampling: Truncated Normal")
|
|
rej_samples, _ = truncated_normal_demo(0, 1, -1, 2, n=5000)
|
|
ax.hist(rej_samples, bins=50, density=True, alpha=0.7, color="#4a90d9",
|
|
label="Samples")
|
|
xs_tn = [-1 + 3 * i / 200 for i in range(201)]
|
|
ys_tn = [normal_pdf(x, 0, 1) for x in xs_tn]
|
|
area = sum(ys_tn) * 3 / 200
|
|
ys_tn_norm = [y / area for y in ys_tn]
|
|
ax.plot(xs_tn, ys_tn_norm, "r-", linewidth=2, label="True PDF (normalized)")
|
|
ax.set_xlabel("x")
|
|
ax.set_ylabel("Density")
|
|
ax.legend()
|
|
|
|
ax = axes[0][2]
|
|
ax.set_title("Monte Carlo: Estimating Pi")
|
|
n_mc_vis = 5000
|
|
mc_x = [random.uniform(-1, 1) for _ in range(n_mc_vis)]
|
|
mc_y = [random.uniform(-1, 1) for _ in range(n_mc_vis)]
|
|
inside_x = [mc_x[i] for i in range(n_mc_vis) if mc_x[i]**2 + mc_y[i]**2 <= 1]
|
|
inside_y = [mc_y[i] for i in range(n_mc_vis) if mc_x[i]**2 + mc_y[i]**2 <= 1]
|
|
outside_x = [mc_x[i] for i in range(n_mc_vis) if mc_x[i]**2 + mc_y[i]**2 > 1]
|
|
outside_y = [mc_y[i] for i in range(n_mc_vis) if mc_x[i]**2 + mc_y[i]**2 > 1]
|
|
ax.scatter(inside_x, inside_y, s=1, c="#4a90d9", alpha=0.5)
|
|
ax.scatter(outside_x, outside_y, s=1, c="#d94a4a", alpha=0.5)
|
|
theta = [2 * math.pi * i / 200 for i in range(201)]
|
|
circle_x = [math.cos(t) for t in theta]
|
|
circle_y = [math.sin(t) for t in theta]
|
|
ax.plot(circle_x, circle_y, "k-", linewidth=1.5)
|
|
ax.set_aspect("equal")
|
|
pi_est = 4 * len(inside_x) / n_mc_vis
|
|
ax.set_xlabel(f"pi ~ {pi_est:.4f}")
|
|
|
|
ax = axes[1][0]
|
|
ax.set_title("MCMC: Bimodal Distribution")
|
|
mcmc_samples, _ = metropolis_hastings(
|
|
bimodal_log_pdf, x0=0.0, n_samples=20000, burn_in=5000, proposal_std=2.0
|
|
)
|
|
ax.hist(mcmc_samples, bins=80, density=True, alpha=0.7, color="#4a90d9",
|
|
label="MCMC samples")
|
|
xs_bm = [-8 + 16 * i / 400 for i in range(401)]
|
|
ys_bm = [math.exp(bimodal_log_pdf(x)) for x in xs_bm]
|
|
area_bm = sum(ys_bm) * 16 / 400
|
|
ys_bm_norm = [y / area_bm for y in ys_bm]
|
|
ax.plot(xs_bm, ys_bm_norm, "r-", linewidth=2, label="True density")
|
|
ax.set_xlabel("x")
|
|
ax.set_ylabel("Density")
|
|
ax.legend()
|
|
|
|
ax = axes[1][1]
|
|
ax.set_title("Gibbs Sampling: 2D Gaussian (rho=0.8)")
|
|
gibbs_vis = gibbs_sampling_2d(0.8, n_samples=3000, burn_in=500)
|
|
gvx = [s[0] for s in gibbs_vis]
|
|
gvy = [s[1] for s in gibbs_vis]
|
|
ax.scatter(gvx, gvy, s=2, alpha=0.3, c="#4a90d9")
|
|
ax.plot(gvx[:100], gvy[:100], "r-", alpha=0.3, linewidth=0.5)
|
|
ax.set_xlabel("x")
|
|
ax.set_ylabel("y")
|
|
ax.set_aspect("equal")
|
|
|
|
ax = axes[1][2]
|
|
ax.set_title("Temperature Scaling")
|
|
temps = [0.1, 0.5, 1.0, 2.0, 5.0]
|
|
bar_width = 0.15
|
|
positions = list(range(len(token_logits)))
|
|
for t_idx, temp in enumerate(temps):
|
|
dist = temperature_distribution(token_logits, temp)
|
|
offset = (t_idx - 2) * bar_width
|
|
bars = [pos + offset for pos in positions]
|
|
ax.bar(bars, dist, bar_width, label=f"T={temp}", alpha=0.8)
|
|
ax.set_xticks(positions)
|
|
ax.set_xticklabels(vocab, rotation=45)
|
|
ax.set_ylabel("Probability")
|
|
ax.legend(fontsize=8)
|
|
|
|
ax = axes[2][0]
|
|
ax.set_title("Top-k vs Top-p Distributions")
|
|
k_dist = top_k_distribution(token_logits, k=3)
|
|
p_dist = top_p_distribution(token_logits, p=0.9)
|
|
full_dist = softmax(token_logits)
|
|
x_pos = list(range(len(token_logits)))
|
|
w = 0.25
|
|
ax.bar([x - w for x in x_pos], full_dist, w, label="Full", alpha=0.8, color="#aaaaaa")
|
|
ax.bar(x_pos, k_dist, w, label="Top-3", alpha=0.8, color="#4a90d9")
|
|
ax.bar([x + w for x in x_pos], p_dist, w, label="Top-p=0.9", alpha=0.8, color="#d94a4a")
|
|
ax.set_xticks(x_pos)
|
|
ax.set_xticklabels(vocab, rotation=45)
|
|
ax.set_ylabel("Probability")
|
|
ax.legend(fontsize=8)
|
|
|
|
ax = axes[2][1]
|
|
ax.set_title("Gumbel-Softmax: Temperature Effect")
|
|
taus = [0.1, 0.5, 1.0, 5.0]
|
|
g_log_probs = [math.log(p) for p in [0.5, 0.3, 0.15, 0.05]]
|
|
n_trials_vis = 500
|
|
for tau in taus:
|
|
max_vals = []
|
|
for _ in range(n_trials_vis):
|
|
soft = gumbel_softmax_sample(g_log_probs, tau)
|
|
max_vals.append(max(soft))
|
|
ax.hist(max_vals, bins=30, alpha=0.5, label=f"tau={tau}", density=True)
|
|
ax.set_xlabel("Max component value")
|
|
ax.set_ylabel("Density")
|
|
ax.legend(fontsize=8)
|
|
|
|
ax = axes[2][2]
|
|
ax.set_title("Stratified vs Standard Sampling")
|
|
n_strat_vis = 20
|
|
standard_pts = sorted([random.random() for _ in range(n_strat_vis)])
|
|
stratified_pts = sorted(stratified_sample_1d(n_strat_vis))
|
|
ax.scatter(standard_pts, [1] * n_strat_vis, s=30, c="#d94a4a", label="Standard",
|
|
zorder=3)
|
|
ax.scatter(stratified_pts, [0] * n_strat_vis, s=30, c="#4a90d9", label="Stratified",
|
|
zorder=3)
|
|
for i in range(n_strat_vis + 1):
|
|
ax.axvline(i / n_strat_vis, color="#cccccc", linewidth=0.5, linestyle="--")
|
|
ax.set_yticks([0, 1])
|
|
ax.set_yticklabels(["Stratified", "Standard"])
|
|
ax.set_xlabel("Sample value")
|
|
ax.legend()
|
|
ax.set_ylim(-0.5, 1.5)
|
|
|
|
plt.tight_layout()
|
|
plt.savefig("sampling_methods.png", dpi=150)
|
|
print(" Saved: sampling_methods.png")
|
|
plt.close()
|
|
|
|
except ImportError:
|
|
print(" matplotlib not available, skipping visualization.")
|
|
|
|
print("\n" + "=" * 65)
|
|
print("All sampling methods complete.")
|
|
print("=" * 65)
|