import pymc as pm
import numpy as np

def gdemo(x_obs, y_obs):
    """
    A simple Normal model with unknown mean and variance.

    Parameters:
    -----------
    x_obs : float or array-like
        First observed data point
    y_obs : float or array-like
        Second observed data point

    Returns:
    --------
    model : pm.Model
        PyMC model object ready for inference
    """
    with pm.Model() as model:
        # Priors
        s = pm.InverseGamma("s", alpha=2, beta=3)  # Variance parameter
        m = pm.Normal("m", mu=5, sigma=1)  # Mean parameter

        # Likelihood
        x = pm.Normal("x", mu=m, sigma=pm.math.sqrt(s), observed=x_obs)
        y = pm.Normal("y", mu=m, sigma=pm.math.sqrt(s), observed=y_obs)

    return model

# Example usage:
if __name__ == "__main__":
    # Example observed data
    x_obs = 4.0
    y_obs = 6.0

    # Create model
    model = gdemo(x_obs, y_obs)

    # Run inference
    with model:
        # Using NUTS (No U-Turn Sampler)
        trace = pm.sample(1000, tune=1000)

    # Summarize posterior
    print(pm.summary(trace, var_names=["s", "m"]))

    # Plot posterior distributions
    pm.plot_posterior(trace, var_names=["s", "m"])