# Complex demo: Load California Housing, simulate probability distributions (normal, binomial), compute Bayes' theorem for price probability, visualize with SciPy/NumPy/Matplotlib, integrate with regression
import pandas as pd
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from scipy.stats import norm, binom
from scipy.special import logsumexp
import matplotlib.pyplot as plt
import seaborn as sns

def advanced_housing_probability_demo():
    # Load and prep
    housing = fetch_california_housing()
    df = pd.DataFrame(housing.data, columns=housing.feature_names)
    df['Target'] = housing.target
    # Feature engineering
    df['RoomsPerPop'] = df['AveRooms'] / (df['Population'] + 1e-5)
    
    # Simulate probability distributions
    # Normal distribution for median income (mean=3.87, std=1.9)
    x = np.linspace(0, 10, 100)
    y_normal = norm.pdf(x, loc=3.87, scale=1.9)
    # Binomial distribution for simulated price category (n=10, p=0.5)
    y_binomial = binom.pmf(np.arange(0, 11), n=10, p=0.5)
    # Bayes' theorem simulation for price probability
    prior = np.array([0.5, 0.5])  # Prior P(High/Low Price)
    likelihood = np.array([0.6, 0.4])  # Simulated P(Feature|Price)
    evidence = np.sum(likelihood * prior)
    posterior = (likelihood * prior) / evidence
    print(f"Posterior Probabilities for Price: {posterior}")
    
    # Regression integration
    X = df.drop('Target', axis=1)
    y = df['Target']
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)
    model = LinearRegression()
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    mse = mean_squared_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    print(f"MSE: {mse:.2f}, R2: {r2:.2f}")
    
    # Visualize distributions
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    axes[0].plot(x, y_normal, label='Normal PDF (Median Income)')
    axes[0].set_title('Probability Distribution for Housing Feature')
    axes[0].set_xlabel('Value')
    axes[0].set_ylabel('Probability Density')
    axes[0].annotate('Peak at Mean', xy=(3.87, norm.pdf(3.87, loc=3.87, scale=1.9)), xytext=(4, 0.2), arrowprops=dict(facecolor='black', shrink=0.05))
    axes[0].legend()
    axes[1].bar(np.arange(len(posterior)), posterior, label='Bayes Posterior')
    axes[1].set_title('Bayes\' Theorem for Price Probability')
    axes[1].set_xlabel('Price Category')
    axes[1].set_ylabel('Posterior Probability')
    axes[1].annotate('High Prob Category', xy=(0, posterior[0]), xytext=(1, posterior[0] + 0.05), arrowprops=dict(facecolor='black', shrink=0.05))
    axes[1].legend()
    plt.tight_layout()
    plt.show()
    
    # Visualize predictions
    plt.figure(figsize=(8, 5))
    plt.scatter(y_test, y_pred, c='blue', alpha=0.5)
    plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'k--', lw=2)
    plt.title('Housing Regression Predictions')
    plt.xlabel('True Values')
    plt.ylabel('Predictions')
    plt.annotate('Best Fit', xy=(y_test.mean(), y_pred.mean()), xytext=(y_test.mean() + 0.5, y_pred.mean() + 0.5), arrowprops=dict(facecolor='black', shrink=0.05))
    plt.show()

advanced_housing_probability_demo()