help@rskworld.in +91 93305 39277
RSK World
  • Home
  • Development
    • Web Development
    • Mobile Apps
    • Software
    • Games
    • Project
  • Technologies
    • Data Science
    • AI Development
    • Cloud Development
    • Blockchain
    • Cyber Security
    • Dev Tools
    • Testing Tools
  • Blog
  • About
  • Contact

Theme Settings

Color Scheme
Display Options
Font Size
100%
Back to Project
RSK World
statsmodels-statistical
RSK World
statsmodels-statistical
Statistical Modeling with Statsmodels
statsmodels-statistical
  • __pycache__
  • data
  • examples
  • notebooks
  • .gitignore458 B
  • CHANGELOG.md4 KB
  • FEATURES.md6.3 KB
  • LICENSE1.2 KB
  • PROJECT_INFO.md2.2 KB
  • PROJECT_SUMMARY.md4.2 KB
  • README.md7.4 KB
  • RELEASE_NOTES_v1.0.0.md6.5 KB
  • UNIQUE_FEATURES.md5.3 KB
  • advanced_time_series.py9.8 KB
  • automated_reporting.py8.3 KB
  • bayesian_statistics.py7.5 KB
  • data_preprocessing.py8.2 KB
  • econometric_modeling.py9.8 KB
  • hypothesis_testing.py12.5 KB
  • index.html10.8 KB
  • model_evaluation.py9.1 KB
  • model_persistence.py6.5 KB
  • model_selection.py9.7 KB
  • panel_data_analysis.py7.3 KB
  • performance_benchmarking.py7.3 KB
  • regression_analysis.py9 KB
  • requirements.txt361 B
  • statistical_diagnostics.py13.8 KB
  • statsmodels-statistical.png284 B
  • time_series_analysis.py10.3 KB
  • visualization_utils.py8.9 KB
bayesian_statistics.pyperformance_benchmarking.pyconfig.yamlREADME.mdstatistical_diagnostics.py
bayesian_statistics.py
Raw Download
Find: Go to:
"""
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"
    )

252 lines•7.5 KB
python
performance_benchmarking.py
Raw Download
Find: Go to:
"""
Performance Benchmarking and Profiling

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

import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from functools import wraps
import warnings
warnings.filterwarnings('ignore')


class PerformanceBenchmark:
    """
    Performance Benchmarking Tools
    
    Author: RSK World
    Website: https://rskworld.in
    Email: help@rskworld.in
    Phone: +91 93305 39277
    """
    
    def __init__(self):
        self.benchmarks = []
    
    def benchmark_function(self, func, *args, **kwargs):
        """
        Benchmark a function execution
        
        Parameters:
        -----------
        func : callable
            Function to benchmark
        *args, **kwargs
            Arguments to pass to function
        """
        start_time = time.time()
        start_memory = self._get_memory_usage()
        
        result = func(*args, **kwargs)
        
        end_time = time.time()
        end_memory = self._get_memory_usage()
        
        execution_time = end_time - start_time
        memory_used = end_memory - start_memory
        
        benchmark = {
            'function': func.__name__,
            'execution_time': execution_time,
            'memory_used': memory_used,
            'timestamp': time.time()
        }
        
        self.benchmarks.append(benchmark)
        
        print(f"Benchmark: {func.__name__}")
        print(f"Execution Time: {execution_time:.4f} seconds")
        print(f"Memory Used: {memory_used:.2f} MB")
        
        return result, benchmark
    
    def compare_models(self, models_dict, X, y):
        """
        Compare performance of multiple models
        
        Parameters:
        -----------
        models_dict : dict
            Dictionary of model names and model objects
        X : array-like
            Independent variables
        y : array-like
            Dependent variable
        """
        results = []
        
        for name, model in models_dict.items():
            print(f"\nBenchmarking: {name}")
            print("-" * 70)
            
            # Training time
            start = time.time()
            if hasattr(model, 'fit'):
                model.fit(X, y)
            training_time = time.time() - start
            
            # Prediction time
            start = time.time()
            if hasattr(model, 'predict'):
                predictions = model.predict(X)
            prediction_time = time.time() - start
            
            # Model metrics
            metrics = {}
            if hasattr(model, 'results'):
                metrics['aic'] = model.results.aic if hasattr(model.results, 'aic') else None
                metrics['bic'] = model.results.bic if hasattr(model.results, 'bic') else None
                metrics['r_squared'] = model.results.rsquared if hasattr(model.results, 'rsquared') else None
            
            results.append({
                'model': name,
                'training_time': training_time,
                'prediction_time': prediction_time,
                'total_time': training_time + prediction_time,
                **metrics
            })
        
        results_df = pd.DataFrame(results)
        
        print("\n" + "=" * 70)
        print("MODEL PERFORMANCE COMPARISON")
        print("=" * 70)
        print(results_df.to_string(index=False))
        
        return results_df
    
    def plot_benchmark_comparison(self, metric='execution_time'):
        """Plot benchmark comparison"""
        if not self.benchmarks:
            print("No benchmarks to plot")
            return
        
        df = pd.DataFrame(self.benchmarks)
        
        plt.figure(figsize=(10, 6))
        plt.bar(df['function'], df[metric])
        plt.xlabel('Function')
        plt.ylabel(metric.replace('_', ' ').title())
        plt.title('Performance Benchmark Comparison')
        plt.xticks(rotation=45, ha='right')
        plt.grid(True, alpha=0.3, axis='y')
        plt.tight_layout()
        plt.show()
    
    def _get_memory_usage(self):
        """Get current memory usage in MB"""
        try:
            import psutil
            process = psutil.Process()
            return process.memory_info().rss / 1024 / 1024  # Convert to MB
        except ImportError:
            return 0
    
    def timeit(self, n_iterations=10):
        """
        Decorator for timing function execution
        
        Parameters:
        -----------
        n_iterations : int
            Number of iterations to average
        """
        def decorator(func):
            @wraps(func)
            def wrapper(*args, **kwargs):
                times = []
                for _ in range(n_iterations):
                    start = time.time()
                    result = func(*args, **kwargs)
                    times.append(time.time() - start)
                
                avg_time = np.mean(times)
                std_time = np.std(times)
                
                print(f"{func.__name__}: {avg_time:.4f} ± {std_time:.4f} seconds (avg of {n_iterations} runs)")
                
                return result
            return wrapper
        return decorator


def benchmark_model_complexity(X, y, model_class, complexity_params):
    """
    Benchmark model performance across different complexities
    
    Parameters:
    -----------
    X : array-like
        Independent variables
    y : array-like
        Dependent variable
    model_class : class
        Model class to instantiate
    complexity_params : dict
        Dictionary of complexity parameter names and values
    """
    results = []
    
    for param_name, param_values in complexity_params.items():
        for param_value in param_values:
            model = model_class()
            setattr(model, param_name, param_value)
            
            start = time.time()
            model.fit(X, y)
            training_time = time.time() - start
            
            start = time.time()
            predictions = model.predict(X)
            prediction_time = time.time() - start
            
            results.append({
                param_name: param_value,
                'training_time': training_time,
                'prediction_time': prediction_time,
                'total_time': training_time + prediction_time
            })
    
    return pd.DataFrame(results)


if __name__ == "__main__":
    # Example usage
    print("Performance Benchmarking Example")
    print("=" * 70)
    
    from regression_analysis import LinearRegressionModel
    import numpy as np
    
    # Generate sample data
    np.random.seed(42)
    X = np.random.randn(1000, 10)
    y = 2 + np.sum(X[:, :3], axis=1) + np.random.randn(1000) * 0.5
    
    # Benchmark
    benchmark = PerformanceBenchmark()
    
    def fit_model():
        model = LinearRegressionModel()
        model.fit(X, y)
        return model
    
    result, bench = benchmark.benchmark_function(fit_model)
    
    # Compare models
    models = {
        'Simple Model': LinearRegressionModel()
    }
    
    comparison = benchmark.compare_models(models, X, y)

249 lines•7.3 KB
python
README.md
Raw Download

README.md

# Statsmodels Statistical Modeling

<!--
Author: RSK World
Website: https://rskworld.in
Email: help@rskworld.in
Phone: +91 93305 39277
Description: Statistical modeling with Statsmodels including regression analysis, time series models, hypothesis testing, and statistical tests.
-->

Statistical modeling with Statsmodels including regression analysis, time series models, hypothesis testing, and statistical tests.

## Description

This project demonstrates Statsmodels, a library for statistical modeling and econometrics in Python. It covers linear and generalized linear models, time series analysis, hypothesis testing, statistical tests, and diagnostic tools. Perfect for statistical analysis and econometric modeling.

## Features

- **Linear and GLM regression** - OLS, GLM with multiple families, comprehensive diagnostics
- **Time series analysis** - ARIMA, SARIMA, exponential smoothing, decomposition, forecasting
- **Advanced time series** - Auto ARIMA selection, SARIMA models, comprehensive stationarity tests
- **Hypothesis testing** - T-tests, ANOVA, chi-square, normality tests, non-parametric tests
- **Statistical diagnostics** - Multicollinearity, heteroscedasticity, autocorrelation, influential points
- **Econometric modeling** - VAR, VARMAX, cointegration tests, impulse response functions, Granger causality
- **Model selection** - Stepwise selection, model comparison, information criteria
- **Model evaluation** - Cross-validation, time series CV, multiple metrics, learning curves
- **Feature selection** - VIF-based removal, correlation filtering
- **Data preprocessing** - Missing value handling, outlier detection/removal, scaling, stationarity transformation
- **Visualization utilities** - Comprehensive plotting functions for all analyses
- **Bayesian statistics** - Bayesian inference, posterior distributions, Bayes factors
- **Panel data analysis** - Fixed effects, random effects, Hausman test
- **Model persistence** - Save/load models, model serialization, metadata management
- **Automated reporting** - Generate comprehensive reports in TXT and HTML formats
- **Performance benchmarking** - Model comparison, execution time profiling, memory usage

## Technologies

- Python 3.8+
- Statsmodels
- Pandas
- NumPy
- Matplotlib
- Seaborn
- Scikit-learn
- SciPy
- Jupyter Notebook

## Installation

```bash
pip install -r requirements.txt
```

## Usage

### Linear Regression

```python
from regression_analysis import LinearRegressionModel

# Create and fit model
model = LinearRegressionModel()
model.fit(X, y)
model.summary()
```

### Time Series Analysis

```python
from time_series_analysis import TimeSeriesModel

# Create and fit time series model
ts_model = TimeSeriesModel()
ts_model.fit(data)
ts_model.forecast(steps=10)
```

### Hypothesis Testing

```python
from hypothesis_testing import StatisticalTests

# Perform statistical tests
tests = StatisticalTests()
tests.t_test(data)
tests.chi_square_test(data)
```

### Model Selection

```python
from model_selection import ModelSelection

# Compare multiple models
selector = ModelSelection()
comparison = selector.compare_models(X, y, models_dict)

# Stepwise feature selection
features, model = selector.stepwise_selection(X, y)
```

### Model Evaluation

```python
from model_evaluation import ModelEvaluation

# Cross-validation
evaluator = ModelEvaluation()
cv_results = evaluator.cross_validate(X, y, model_func, cv_folds=5)

# Calculate metrics
metrics = evaluator.calculate_metrics(y_true, y_pred)
```

### Advanced Time Series

```python
from advanced_time_series import SARIMAModel, AutoARIMA

# SARIMA model
sarima = SARIMAModel()
sarima.fit(data, order=(1,1,1), seasonal_order=(1,1,1,12))

# Auto ARIMA selection
auto_arima = AutoARIMA()
best_model = auto_arima.auto_select(data)
```

### Data Preprocessing

```python
from data_preprocessing import DataPreprocessor

# Handle missing values and outliers
preprocessor = DataPreprocessor()
cleaned_data = preprocessor.remove_outliers(data)
scaled_data = preprocessor.scale_data(data, method='standard')
```

### Visualization

```python
from visualization_utils import StatisticalVisualizations

# Create comprehensive plots
viz = StatisticalVisualizations()
viz.plot_correlation_matrix(data)
viz.plot_residual_analysis(residuals, fitted_values)
```

### Bayesian Statistics

```python
from bayesian_statistics import BayesianAnalysis

# Bayesian t-test
result = BayesianAnalysis.bayesian_ttest(sample1, sample2)

# Bayesian linear regression
bayesian_result = BayesianAnalysis.bayesian_linear_regression(X, y)
```

### Panel Data Analysis

```python
from panel_data_analysis import PanelDataAnalysis

# Prepare and analyze panel data
panel = PanelDataAnalysis()
panel.prepare_panel_data(df, 'entity', 'time', ['X1', 'X2', 'y'])
fe_model = panel.fixed_effects_regression('y', ['X1', 'X2'])
```

### Model Persistence

```python
from model_persistence import ModelPersistence

# Save and load models
persistence = ModelPersistence()
persistence.save_model(model, 'my_model', metadata={'r_squared': 0.95})
loaded_model, metadata = persistence.load_model('saved_models/my_model.pkl')
```

### Automated Reporting

```python
from automated_reporting import AutomatedReport

# Generate comprehensive reports
reporter = AutomatedReport()
reporter.generate_regression_report(model, X, y)
reporter.save_report('analysis_report', format='html')
```

### Performance Benchmarking

```python
from performance_benchmarking import PerformanceBenchmark

# Benchmark model performance
benchmark = PerformanceBenchmark()
comparison = benchmark.compare_models(models_dict, X, y)
```

## Project Structure

```
statsmodels-statistical/
├── README.md
├── requirements.txt
├── LICENSE
├── index.html
├── regression_analysis.py # Linear and GLM regression
├── time_series_analysis.py # Basic time series models
├── advanced_time_series.py # SARIMA, Auto ARIMA
├── hypothesis_testing.py # Statistical tests
├── statistical_diagnostics.py # Model diagnostics
├── econometric_modeling.py # VAR, cointegration
├── model_selection.py # Model comparison, stepwise selection
├── model_evaluation.py # Cross-validation, metrics
├── data_preprocessing.py # Data cleaning, scaling
├── visualization_utils.py # Advanced plotting
├── bayesian_statistics.py # Bayesian inference
├── panel_data_analysis.py # Panel data models
├── model_persistence.py # Model saving/loading
├── automated_reporting.py # Report generation
├── performance_benchmarking.py # Performance profiling
├── notebooks/
│ ├── 01_linear_regression.ipynb
│ ├── 02_time_series.ipynb
│ ├── 03_hypothesis_testing.ipynb
│ └── 04_econometric_modeling.ipynb
├── data/
│ └── sample_data.csv
└── examples/
├── regression_example.py
├── time_series_example.py
└── hypothesis_testing_example.py
```

## Author

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

## License

This project is provided as educational material for statistical modeling and analysis.

statistical_diagnostics.py
Raw Download
Find: Go to:
"""
Statistical Diagnostics and Model Validation

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
import seaborn as sns
from statsmodels.stats.diagnostic import (
    het_breuschpagan, het_white, acorr_ljungbox,
    linear_harvey_collier, linear_rainbow
)
from statsmodels.stats.outliers_influence import (
    variance_inflation_factor, influence_plot
)
from statsmodels.stats.stattools import durbin_watson
from scipy import stats
import warnings
warnings.filterwarnings('ignore')


class ModelDiagnostics:
    """
    Comprehensive Model Diagnostics
    
    Author: RSK World
    Website: https://rskworld.in
    Email: help@rskworld.in
    Phone: +91 93305 39277
    """
    
    def __init__(self, model_results):
        """
        Initialize diagnostics with model results
        
        Parameters:
        -----------
        model_results : RegressionResults
            Fitted model results from statsmodels
        """
        self.results = model_results
        self.residuals = model_results.resid
        self.fitted_values = model_results.fittedvalues
    
    def check_linearity(self):
        """Check linearity assumption using Rainbow test"""
        try:
            f_stat, p_value = linear_rainbow(self.results)
            print("Rainbow Test for Linearity:")
            print(f"F-statistic: {f_stat:.4f}")
            print(f"p-value: {p_value:.4f}")
            
            if p_value < 0.05:
                print("Warning: Non-linearity detected (p < 0.05)")
            else:
                print("Linearity assumption appears valid")
            
            return {'f_statistic': f_stat, 'p_value': p_value}
        except Exception as e:
            print(f"Error in linearity test: {e}")
            return None
    
    def check_heteroscedasticity(self, test='both'):
        """
        Check for heteroscedasticity
        
        Parameters:
        -----------
        test : str
            'breusch-pagan', 'white', or 'both'
        """
        results = {}
        
        if test in ['breusch-pagan', 'both']:
            try:
                lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(
                    self.residuals, self.results.model.exog
                )
                print("\nBreusch-Pagan Test for Heteroscedasticity:")
                print(f"LM Statistic: {lm:.4f}")
                print(f"LM p-value: {lm_pvalue:.4f}")
                print(f"F Statistic: {fvalue:.4f}")
                print(f"F p-value: {f_pvalue:.4f}")
                
                if f_pvalue < 0.05:
                    print("Warning: Heteroscedasticity detected (p < 0.05)")
                else:
                    print("No significant heteroscedasticity detected")
                
                results['breusch_pagan'] = {
                    'lm': lm, 'lm_pvalue': lm_pvalue,
                    'fvalue': fvalue, 'f_pvalue': f_pvalue
                }
            except Exception as e:
                print(f"Error in Breusch-Pagan test: {e}")
        
        if test in ['white', 'both']:
            try:
                lm, lm_pvalue, fvalue, f_pvalue = het_white(
                    self.residuals, self.results.model.exog
                )
                print("\nWhite Test for Heteroscedasticity:")
                print(f"LM Statistic: {lm:.4f}")
                print(f"LM p-value: {lm_pvalue:.4f}")
                print(f"F Statistic: {fvalue:.4f}")
                print(f"F p-value: {f_pvalue:.4f}")
                
                if f_pvalue < 0.05:
                    print("Warning: Heteroscedasticity detected (p < 0.05)")
                else:
                    print("No significant heteroscedasticity detected")
                
                results['white'] = {
                    'lm': lm, 'lm_pvalue': lm_pvalue,
                    'fvalue': fvalue, 'f_pvalue': f_pvalue
                }
            except Exception as e:
                print(f"Error in White test: {e}")
        
        return results
    
    def check_autocorrelation(self, lags=10):
        """
        Check for autocorrelation
        
        Parameters:
        -----------
        lags : int
            Number of lags to test
        """
        # Durbin-Watson test
        dw = durbin_watson(self.residuals)
        print("\nDurbin-Watson Test for Autocorrelation:")
        print(f"Durbin-Watson Statistic: {dw:.4f}")
        
        if dw < 1.5:
            print("Warning: Positive autocorrelation detected")
        elif dw > 2.5:
            print("Warning: Negative autocorrelation detected")
        else:
            print("No significant autocorrelation detected")
        
        # Ljung-Box test
        lb_test = acorr_ljungbox(self.residuals, lags=lags, return_df=True)
        print(f"\nLjung-Box Test (lags={lags}):")
        print(lb_test)
        
        if lb_test['lb_pvalue'].iloc[-1] < 0.05:
            print("Warning: Residual autocorrelation detected (p < 0.05)")
        else:
            print("No significant residual autocorrelation detected")
        
        return {'durbin_watson': dw, 'ljung_box': lb_test}
    
    def check_multicollinearity(self):
        """Check for multicollinearity using VIF"""
        try:
            exog = self.results.model.exog
            vif_data = pd.DataFrame()
            vif_data["Variable"] = range(exog.shape[1])
            vif_data["VIF"] = [variance_inflation_factor(exog, i) 
                              for i in range(exog.shape[1])]
            
            print("\nVariance Inflation Factor (VIF):")
            print(vif_data)
            print("\nVIF > 10 indicates multicollinearity")
            
            high_vif = vif_data[vif_data['VIF'] > 10]
            if len(high_vif) > 0:
                print(f"\nWarning: {len(high_vif)} variable(s) with VIF > 10")
            else:
                print("\nNo significant multicollinearity detected")
            
            return vif_data
        except Exception as e:
            print(f"Error in VIF calculation: {e}")
            return None
    
    def check_normality(self):
        """Check normality of residuals"""
        print("\nNormality Tests for Residuals:")
        print("=" * 50)
        
        # Shapiro-Wilk test
        shapiro_stat, shapiro_p = stats.shapiro(self.residuals)
        print(f"\nShapiro-Wilk Test:")
        print(f"Statistic: {shapiro_stat:.4f}")
        print(f"p-value: {shapiro_p:.4f}")
        
        # Jarque-Bera test
        from statsmodels.stats.diagnostic import jarque_bera
        jb_stat, jb_p, _, _ = jarque_bera(self.residuals)
        print(f"\nJarque-Bera Test:")
        print(f"Statistic: {jb_stat:.4f}")
        print(f"p-value: {jb_p:.4f}")
        
        alpha = 0.05
        is_normal = shapiro_p > alpha and jb_p > alpha
        
        if not is_normal:
            print("\nWarning: Residuals may not be normally distributed")
        else:
            print("\nResiduals appear to be normally distributed")
        
        return {
            'shapiro': {'statistic': shapiro_stat, 'p_value': shapiro_p},
            'jarque_bera': {'statistic': jb_stat, 'p_value': jb_p},
            'is_normal': is_normal
        }
    
    def check_influential_points(self, threshold=2):
        """
        Check for influential observations
        
        Parameters:
        -----------
        threshold : float
            Threshold for Cook's distance
        """
        try:
            influence = self.results.get_influence()
            cooks_d = influence.cooks_distance[0]
            
            print(f"\nInfluential Points Analysis (Cook's Distance > {threshold}):")
            influential = np.where(cooks_d > threshold)[0]
            
            if len(influential) > 0:
                print(f"Found {len(influential)} influential observation(s):")
                print(f"Indices: {influential}")
                print(f"Cook's D values: {cooks_d[influential]}")
            else:
                print("No influential observations detected")
            
            # Plot influence
            fig, axes = plt.subplots(1, 2, figsize=(12, 5))
            
            axes[0].stem(range(len(cooks_d)), cooks_d)
            axes[0].axhline(y=threshold, color='r', linestyle='--', 
                          label=f'Threshold ({threshold})')
            axes[0].set_xlabel('Observation')
            axes[0].set_ylabel("Cook's Distance")
            axes[0].set_title("Cook's Distance")
            axes[0].legend()
            axes[0].grid(True, alpha=0.3)
            
            # Leverage plot
            leverage = influence.hat_matrix_diag
            axes[1].scatter(leverage, cooks_d, alpha=0.6)
            axes[1].axhline(y=threshold, color='r', linestyle='--')
            axes[1].set_xlabel('Leverage')
            axes[1].set_ylabel("Cook's Distance")
            axes[1].set_title('Influence Plot')
            axes[1].grid(True, alpha=0.3)
            
            plt.tight_layout()
            plt.show()
            
            return {
                'cooks_distance': cooks_d,
                'influential_indices': influential,
                'leverage': leverage
            }
        except Exception as e:
            print(f"Error in influential points analysis: {e}")
            return None
    
    def comprehensive_diagnostics(self):
        """Run all diagnostic tests"""
        print("=" * 70)
        print("COMPREHENSIVE MODEL DIAGNOSTICS")
        print("=" * 70)
        
        diagnostics = {}
        
        # Linearity
        print("\n" + "-" * 70)
        diagnostics['linearity'] = self.check_linearity()
        
        # Heteroscedasticity
        print("\n" + "-" * 70)
        diagnostics['heteroscedasticity'] = self.check_heteroscedasticity(test='both')
        
        # Autocorrelation
        print("\n" + "-" * 70)
        diagnostics['autocorrelation'] = self.check_autocorrelation()
        
        # Multicollinearity
        print("\n" + "-" * 70)
        diagnostics['multicollinearity'] = self.check_multicollinearity()
        
        # Normality
        print("\n" + "-" * 70)
        diagnostics['normality'] = self.check_normality()
        
        # Influential points
        print("\n" + "-" * 70)
        diagnostics['influential_points'] = self.check_influential_points()
        
        return diagnostics
    
    def plot_diagnostics(self):
        """Create comprehensive diagnostic plots"""
        fig = plt.figure(figsize=(15, 10))
        
        # Residuals vs Fitted
        ax1 = plt.subplot(2, 3, 1)
        ax1.scatter(self.fitted_values, self.residuals, alpha=0.6)
        ax1.axhline(y=0, color='r', linestyle='--')
        ax1.set_xlabel('Fitted Values')
        ax1.set_ylabel('Residuals')
        ax1.set_title('Residuals vs Fitted')
        ax1.grid(True, alpha=0.3)
        
        # Q-Q Plot
        ax2 = plt.subplot(2, 3, 2)
        stats.probplot(self.residuals, dist="norm", plot=ax2)
        ax2.set_title('Q-Q Plot')
        ax2.grid(True, alpha=0.3)
        
        # Scale-Location Plot
        ax3 = plt.subplot(2, 3, 3)
        sqrt_abs_residuals = np.sqrt(np.abs(self.residuals))
        ax3.scatter(self.fitted_values, sqrt_abs_residuals, alpha=0.6)
        ax3.set_xlabel('Fitted Values')
        ax3.set_ylabel('√|Standardized Residuals|')
        ax3.set_title('Scale-Location Plot')
        ax3.grid(True, alpha=0.3)
        
        # Residuals Histogram
        ax4 = plt.subplot(2, 3, 4)
        ax4.hist(self.residuals, bins=30, edgecolor='black', alpha=0.7)
        ax4.set_xlabel('Residuals')
        ax4.set_ylabel('Frequency')
        ax4.set_title('Residuals Distribution')
        ax4.grid(True, alpha=0.3)
        
        # Leverage Plot
        ax5 = plt.subplot(2, 3, 5)
        try:
            influence = self.results.get_influence()
            leverage = influence.hat_matrix_diag
            cooks_d = influence.cooks_distance[0]
            ax5.scatter(leverage, cooks_d, alpha=0.6)
            ax5.set_xlabel('Leverage')
            ax5.set_ylabel("Cook's Distance")
            ax5.set_title('Influence Plot')
            ax5.grid(True, alpha=0.3)
        except:
            ax5.text(0.5, 0.5, 'Leverage data not available', 
                    ha='center', va='center')
        
        # ACF of Residuals
        ax6 = plt.subplot(2, 3, 6)
        from statsmodels.tsa.stattools import acf
        acf_values = acf(self.residuals, nlags=20)
        ax6.stem(range(len(acf_values)), acf_values)
        ax6.axhline(y=1.96/np.sqrt(len(self.residuals)), 
              color='r', linestyle='--', alpha=0.5)
        ax6.axhline(y=-1.96/np.sqrt(len(self.residuals)), 
                   color='r', linestyle='--', alpha=0.5)
        ax6.set_xlabel('Lag')
        ax6.set_ylabel('ACF')
        ax6.set_title('ACF of Residuals')
        ax6.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()


if __name__ == "__main__":
    # Example usage
    print("Statistical Diagnostics Example")
    print("=" * 50)
    
    from regression_analysis import LinearRegressionModel
    
    # Generate sample data
    np.random.seed(42)
    n = 100
    X = np.random.randn(n, 3)
    y = 2 + 1.5 * X[:, 0] + 0.8 * X[:, 1] - 0.5 * X[:, 2] + np.random.randn(n) * 0.5
    
    # Fit model
    model = LinearRegressionModel()
    model.fit(X, y)
    
    # Create diagnostics
    diagnostics = ModelDiagnostics(model.results)
    
    # Run comprehensive diagnostics
    diagnostics.comprehensive_diagnostics()
    
    # Plot diagnostics
    diagnostics.plot_diagnostics()

401 lines•13.8 KB
python
🚀 Support RSK World

Subscribe to our YouTube channel for latest tutorials & updates!



Click subscribe & support our work ❤️

About RSK World

Founded by Molla Samser, with Designer & Tester Rima Khatun, RSK World is your one-stop destination for free programming resources, source code, and development tools.

Founder: Molla Samser
Designer & Tester: Rima Khatun

Development

  • Game Development
  • Web Development
  • Mobile Development
  • AI Development
  • Development Tools

Legal

  • Terms & Conditions
  • Privacy Policy
  • Disclaimer

Contact Info

Nutanhat, Mongolkote
Purba Burdwan, West Bengal
India, 713147

+91 93305 39277

hello@rskworld.in
support@rskworld.in

© 2026 RSK World. All rights reserved.

Content used for educational purposes only. View Disclaimer