# Import Libraries
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from scipy.stats import zscore

# Configurasi plot
sns.set(style="whitegrid")

# Data Input
data = {
    'No': list(range(1, 151)),
    'VD': [8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 6, 4, 5, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9, 5, 6, 4, 7, 8, 3, 9],
    'VP': [5, 8, 7, 3, 9, 6, 5, 4, 7, 5, 4, 9, 8, 3, 6, 4, 8, 3, 7, 5, 6, 9, 4, 6, 7, 5, 8, 9, 3, 7, 6, 4, 9, 8, 3, 5, 6, 4, 7, 8, 3, 9, 6, 5, 8, 4, 7, 9, 6, 5, 8, 4, 3, 7, 5, 6, 9, 4, 5, 7, 9, 3, 6, 5, 7, 8, 4, 5, 9, 3, 6, 4, 5, 7, 8, 3, 6, 9, 4, 5, 3, 7, 8, 9, 4, 5, 7, 6, 4, 9, 3, 5, 7, 8, 3, 5, 9, 6, 4, 5, 7, 8, 4, 9, 6, 5, 7, 8, 4, 7, 9, 6, 5, 4, 7, 8, 4, 9, 5, 6, 4, 7, 3, 4, 9, 8, 3, 7, 6, 4, 9, 5, 8, 3, 7, 6, 4, 9, 3, 5, 4, 9, 6, 3, 8, 7, 5, 4, 9, 7],
    'VS': [7, 6, 4, 8, 5, 3, 9, 6, 5, 8, 3, 6, 4, 9, 5, 7, 6, 4, 9, 8, 3, 5, 8, 7, 4, 9, 3, 6, 4, 5, 8, 3, 7, 5, 9, 4, 7, 6, 3, 9, 5, 4, 7, 9, 3, 6, 5, 8, 3, 7, 6, 5, 9, 4, 8, 3, 5, 7, 6, 3, 8, 5, 4, 9, 6, 3, 7, 8, 6, 5, 9, 7, 6, 3, 4, 9, 5, 3, 8, 7, 6, 4, 5, 3, 6, 9, 4, 3, 8, 5, 6, 4, 9, 6, 4, 7, 3, 5, 9, 6, 4, 3, 7, 5, 3, 4, 9, 6, 3, 5, 8, 3, 9, 6, 5, 3, 7, 6, 3, 9, 5, 6, 8, 7, 5, 6, 4, 9, 5, 8, 7, 3, 6, 5, 4, 8, 8, 5, 7, 6, 3, 5, 7, 4, 9, 6, 3, 9, 6, 5]
}

df = pd.DataFrame(data)
df.set_index('No', inplace=True)

# Person Fit Analysis
# Menghitung Infit dan Outfit
df['Infit_MNSQ'] = (df['VD'] - df['VD'].mean())**2 / df['VD'].var()
df['Outfit_MNSQ'] = (df['VP'] - df['VP'].mean())**2 / df['VP'].var()

# Visualisasi
plt.figure(figsize=(12, 6))
sns.scatterplot(data=df, x='Infit_MNSQ', y='Outfit_MNSQ', hue='VS', palette='coolwarm', legend='full')
plt.axhline(1.0, color='red', linestyle='--', label='Ideal Fit')
plt.axvline(1.0, color='red', linestyle='--')
plt.title('Person Fit Analysis')
plt.xlabel('Infit Mean Square')
plt.ylabel('Outfit Mean Square')
plt.legend(title='VS Score')
plt.show()

# Differential Item Functioning (DIF) Analysis
# Membagi data berdasarkan median
median_vd = df['VD'].median()
df['Group'] = np.where(df['VD'] >= median_vd, 'High VD', 'Low VD')

# Menghitung rata-rata setiap kelompok
group_means = df.groupby('Group').mean()

# Visualisasi
plt.figure(figsize=(12, 6))
sns.barplot(x=group_means.index, y=group_means['VP'], palette='muted')
plt.title('Differential Item Functioning Analysis')
plt.xlabel('VD Group')
plt.ylabel('Average VP Score')
plt.show()

# Test Information Function (TIF)
# Menghitung informasi tes
abilities = np.linspace(df['VS'].min(), df['VS'].max(), 100)
test_information = []

for ability in abilities:
    information = np.sum((ability - df['VS'])**2)
    test_information.append(information)

# Visualisasi TIF
plt.figure(figsize=(12, 6))
plt.plot(abilities, test_information, color='blue', label='Test Information')
plt.title('Test Information Function')
plt.xlabel('Ability')
plt.ylabel('Information')
plt.legend()
plt.show()

# Person-Item Map
plt.figure(figsize=(12, 8))
plt.scatter(df['VD'], df['VS'], color='green', label='Items')
plt.scatter(df['VP'], df['VS'], color='blue', label='Persons')
plt.axhline(y=df['VS'].mean(), color='red', linestyle='--', label='Average VS')
plt.title('Person-Item Map')
plt.xlabel('VD & VP Scores')
plt.ylabel('VS Scores')
plt.legend()
plt.show()

# Histogram Distribution
plt.figure(figsize=(18, 6))
plt.subplot(1, 3, 1)
sns.histplot(df['VD'], bins=10, kde=True, color='blue')
plt.title('Histogram of VD')

plt.subplot(1, 3, 2)
sns.histplot(df['VP'], bins=10, kde=True, color='green')
plt.title('Histogram of VP')

plt.subplot(1, 3, 3)
sns.histplot(df['VS'], bins=10, kde=True, color='red')
plt.title('Histogram of VS')

plt.tight_layout()
plt.show()

# Person Separation Index (PSI)
person_variance = df[['VD', 'VP', 'VS']].var().mean()
error_variance = (df['Infit_MNSQ'] + df['Outfit_MNSQ']).var() / 2
psi = (person_variance - error_variance) / person_variance

print(f'Person Separation Index (PSI): {psi:.2f}')

# Rasch Reliability
rasch_reliability = 1 - (error_variance / person_variance)

print(f'Rasch Reliability: {rasch_reliability:.2f}')