Phân phối chuẩn hóa – Normality

* Kiểm tra giả định phân phối chuẩn trong SPSS
* Giả sử biến cần kiểm tra là 'income' (thu nhập)

* Phương pháp 1: Tạo histogram với đường cong chuẩn
FREQUENCIES VARIABLES=income
  /FORMAT=NOTABLE
  /HISTOGRAM=NORMAL
  /ORDER=ANALYSIS.

* Phương pháp 2: Tính toán Skewness và Kurtosis
DESCRIPTIVES VARIABLES=income
  /STATISTICS=MEAN STDDEV VARIANCE SKEWNESS KURTOSIS.

* Phương pháp 3: Kiểm định K-S và Shapiro-Wilk
EXAMINE VARIABLES=income
  /PLOT=HISTOGRAM NPPLOT
  /STATISTICS=DESCRIPTIVES
  /CINTERVAL=95
  /MISSING=LISTWISE
  /NOTOTAL.

* Phương pháp 4: Kiểm tra chi tiết với Q-Q plot
PPLOT
  /VARIABLES=income
  /NOLOG
  /NOSTANDARDIZE
  /TYPE=Q-Q
  /FRACTION=BLOM
  /TIES=MEAN
  /DIST=NORMAL.

* Ví dụ biến đổi dữ liệu nếu không có phân phối chuẩn
* Biến đổi logarithm
COMPUTE income_log = LN(income).
EXECUTE.

* Biến đổi căn bậc hai
COMPUTE income_sqrt = SQRT(income).
EXECUTE.

* Biến đổi nghịch đảo
COMPUTE income_reciprocal = 1/income.
EXECUTE.

* Kiểm tra lại phân phối chuẩn sau biến đổi
EXAMINE VARIABLES=income_log income_sqrt income_reciprocal
  /PLOT=HISTOGRAM NPPLOT
  /STATISTICS=DESCRIPTIVES
  /CINTERVAL=95
  /MISSING=LISTWISE
  /NOTOTAL.

* Kiểm tra giả định phân phối chuẩn trong Stata
* Giả sử biến cần kiểm tra là 'income' (thu nhập)

* Mở dữ liệu
use "data.dta", clear

* Phương pháp 1: Tạo histogram với đường cong chuẩn
histogram income, normal ///
    title("Histogram của thu nhập với đường cong chuẩn") ///
    xtitle("Thu nhập") ytitle("Mật độ") ///
    scheme(s1mono)

* Phương pháp 2: Tính toán thống kê mô tả
summarize income, detail

* Tính toán Skewness và Kurtosis chi tiết
sktest income

* Phương pháp 3: Kiểm định Shapiro-Wilk
swilk income

* Kiểm định Kolmogorov-Smirnov
ksmirnov income = normal((income-r(mean))/r(sd))

* Tạo Q-Q plot
qnorm income, ///
    title("Q-Q Plot cho thu nhập") ///
    ytitle("Phân vị mẫu") xtitle("Phân vị lý thuyết")

* P-P plot
pnorm income, ///
    title("P-P Plot cho thu nhập")

* Kiểm định Jarque-Bera
jb income

* Ví dụ biến đổi dữ liệu nếu không có phân phối chuẩn
* Biến đổi logarithm
generate income_log = ln(income)
label variable income_log "Log của thu nhập"

* Biến đổi căn bậc hai
generate income_sqrt = sqrt(income)
label variable income_sqrt "Căn bậc hai của thu nhập"

* Biến đổi nghịch đảo
generate income_reciprocal = 1/income
label variable income_reciprocal "Nghịch đảo của thu nhập"

* Biến đổi Box-Cox
gladder income

* Kiểm tra lại phân phối chuẩn sau biến đổi
foreach var of varlist income_log income_sqrt income_reciprocal {
    display "Kiểm tra phân phối chuẩn cho `var'"
    swilk `var'
    sktest `var'
    histogram `var', normal title("Histogram của `var'")
}

* So sánh các phương pháp biến đổi
graph combine ///
    (histogram income, normal title("Gốc")) ///
    (histogram income_log, normal title("Log")) ///
    (histogram income_sqrt, normal title("Sqrt")) ///
    (histogram income_reciprocal, normal title("Reciprocal")), ///
    title("So sánh các phương pháp biến đổi")

# Kiểm tra giả định phân phối chuẩn trong R
# Cài đặt và tải các package cần thiết
if (!require(car)) install.packages("car")
if (!require(nortest)) install.packages("nortest")
if (!require(moments)) install.packages("moments")
if (!require(ggplot2)) install.packages("ggplot2")
if (!require(gridExtra)) install.packages("gridExtra")

library(car)
library(nortest)
library(moments)
library(ggplot2)
library(gridExtra)

# Giả sử có dữ liệu thu nhập
# Tạo dữ liệu mẫu cho ví dụ
set.seed(123)
income <- c(rnorm(150, 50000, 15000), rexp(50, 1/30000))

# Phương pháp 1: Tạo histogram với đường cong chuẩn
hist(income, 
     probability = TRUE, 
     main = "Histogram của thu nhập với đường cong chuẩn",
     xlab = "Thu nhập", 
     ylab = "Mật độ",
     col = "lightblue",
     border = "black")

# Thêm đường cong chuẩn
curve(dnorm(x, mean = mean(income), sd = sd(income)), 
      add = TRUE, 
      col = "red", 
      lwd = 2)

# Sử dụng ggplot2 để tạo biểu đồ đẹp hơn
df <- data.frame(income = income)
p1 <- ggplot(df, aes(x = income)) +
  geom_histogram(aes(y = ..density..), bins = 30, alpha = 0.7, fill = "lightblue") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(income), sd = sd(income)),
                color = "red", size = 1) +
  labs(title = "Histogram với đường cong chuẩn", x = "Thu nhập", y = "Mật độ")

# Phương pháp 2: Tính toán Skewness và Kurtosis
cat("Thống kê mô tả:\n")
cat("Trung bình:", mean(income), "\n")
cat("Độ lệch chuẩn:", sd(income), "\n")
cat("Skewness:", skewness(income), "\n")
cat("Kurtosis:", kurtosis(income), "\n")

# Phương pháp 3: Các kiểm định phân phối chuẩn
cat("\nKiểm định Shapiro-Wilk:\n")
shapiro_test <- shapiro.test(income)
print(shapiro_test)

cat("\nKiểm định Kolmogorov-Smirnov:\n")
ks_test <- ks.test(income, "pnorm", mean = mean(income), sd = sd(income))
print(ks_test)

cat("\nKiểm định Anderson-Darling:\n")
ad_test <- ad.test(income)
print(ad_test)

cat("\nKiểm định Jarque-Bera:\n")
jb_test <- jarque.test(income)
print(jb_test)

# Tạo Q-Q plot
p2 <- ggplot(df, aes(sample = income)) +
  stat_qq() +
  stat_qq_line(color = "red") +
  labs(title = "Q-Q Plot", x = "Phân vị lý thuyết", y = "Phân vị mẫu")

# P-P plot
p_values <- pnorm(sort(income), mean = mean(income), sd = sd(income))
theoretical_p <- (1:length(income)) / length(income)
pp_data <- data.frame(theoretical = theoretical_p, observed = p_values)

p3 <- ggplot(pp_data, aes(x = theoretical, y = observed)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red") +
  labs(title = "P-P Plot", x = "Xác suất lý thuyết", y = "Xác suất quan sát")

# Kết hợp các biểu đồ
grid.arrange(p1, p2, p3, ncol = 2)

# Ví dụ biến đổi dữ liệu nếu không có phân phối chuẩn
# Biến đổi logarithm (chỉ với giá trị dương)
income_positive <- income[income > 0]
income_log <- log(income_positive)

# Biến đổi căn bậc hai
income_sqrt <- sqrt(income_positive)

# Biến đổi nghịch đảo
income_reciprocal <- 1/income_positive

# Biến đổi Box-Cox
library(car)
bc_result <- powerTransform(income_positive)
optimal_lambda <- bc_result$lambda
income_boxcox <- bcPower(income_positive, optimal_lambda)

# Hàm kiểm tra phân phối chuẩn
check_normality <- function(data, name) {
  cat("\n", "="*50, "\n")
  cat("Kiểm tra phân phối chuẩn cho:", name, "\n")
  cat("="*50, "\n")
  
  # Thống kê mô tả
  cat("Skewness:", skewness(data), "\n")
  cat("Kurtosis:", kurtosis(data), "\n")
  
  # Kiểm định
  if (length(data) <= 5000) {
    shapiro_result <- shapiro.test(data)
    cat("Shapiro-Wilk p-value:", shapiro_result$p.value, "\n")
  }
  
  ad_result <- ad.test(data)
  cat("Anderson-Darling p-value:", ad_result$p.value, "\n")
  
  jb_result <- jarque.test(data)
  cat("Jarque-Bera p-value:", jb_result$p.value, "\n")
  
  # Tạo biểu đồ
  par(mfrow = c(1, 2))
  hist(data, probability = TRUE, main = paste("Histogram -", name), 
       xlab = name, col = "lightblue")
  curve(dnorm(x, mean = mean(data), sd = sd(data)), add = TRUE, col = "red", lwd = 2)
  
  qqnorm(data, main = paste("Q-Q Plot -", name))
  qqline(data, col = "red")
  
  par(mfrow = c(1, 1))
}

# Kiểm tra tất cả các biến đổi
check_normality(income, "Thu nhập gốc")
check_normality(income_log, "Log thu nhập")
check_normality(income_sqrt, "Căn bậc hai thu nhập")
check_normality(income_reciprocal, "Nghịch đảo thu nhập")
check_normality(income_boxcox, "Box-Cox thu nhập")

# Tạo bảng tóm tắt kết quả
transformations <- c("Gốc", "Log", "Sqrt", "Reciprocal", "Box-Cox")
datasets <- list(income, income_log, income_sqrt, income_reciprocal, income_boxcox)

results_summary <- data.frame(
  Transformation = transformations,
  Skewness = sapply(datasets, skewness),
  Kurtosis = sapply(datasets, kurtosis),
  Shapiro_p = sapply(datasets, function(x) {
    if (length(x) <= 5000) shapiro.test(x)$p.value else NA
  }),
  AD_p = sapply(datasets, function(x) ad.test(x)$p.value),
  JB_p = sapply(datasets, function(x) jarque.test(x)$p.value)
)

print(results_summary)

# Kiểm tra giả định phân phối chuẩn trong Python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import jarque_bera, normaltest, anderson
import warnings
warnings.filterwarnings('ignore')

# Thiết lập matplotlib để hiển thị tiếng Việt
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['figure.figsize'] = (12, 8)

# Tạo dữ liệu mẫu (kết hợp phân phối chuẩn và không chuẩn)
np.random.seed(123)
income = np.concatenate([
    np.random.normal(50000, 15000, 150),  # Phân phối chuẩn
    np.random.exponential(30000, 50)      # Phân phối không chuẩn
])

# Tạo DataFrame
df = pd.DataFrame({'income': income})

print("KIỂM TRA GIẢ ĐỊNH PHÂN PHỐI CHUẨN")
print("="*60)

# Phương pháp 1: Tạo histogram với đường cong chuẩn
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# Histogram với đường cong chuẩn
axes[0, 0].hist(income, bins=30, density=True, alpha=0.7, color='lightblue', edgecolor='black')
mu, sigma = np.mean(income), np.std(income)
x = np.linspace(np.min(income), np.max(income), 100)
axes[0, 0].plot(x, stats.norm.pdf(x, mu, sigma), 'r-', linewidth=2, label='Đường cong chuẩn')
axes[0, 0].set_title('Histogram với đường cong chuẩn')
axes[0, 0].set_xlabel('Thu nhập')
axes[0, 0].set_ylabel('Mật độ')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Phương pháp 2: Tính toán Skewness và Kurtosis
skewness_value = stats.skew(income)
kurtosis_value = stats.kurtosis(income)

print(f"Thống kê mô tả:")
print(f"Trung bình: {np.mean(income):.2f}")
print(f"Độ lệch chuẩn: {np.std(income):.2f}")
print(f"Skewness: {skewness_value:.4f}")
print(f"Kurtosis: {kurtosis_value:.4f}")

# Phương pháp 3: Các kiểm định phân phối chuẩn
print(f"\nCác kiểm định phân phối chuẩn:")
print("-" * 40)

# Kiểm định Shapiro-Wilk (chỉ cho mẫu nhỏ)
if len(income) <= 5000: shapiro_stat, shapiro_p = stats.shapiro(income) print(f"Shapiro-Wilk:") print(f" Statistic: {shapiro_stat:.4f}") print(f" p-value: {shapiro_p:.4f}") print(f" Kết luận: {'Phân phối chuẩn' if shapiro_p > 0.05 else 'Không phân phối chuẩn'}")

# Kiểm định Kolmogorov-Smirnov
ks_stat, ks_p = stats.kstest(income, lambda x: stats.norm.cdf(x, mu, sigma))
print(f"\nKolmogorov-Smirnov:")
print(f"  Statistic: {ks_stat:.4f}")
print(f"  p-value: {ks_p:.4f}")
print(f"  Kết luận: {'Phân phối chuẩn' if ks_p > 0.05 else 'Không phân phối chuẩn'}")

# Kiểm định Anderson-Darling
ad_stat, ad_critical, ad_significance = anderson(income, dist='norm')
print(f"\nAnderson-Darling:")
print(f"  Statistic: {ad_stat:.4f}")
for i, sl in enumerate(ad_significance):
    print(f"  Significance level {sl}%: {ad_critical[i]:.4f}")

# Kiểm định Jarque-Bera
jb_stat, jb_p = jarque_bera(income)
print(f"\nJarque-Bera:")
print(f"  Statistic: {jb_stat:.4f}")
print(f"  p-value: {jb_p:.4f}")
print(f"  Kết luận: {'Phân phối chuẩn' if jb_p > 0.05 else 'Không phân phối chuẩn'}")

# Kiểm định D'Agostino-Pearson
dp_stat, dp_p = normaltest(income)
print(f"\nD'Agostino-Pearson:")
print(f"  Statistic: {dp_stat:.4f}")
print(f"  p-value: {dp_p:.4f}")
print(f"  Kết luận: {'Phân phối chuẩn' if dp_p > 0.05 else 'Không phân phối chuẩn'}")

# Tạo Q-Q plot
stats.probplot(income, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot')
axes[0, 1].grid(True, alpha=0.3)

# Tạo P-P plot
sorted_data = np.sort(income)
p_theoretical = np.arange(1, len(sorted_data) + 1) / len(sorted_data)
p_observed = stats.norm.cdf(sorted_data, mu, sigma)
axes[1, 0].plot(p_theoretical, p_observed, 'bo', markersize=3)
axes[1, 0].plot([0, 1], [0, 1], 'r-', linewidth=2)
axes[1, 0].set_title('P-P Plot')
axes[1, 0].set_xlabel('Xác suất lý thuyết')
axes[1, 0].set_ylabel('Xác suất quan sát')
axes[1, 0].grid(True, alpha=0.3)

# Boxplot để phát hiện outliers
axes[1, 1].boxplot(income, vert=True)
axes[1, 1].set_title('Boxplot để phát hiện outliers')
axes[1, 1].set_ylabel('Thu nhập')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Ví dụ biến đổi dữ liệu nếu không có phân phối chuẩn
print(f"\nBIẾN ĐỔI DỮ LIỆU")
print("="*60)

# Chỉ lấy giá trị dương cho các phép biến đổi
income_positive = income[income > 0]

# Các phương pháp biến đổi
transformations = {
    'Gốc': income_positive,
    'Log': np.log(income_positive),
    'Sqrt': np.sqrt(income_positive),
    'Reciprocal': 1/income_positive,
    'Square': income_positive**2
}

# Hàm kiểm tra phân phối chuẩn
def check_normality(data, name):
    print(f"\n{name}:")
    print("-" * 30)
    
    # Thống kê mô tả
    skew = stats.skew(data)
    kurt = stats.kurtosis(data)
    print(f"Skewness: {skew:.4f}")
    print(f"Kurtosis: {kurt:.4f}")
    
    # Kiểm định
    if len(data) <= 5000: shapiro_stat, shapiro_p = stats.shapiro(data) print(f"Shapiro-Wilk p-value: {shapiro_p:.4f}") shapiro_result = "Chuẩn" if shapiro_p > 0.05 else "Không chuẩn"
        print(f"Shapiro-Wilk: {shapiro_result}")
    
    jb_stat, jb_p = jarque_bera(data)
    print(f"Jarque-Bera p-value: {jb_p:.4f}")
    jb_result = "Chuẩn" if jb_p > 0.05 else "Không chuẩn"
    print(f"Jarque-Bera: {jb_result}")
    
    return {
        'name': name,
        'skewness': skew,
        'kurtosis': kurt,
        'shapiro_p': shapiro_p if len(data) <= 5000 else np.nan,
        'jb_p': jb_p
    }

# Kiểm tra tất cả các biến đổi
results = []
for name, data in transformations.items():
    result = check_normality(data, name)
    results.append(result)

# Tạo bảng tóm tắt kết quả
results_df = pd.DataFrame(results)
print(f"\nBẢNG TÓM TẮT KẾT QUẢ:")
print("="*80)
print(results_df.to_string(index=False, float_format='%.4f'))

# Tạo biểu đồ so sánh các phương pháp biến đổi
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes = axes.flatten()

for i, (name, data) in enumerate(transformations.items()):
    if i < len(axes):
        axes[i].hist(data, bins=30, density=True, alpha=0.7, color='lightblue', edgecolor='black')
        
        # Thêm đường cong chuẩn
        mu_trans, sigma_trans = np.mean(data), np.std(data)
        x_trans = np.linspace(np.min(data), np.max(data), 100)
        axes[i].plot(x_trans, stats.norm.pdf(x_trans, mu_trans, sigma_trans), 
                    'r-', linewidth=2, label='Đường cong chuẩn')
        
        axes[i].set_title(f'{name} (Skew: {stats.skew(data):.3f})')
        axes[i].set_xlabel('Giá trị')
        axes[i].set_ylabel('Mật độ')
        axes[i].legend()
        axes[i].grid(True, alpha=0.3)

# Ẩn subplot thừa
if len(transformations) < len(axes):
    axes[-1].set_visible(False)

plt.tight_layout()
plt.show()

# Tìm phương pháp biến đổi tốt nhất
best_transformation = None
best_score = float('inf')

for result in results:
    if not np.isnan(result['shapiro_p']):
        # Điểm số dựa trên kết hợp của các tiêu chí
        score = abs(result['skewness']) + abs(result['kurtosis']) - result['shapiro_p']
        if score < best_score:
            best_score = score
            best_transformation = result['name']

if best_transformation:
    print(f"\nPHƯƠNG PHÁP BIẾN ĐỔI TỐT NHẤT: {best_transformation}")
    print("="*60)
    
    # Hiển thị so sánh trước và sau biến đổi
    original_data = income_positive
    transformed_data = transformations[best_transformation]
    
    fig, axes = plt.subplots(2, 2, figsize=(15, 10))
    
    # Histogram gốc
    axes[0, 0].hist(original_data, bins=30, density=True, alpha=0.7, color='lightcoral')
    axes[0, 0].set_title('Dữ liệu gốc')
    axes[0, 0].set_xlabel('Thu nhập')
    axes[0, 0].set_ylabel('Mật độ')
    
    # Q-Q plot gốc
    stats.probplot(original_data, dist="norm", plot=axes[0, 1])
    axes[0, 1].set_title('Q-Q Plot - Dữ liệu gốc')
    
    # Histogram sau biến đổi
    axes[1, 0].hist(transformed_data, bins=30, density=True, alpha=0.7, color='lightgreen')
    mu_best, sigma_best = np.mean(transformed_data), np.std(transformed_data)
    x_best = np.linspace(np.min(transformed_data), np.max(transformed_data), 100)
    axes[1, 0].plot(x_best, stats.norm.pdf(x_best, mu_best, sigma_best), 
                   'r-', linewidth=2, label='Đường cong chuẩn')
    axes[1, 0].set_title(f'Sau biến đổi {best_transformation}')
    axes[1, 0].set_xlabel('Giá trị biến đổi')
    axes[1, 0].set_ylabel('Mật độ')
    axes[1, 0].legend()
    
    # Q-Q plot sau biến đổi
    stats.probplot(transformed_data, dist="norm", plot=axes[1, 1])
    axes[1, 1].set_title(f'Q-Q Plot - Sau biến đổi {best_transformation}')
    
    plt.tight_layout()
    plt.show()

print(f"\nKẾT LUẬN:")
print("="*60)
print("1. Luôn kiểm tra phân phối chuẩn trước khi áp dụng kiểm định tham số")
print("2. Kết hợp nhiều phương pháp để có đánh giá toàn diện")
print("3. Xem xét biến đổi dữ liệu nếu cần thiết")
print("4. Với mẫu lớn, tập trung vào phân tích đồ thị hơn là kiểm định thống kê")
print("5. Nhớ rằng phân phối chuẩn là giả định lý tưởng, không phải yêu cầu tuyệt đối")