using StatsBase, Statistics

data = [23, 45, 12, 67, 34, 89, 56, 78, 43, 21, 65, 38]

# 基本统计量
mean(data)              # 均值
median(data)            # 中位数
std(data)               # 标准差
var(data)               # 方差
minimum(data)           # 最小值
maximum(data)           # 最大值
extrema(data)           # (最小, 最大)

# 分位数
quantile(data, [0.25, 0.5, 0.75])   # 四分位数
iqr(data)                            # 四分位距

# 高阶统计
skewness(data)          # 偏度（对称性）
kurtosis(data)          # 峰度（尾部厚度）

# 汇总统计
summarystats(data)      # 一次性输出所有基本统计量

# 频率统计
categorical_data = ["A", "B", "A", "C", "B", "A", "A", "C"]
counts = countmap(categorical_data)
proportions(categorical_data)   # 比例

using StatsBase

data = [10, 20, 30, 40, 50]
weights = [1.0, 2.0, 3.0, 2.0, 1.0]

# 加权均值
wmean = mean(data, weights)

# 加权中位数
wmedian = median(data, AnalyticWeights(weights))

# 加权标准差
wstd = std(data, weights)

# 加权方差
wvar = var(data, weights)

# 不同权重类型
aw = AnalyticWeights([1, 2, 3, 4, 5])     # 解释为重复次数的连续推广
fw = FrequencyWeights([1, 2, 3, 4, 5])    # 解释为频率计数
pw = ProbabilityWeights([0.1, 0.2, 0.3, 0.2, 0.2])  # 解释为概率权重

using HypothesisTests

# 单样本 t 检验
data = randn(50) .+ 0.5    # 均值为 0.5 的样本
t_test = OneSampleTTest(data, 0.0)    # 检验均值是否为 0
pvalue(t_test)                        # p 值
confint(t_test, 0.95)                # 95% 置信区间
show(t_test)

# 双样本 t 检验
group_a = randn(30) .+ 1.0
group_b = randn(35) .+ 1.5
t_test2 = UnequalVarianceTTest(group_a, group_b)   # Welch t 检验
pvalue(t_test2)

# 配对 t 检验
before = randn(20) .+ 5.0
after = before .+ randn(20) .* 0.5 .+ 0.3
paired_test = OneSampleTTest(after - before)
pvalue(paired_test)

# 卡方检验
observed = [50, 30, 20]
expected = [1/3, 1/3, 1/3]
chi2_test = ChisqTest(observed, expected)
pvalue(chi2_test)

# Fisher 精确检验（2×2 列联表）
table = [3 1; 1 3]
fisher_test = FisherExactTest(table...)
pvalue(fisher_test)

# Mann-Whitney U 检验（非参数）
u_test = MannWhitneyUTest(randn(20), randn(25) .+ 0.5)
pvalue(u_test)

using GLM, DataFrames

# 线性回归
n = 100
X = randn(n, 3)
β_true = [2.0, -1.0, 0.5, 3.0]   # 截距 + 3 个系数
y = β_true[1] .+ X * β_true[2:4] .+ randn(n) * 0.5

df = DataFrame(X, [:x1, :x2, :x3])
df.y = y

# 拟合线性模型
model = lm(@formula(y ~ x1 + x2 + x3), df)

# 查看结果
coef(model)           # 回归系数
stderror(model)       # 标准误
confint(model)        # 置信区间
r2(model)             # R²
adjr2(model)          # 调整 R²
deviance(model)       # 残差平方和
dof(model)            # 自由度
predict(model)        # 拟合值
residuals(model)      # 残差

# 使用 DataFrame 构建模型
formula_str = @formula(y ~ 1 + x1 + x2 + x3)
glm_model = lm(formula_str, df)

using GLM, DataFrames

# 二分类数据
n = 200
X = randn(n, 2)
prob = 1 ./ (1 .+ exp.(-(X[:,1] .+ 2X[:,2] .- 1)))    # sigmoid
y = Float64.(rand(n) .< prob)

df = DataFrame(x1=X[:,1], x2=X[:,2], y=y)

# 逻辑回归
logit_model = glm(@formula(y ~ x1 + x2), df, Binomial(), LogitLink())

# 预测概率
predicted_probs = predict(logit_model)
predicted_class = predicted_probs .> 0.5

# 准确率
accuracy = mean(predicted_class .== y)
println("准确率: $(round(accuracy*100, digits=1))%")

using Turing, Distributions

# 简单贝叶斯线性回归
@model function bayesian_linear(x, y)
    # 先验分布
    α ~ Normal(0, 10)       # 截距先验
    β ~ Normal(0, 10)       # 斜率先验
    σ ~ Exponential(1)      # 噪声标准差先验

    # 似然
    for i in eachindex(y)
        y[i] ~ Normal(α + β * x[i], σ)
    end
end

# 生成数据
x_data = collect(1.0:10.0)
y_data = 2.5 .* x_data .+ 1.3 .+ randn(10) .* 0.5

# MCMC 采样
model = bayesian_linear(x_data, y_data)
chain = sample(model, NUTS(1000, 0.65), 2000)

# 查看结果
display(chain)

using Turing

@model function my_model(data)
    μ ~ Normal(0, 10)
    σ ~ Exponential(1)
    for i in eachindex(data)
        data[i] ~ Normal(μ, σ)
    end
end

data = randn(100) .+ 3.0
model = my_model(data)

# NUTS（No-U-Turn Sampler）— 推荐
chain_nuts = sample(model, NUTS(1000, 0.65), 2000)

# HMC（Hamiltonian Monte Carlo）
chain_hmc = sample(model, HMC(0.05, 10), 2000)

# Metropolis-Hastings（经典方法）
chain_mh = sample(model, MH(), 10000)

# SMC（Sequential Monte Carlo）
chain_smc = sample(model, SMC(), 2000)

# 多链并行
chain_multi = sample(model, NUTS(1000, 0.65), MCMCThreads(), 2000, 4)   # 4 条链

using Turing, MCMCChains, StatsPlots

# 采样后分析链
chain = sample(model, NUTS(1000, 0.65), 2000)

# 基本统计
mean(chain)             # 后验均值
median(chain)           # 后验中位数
std(chain)              # 后验标准差
quantile(chain)         # 后验分位数

# HPDI（最高后验密度区间）
hpdi(chain; alpha=0.05)  # 95% HPDI

# 可视化
plot(chain)                 # 迹图 + 密度图
corner(chain)               # 相关矩阵图
autocorplot(chain)          # 自相关图

# 诊断
rhat(chain)                 # R-hat 统计量（应接近 1.0）
ess(chain)                  # 有效样本量

using Turing, Distributions

# 分层贝叶斯模型示例：多学校数据
# 每个学校有自己的效果，但所有学校共享一个超参数
@model function hierarchical_schools(J, y, σ)
    # 超先验
    μ ~ Normal(0, 10)          # 总体均值
    τ ~ Exponential(1)         # 总体标准差

    # 学校级别参数
    θ = Vector{Real}(undef, J)
    for j in 1:J
        θ[j] ~ Normal(μ, τ)   # 学校效果先验
        y[j] ~ Normal(θ[j], σ[j])   # 似然
    end
end

# 8 所学校的数据（经典 8-schools 问题）
J = 8
y = [28, 8, -3, 7, -1, 1, 18, 12]      # 效果估计
σ = [15, 10, 16, 11, 9, 11, 10, 18]    # 标准误

model = hierarchical_schools(J, y, σ)
chain = sample(model, NUTS(2000, 0.8), 4000)

# 查看收缩效应（shrinkage）
display(chain)

using Turing, Statistics

# AIC（Akaike Information Criterion）
# AIC = 2k - 2ln(L̂)
function aic(model, chain, data)
    k = length(chain)       # 参数个数
    # 使用后验均值计算对数似然
    params = mean(chain)
    log_lik = loglikelihood(model, params, data)
    return 2 * k - 2 * log_lik
end

# BIC（Bayesian Information Criterion）
# BIC = k·ln(n) - 2ln(L̂)
function bic(model, chain, data)
    k = length(chain)
    n = length(data)
    params = mean(chain)
    log_lik = loglikelihood(model, params, data)
    return k * log(n) - 2 * log_lik
end

# WAIC（Widely Applicable Information Criterion）
# 更适合贝叶斯模型
using StatsFuns
function waic(chain_loglik)
    # chain_loglik: n_samples × n_data 矩阵
    lpd = logsumexp(chain_loglik, dims=1) .- log(size(chain_loglik, 1))
    p_waic = var(chain_loglik, dims=1)
    return -2 * (sum(lpd) - sum(p_waic))
end

using HypothesisTests, Distributions, Turing

# 传统方法：频率派 A/B 测试
n_a, conversions_a = 1000, 120     # A 组
n_b, conversions_b = 1000, 150     # B 组

# 比例检验
prop_test = BinomialTest(conversions_a + conversions_b, n_a + n_b,
    conversions_a / n_a)

# 两比例 z 检验
p_a = conversions_a / n_a
p_b = conversions_b / n_b
p_pool = (conversions_a + conversions_b) / (n_a + n_b)
se = sqrt(p_pool * (1 - p_pool) * (1/n_a + 1/n_b))
z = (p_b - p_a) / se
p_value = 2 * (1 - cdf(Normal(), abs(z)))

println("A 转化率: $(round(p_a*100, digits=1))%")
println("B 转化率: $(round(p_b*100, digits=1))%")
println("p 值: $(round(p_value, digits=4))")

# 贝叶斯方法
@model function ab_test_model(n_a, c_a, n_b, c_b)
    θ_a ~ Beta(1, 1)    # 均匀先验
    θ_b ~ Beta(1, 1)
    c_a ~ Binomial(n_a, θ_a)
    c_b ~ Binomial(n_b, θ_b)
    # 计算差异
    Δθ = θ_b - θ_a
end

model = ab_test_model(n_a, conversions_a, n_b, conversions_b)
chain = sample(model, NUTS(1000, 0.65), 4000)

# P(B > A)
Δ_samples = chain[:θ_b] .- chain[:θ_a]
prob_b_better = mean(Δ_samples .> 0)
println("B 优于 A 的概率: $(round(prob_b_better*100, digits=1))%")

using Turing, Plots

# 指数衰减模型参数估计
@model function decay_model(t, y)
    # 先验
    A ~ LogNormal(log(10), 0.5)     # 初始振幅
    λ ~ LogNormal(log(0.5), 0.3)    # 衰减常数
    σ ~ Exponential(1)               # 噪声

    # 似然
    for i in eachindex(y)
        μ = A * exp(-λ * t[i])
        y[i] ~ Normal(μ, σ)
    end
end

# 生成数据
t_data = collect(0.0:0.5:10.0)
y_data = 10.0 .* exp.(-0.3 .* t_data) .+ randn(length(t_data)) .* 0.5

# 采样
model = decay_model(t_data, y_data)
chain = sample(model, NUTS(2000, 0.8), 4000)

# 后验预测
t_pred = 0:0.1:12
predictions = [chain[:A][i] .* exp.(-chain[:λ][i] .* t_pred) for i in 1:size(chain, 1)]
pred_mat = reduce(hcat, predictions)'
median_pred = median(pred_mat, dims=1)[:]
lower = [quantile(pred_mat[:,j], 0.025) for j in 1:length(t_pred)]
upper = [quantile(pred_mat[:,j], 0.975) for j in 1:length(t_pred)]

using StatsPlots
plot(t_pred, median_pred,    ribbon=(median_pred .- lower, upper .- median_pred),
    label="后验预测", fillalpha=0.3)
scatter!(t_data, y_data, label="数据", color=:red)
xlabel!("时间"); ylabel!("振幅"); title!("贝叶斯参数估计")