Parameter Uncertainty¶

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Height
axes[0].plot(t_span, height, "--", label="True", linewidth=2, alpha=0.7)
axes[0].plot(t_span, height_noisy, "o", label="Observed", alpha=0.6)
axes[0].set_xlabel("Time (s)")
axes[0].set_ylabel("Height (m)")
axes[0].set_title("Ball Height")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Velocity
axes[1].plot(t_span, velocity, "--", label="True", linewidth=2, alpha=0.7)
axes[1].plot(t_span, velocity_noisy, "o", label="Observed", alpha=0.6)
axes[1].set_xlabel("Time (s)")
axes[1].set_ylabel("Velocity (m/s)")
axes[1].set_title("Ball Velocity")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# DiffSL model for falling ball
dsl_model = """
in_i {g = 1, h = 1 }
u_i {x = h, v = 0}
F_i {v, -g}
stop {x}
"""

print("DiffSL Model:")
print(dsl_model)
print("\nExplanation:")
print("  x = height, v = velocity")
print("  dx/dt = v (velocity determines height change)")
print("  dv/dt = -g (gravity accelerates downward)")
print("  stop {x} (terminate when height reaches zero)")

# Plot traces
fig, axes = parameter_traces(
    mcmc_result.samples,
    param_names=["g (m/s²)", "h (m)"],
    true_values=[g_true, h_true],
    show=False,
)

plt.show()

print("\nWhat to look for in traces:")
print("  - Good mixing (wiggly, no trends)")
print("  - Stationary (mean stays constant)")
print("  - No long excursions")
print("  - Rapid exploration of parameter space")

# Plot distributions
fig, axes = parameter_distributions(
    mcmc_result.samples,
    param_names=["g (m/s²)", "h (m)"],
    true_values=[g_true, h_true],
    show=False,
)

plt.show()

# Burn-in: discard first 20% of samples
burn_in = int(0.2 * len(mcmc_result.samples))
samples_burned = mcmc_result.samples[burn_in:]

# Calculate statistics
g_samples = samples_burned[:, 0]
h_samples = samples_burned[:, 1]

# Means and standard deviations
g_mean = np.mean(g_samples)
g_std = np.std(g_samples)
h_mean = np.mean(h_samples)
h_std = np.std(h_samples)

# 95% credible intervals
g_ci = np.percentile(g_samples, [2.5, 97.5])
h_ci = np.percentile(h_samples, [2.5, 97.5])

print("\n" + "=" * 60)
print("POSTERIOR STATISTICS (after burn-in)")
print("=" * 60)
print(f"Samples used: {len(samples_burned):,}")
print("\nGravitational acceleration (g):")
print(f"  True value:  {g_true:.4f} m/s²")
print(f"  MAP:         {g_map:.4f} m/s²")
print(f"  Posterior:   {g_mean:.4f} ± {g_std:.4f} m/s²")
print(f"  95% CI:      [{g_ci[0]:.4f}, {g_ci[1]:.4f}]")
print("\nInitial height (h):")
print(f"  True value:  {h_true:.4f} m")
print(f"  MAP:         {h_map:.4f} m")
print(f"  Posterior:   {h_mean:.4f} ± {h_std:.4f} m")
print(f"  95% CI:      [{h_ci[0]:.4f}, {h_ci[1]:.4f}]")

# Check if true values are in credible intervals
g_in_ci = g_ci[0] <= g_true <= g_ci[1]
h_in_ci = h_ci[0] <= h_true <= h_ci[1]

print("\nTrue values in 95% CI:")
print(f"  g: {'✓' if g_in_ci else '✗'}")
print(f"  h: {'✓' if h_in_ci else '✗'}")

# Scatter plot of joint distribution
fig, ax = plt.subplots(figsize=(8, 8))

# Hexbin plot for large sample sizes
hb = ax.hexbin(g_samples, h_samples, gridsize=50, cmap="Blues", mincnt=1)
ax.plot(g_true, h_true, "r*", markersize=20, label="True values", zorder=5)
ax.plot(g_map, h_map, "go", markersize=12, label="MAP estimate", zorder=5)
ax.plot(g_mean, h_mean, "mo", markersize=12, label="Posterior mean", zorder=5)

ax.set_xlabel("g (m/s²)", fontsize=12)
ax.set_ylabel("h (m)", fontsize=12)
ax.set_title("Joint Posterior Distribution", fontsize=14)
ax.legend()
ax.grid(True, alpha=0.3)

plt.colorbar(hb, ax=ax, label="Sample density")
plt.tight_layout()
plt.show()

# Calculate correlation
correlation = np.corrcoef(g_samples, h_samples)[0, 1]
print(f"\nParameter correlation: {correlation:.4f}")

if abs(correlation) > 0.7:
    print("  High correlation - parameters are not independently identifiable")
elif abs(correlation) > 0.3:
    print("  Moderate correlation")
else:
    print("  Low correlation - parameters are nearly independent")

# Take random subset of samples for predictions
n_pred_samples = 200
sample_indices = np.random.choice(len(samples_burned), n_pred_samples, replace=False)

# Generate predictions for each sample
predictions_height = []
predictions_velocity = []

for idx in sample_indices:
    g_sample, h_sample = samples_burned[idx]
    h_pred, v_pred = ball_states(t_span, g_sample, h_sample)
    predictions_height.append(h_pred)
    predictions_velocity.append(v_pred)

predictions_height = np.array(predictions_height)
predictions_velocity = np.array(predictions_velocity)

# Calculate percentiles
h_median = np.median(predictions_height, axis=0)
h_lower = np.percentile(predictions_height, 2.5, axis=0)
h_upper = np.percentile(predictions_height, 97.5, axis=0)

v_median = np.median(predictions_velocity, axis=0)
v_lower = np.percentile(predictions_velocity, 2.5, axis=0)
v_upper = np.percentile(predictions_velocity, 97.5, axis=0)

# Plot
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Height
axes[0].fill_between(t_span, h_lower, h_upper, alpha=0.3, label="95% CI")
axes[0].plot(t_span, h_median, "-", linewidth=2, label="Median prediction")
axes[0].plot(t_span, height, "--", linewidth=2, label="True", alpha=0.7)
axes[0].plot(t_span, height_noisy, "o", label="Observed", alpha=0.4, markersize=4)
axes[0].set_xlabel("Time (s)")
axes[0].set_ylabel("Height (m)")
axes[0].set_title("Posterior Predictive: Height")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Velocity
axes[1].fill_between(t_span, v_lower, v_upper, alpha=0.3, label="95% CI")
axes[1].plot(t_span, v_median, "-", linewidth=2, label="Median prediction")
axes[1].plot(t_span, velocity, "--", linewidth=2, label="True", alpha=0.7)
axes[1].plot(t_span, velocity_noisy, "o", label="Observed", alpha=0.4, markersize=4)
axes[1].set_xlabel("Time (s)")
axes[1].set_ylabel("Velocity (m/s)")
axes[1].set_title("Posterior Predictive: Velocity")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("The shaded region shows the 95% credible interval for predictions")
print("This accounts for both parameter uncertainty and observation noise")