"""
Bayesian Statistical Analysis

Author: RSK World
Website: https://rskworld.in
Email: help@rskworld.in
Phone: +91 93305 39277
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')


class BayesianAnalysis:
    """
    Bayesian Statistical Analysis Tools
    
    Author: RSK World
    Website: https://rskworld.in
    Email: help@rskworld.in
    Phone: +91 93305 39277
    """
    
    @staticmethod
    def bayesian_ttest(sample1, sample2, prior_mean=0, prior_std=1):
        """
        Bayesian t-test
        
        Parameters:
        -----------
        sample1 : array-like
            First sample
        sample2 : array-like
            Second sample
        prior_mean : float
            Prior mean
        prior_std : float
            Prior standard deviation
        """
        n1, n2 = len(sample1), len(sample2)
        mean1, mean2 = np.mean(sample1), np.mean(sample2)
        var1, var2 = np.var(sample1, ddof=1), np.var(sample2, ddof=1)
        
        # Pooled variance
        pooled_var = ((n1 - 1) * var1 + (n2 - 1) * var2) / (n1 + n2 - 2)
        pooled_std = np.sqrt(pooled_var)
        
        # Standard error
        se = pooled_std * np.sqrt(1/n1 + 1/n2)
        
        # Posterior distribution parameters
        posterior_mean = (prior_mean / prior_std**2 + (mean1 - mean2) / se**2) / \
                        (1 / prior_std**2 + 1 / se**2)
        posterior_var = 1 / (1 / prior_std**2 + 1 / se**2)
        posterior_std = np.sqrt(posterior_var)
        
        # Credible interval (95%)
        ci_lower = posterior_mean - 1.96 * posterior_std
        ci_upper = posterior_mean + 1.96 * posterior_std
        
        print("Bayesian t-test Results:")
        print("=" * 70)
        print(f"Posterior Mean: {posterior_mean:.4f}")
        print(f"Posterior Std: {posterior_std:.4f}")
        print(f"95% Credible Interval: [{ci_lower:.4f}, {ci_upper:.4f}]")
        
        # Probability that difference > 0
        prob_positive = 1 - stats.norm.cdf(0, posterior_mean, posterior_std)
        print(f"Probability (difference > 0): {prob_positive:.4f}")
        
        return {
            'posterior_mean': posterior_mean,
            'posterior_std': posterior_std,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper,
            'prob_positive': prob_positive
        }
    
    @staticmethod
    def bayesian_linear_regression(X, y, prior_precision=1.0):
        """
        Bayesian Linear Regression
        
        Parameters:
        -----------
        X : array-like
            Independent variables
        y : array-like
            Dependent variable
        prior_precision : float
            Prior precision (inverse variance)
        """
        from statsmodels.api import add_constant, OLS
        
        X_with_const = add_constant(X)
        n, p = X_with_const.shape
        
        # OLS estimates
        ols_model = OLS(y, X_with_const).fit()
        beta_ols = ols_model.params
        sigma_sq = ols_model.mse_resid
        
        # Prior covariance
        prior_cov = np.eye(p) / prior_precision
        
        # Posterior covariance
        XtX = X_with_const.T @ X_with_const
        posterior_cov = np.linalg.inv(prior_precision * np.eye(p) + XtX / sigma_sq)
        
        # Posterior mean
        posterior_mean = posterior_cov @ (XtX @ beta_ols / sigma_sq)
        
        # Credible intervals
        posterior_std = np.sqrt(np.diag(posterior_cov))
        ci_lower = posterior_mean - 1.96 * posterior_std
        ci_upper = posterior_mean + 1.96 * posterior_std
        
        print("Bayesian Linear Regression Results:")
        print("=" * 70)
        results_df = pd.DataFrame({
            'Posterior Mean': posterior_mean,
            'Posterior Std': posterior_std,
            'CI Lower': ci_lower,
            'CI Upper': ci_upper
        }, index=[f'beta_{i}' for i in range(p)])
        print(results_df)
        
        return {
            'posterior_mean': posterior_mean,
            'posterior_cov': posterior_cov,
            'posterior_std': posterior_std,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper
        }
    
    @staticmethod
    def plot_posterior_distribution(mean, std, true_value=None, title="Posterior Distribution"):
        """Plot posterior distribution"""
        x = np.linspace(mean - 4*std, mean + 4*std, 1000)
        y = stats.norm.pdf(x, mean, std)
        
        plt.figure(figsize=(10, 6))
        plt.plot(x, y, 'b-', linewidth=2, label='Posterior')
        plt.axvline(mean, color='r', linestyle='--', label=f'Mean: {mean:.4f}')
        
        if true_value is not None:
            plt.axvline(true_value, color='g', linestyle='--', label=f'True: {true_value:.4f}')
        
        # Credible interval
        ci_lower = mean - 1.96 * std
        ci_upper = mean + 1.96 * std
        plt.axvspan(ci_lower, ci_upper, alpha=0.2, color='blue', label='95% CI')
        
        plt.xlabel('Parameter Value')
        plt.ylabel('Density')
        plt.title(title)
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()


class BayesianModelComparison:
    """
    Bayesian Model Comparison
    
    Author: RSK World
    Website: https://rskworld.in
    Email: help@rskworld.in
    Phone: +91 93305 39277
    """
    
    @staticmethod
    def bayes_factor(model1_loglik, model2_loglik):
        """
        Calculate Bayes Factor
        
        Parameters:
        -----------
        model1_loglik : float
            Log-likelihood of model 1
        model2_loglik : float
            Log-likelihood of model 2
        """
        bf = np.exp(model1_loglik - model2_loglik)
        
        print("Bayes Factor:")
        print("=" * 70)
        print(f"BF (Model 1 / Model 2): {bf:.4f}")
        
        if bf > 1:
            print(f"Model 1 is {bf:.2f}x more likely than Model 2")
        else:
            print(f"Model 2 is {1/bf:.2f}x more likely than Model 1")
        
        # Interpretation
        if bf > 100:
            strength = "Very strong evidence"
        elif bf > 10:
            strength = "Strong evidence"
        elif bf > 3:
            strength = "Moderate evidence"
        else:
            strength = "Weak evidence"
        
        print(f"Evidence strength: {strength}")
        
        return bf
    
    @staticmethod
    def bayesian_information_criterion(loglik, n_params, n_obs):
        """
        Calculate Bayesian Information Criterion (BIC)
        
        Parameters:
        -----------
        loglik : float
            Log-likelihood
        n_params : int
            Number of parameters
        n_obs : int
            Number of observations
        """
        bic = -2 * loglik + n_params * np.log(n_obs)
        return bic


if __name__ == "__main__":
    # Example usage
    print("Bayesian Statistics Example")
    print("=" * 70)
    
    # Generate sample data
    np.random.seed(42)
    sample1 = np.random.normal(100, 15, 30)
    sample2 = np.random.normal(105, 15, 30)
    
    # Bayesian t-test
    result = BayesianAnalysis.bayesian_ttest(sample1, sample2)
    
    # Plot posterior
    BayesianAnalysis.plot_posterior_distribution(
        result['posterior_mean'], 
        result['posterior_std'],
        title="Posterior Distribution of Mean Difference"
    )