# 单个随机数
r = rand()              # [0, 1) 均匀分布的 Float64

# 指定类型
ri = rand(Int)          # 随机 Int64
rf = rand(Float32)      # 随机 Float32

# 随机数组
arr1d = rand(5)         # 长度为5的向量
arr2d = rand(3, 4)      # 3×4 矩阵
arr3d = rand(2, 3, 4)   # 2×3×4 三维数组

# 指定范围
r1 = rand(1:10)         # 1到10的随机整数
r2 = rand(1:10, 3)      # 3个1到10的随机整数
r3 = rand(0.0:0.1:1.0)  # 从指定步长的范围抽样

# 从集合中随机选择
colors = ["红", "绿", "蓝", "黄"]
println(rand(colors))      # 随机一个
println(rand(colors, 3))   # 随机3个（有放回）
println(rand(colors, 2; replace=false))  # 无放回

# 标准正态分布 N(0,1)
r = randn()             # 单个值
arr = randn(1000)       # 1000个样本
mat = randn(10, 10)     # 10×10 矩阵

# 指定类型
r32 = randn(Float32, 100)

using Random

# 设置种子确保结果可重复
Random.seed!(42)
a = rand(5)
println(a)

Random.seed!(42)
b = rand(5)
println(a == b)  # true

# 种子可以是任意整数
Random.seed!(12345)
println(rand(3))

# 不指定种子 = 随机种子（不可重复）
Random.seed!()
println(rand(3))  # 每次运行不同

using Pkg
Pkg.add("Distributions")

using Distributions

# 离散分布
bern = Bernoulli(0.5)      # 伯努利分布
binom = Binomial(10, 0.3)  # 二项分布
pois = Poisson(5.0)        # 泊松分布
cat = Categorical([0.2, 0.3, 0.5])  # 分类分布

# 连续分布
norm = Normal(0.0, 1.0)    # 正态分布
unif = Uniform(0.0, 1.0)   # 均匀分布
expo = Exponential(2.0)    # 指数分布
beta = Beta(2.0, 5.0)      # Beta 分布
gamma = Gamma(2.0, 1.0)    # Gamma 分布
cauch = Cauchy(0.0, 1.0)   # 柯西分布
logn = LogNormal(0.0, 1.0) # 对数正态分布

d = Normal(100.0, 15.0)  # 均值100，标准差15

# 统计量
println(mean(d))          # 100.0
println(var(d))           # 225.0
println(std(d))           # 15.0
println(median(d))        # 100.0
println(mode(d))          # 100.0
println(skewness(d))      # 0.0
println(kurtosis(d))      # 0.0

# 概率密度/质量函数
println(pdf(d, 100.0))    # 在x=100处的概率密度
println(pdf(d, 115.0))    # 在x=115处的概率密度

# 累积分布函数
println(cdf(d, 100.0))    # P(X ≤ 100) = 0.5
println(cdf(d, 115.0))    # P(X ≤ 115) ≈ 0.8413

# 分位数函数（逆CDF）
println(quantile(d, 0.5))  # 中位数
println(quantile(d, 0.95)) # 95%分位数

# 对数概率（数值更稳定）
println(logpdf(d, 100.0))
println(logcdf(d, 100.0))

# 采样
samples = rand(d, 1000)     # 1000个样本
println("样本均值: ", mean(samples))
println("样本标准差: ", std(samples))

using Random

# 离散布朗运动
function brownian_motion(T, N; seed=nothing)
    seed !== nothing && Random.seed!(seed)
    dt = T / N
    dW = randn(N) .* sqrt(dt)
    W = cumsum(dW)
    W = prepend!(W, 0.0)  # W(0) = 0
    t = range(0, T; length=N+1)
    return t, W
end

t, W = brownian_motion(1.0, 1000; seed=42)
println("W(T) = $(W[end])")
println("理论方差: T = 1.0, 实际方差: $(var(W))")

function geometric_brownian_motion(S0, μ, σ, T, N; seed=nothing)
    seed !== nothing && Random.seed!(seed)
    dt = T / N
    t = range(0, T; length=N+1)
    
    # GBM: S(t) = S0 * exp((μ - σ²/2)t + σ*W(t))
    W = cumsum(randn(N) .* sqrt(dt))
    W = prepend!(W, 0.0)
    
    S = S0 .* exp.((μ - σ^2/2) .* t .+ σ .* W)
    return t, S
end

t, prices = geometric_brownian_motion(100.0, 0.05, 0.2, 1.0, 252)
println("初始价格: $(prices[1])")
println("最终价格: $(round(prices[end], digits=2))")

function poisson_process(λ, T; seed=nothing)
    seed !== nothing && Random.seed!(seed)
    
    events = Float64[]
    t = 0.0
    while t < T
        t += rand(Exponential(1/λ))
        t < T && push!(events, t)
    end
    return events
end

events = poisson_process(5.0, 10.0; seed=42)
println("事件数: $(length(events))")
println("理论期望: $(5.0 * 10.0)")

# 估算积分 ∫₀¹ sin(x) dx
function monte_carlo_integral(f, a, b, n; seed=nothing)
    seed !== nothing && Random.seed!(seed)
    samples = a .+ (b - a) .* rand(n)
    estimates = f.(samples)
    return (b - a) * mean(estimates), (b - a) * std(estimates) / sqrt(n)
end

# 理论值: 1 - cos(1) ≈ 0.4597
result, error = monte_carlo_integral(sin, 0, 1, 1_000_000; seed=42)
println("估计: $(round(result, digits=6)) ± $(round(error, digits=6))")
println("真值: $(round(1 - cos(1), digits=6))")

function estimate_pi(n; seed=nothing)
    seed !== nothing && Random.seed!(seed)
    inside = 0
    for _ in 1:n
        x, y = rand(2)
        if x^2 + y^2 <= 1
            inside += 1
        end
    end
    return 4 * inside / n
end

for n in [1000, 10000, 100000, 1000000]
    pi_est = estimate_pi(n; seed=42)
    println("n=$n: π ≈ $(round(pi_est, digits=6))")
end

using Random

# 默认生成器
rng = MersenneTwister(1234)
rand(rng, 5)

# 独立的 RNG 实例
rng1 = MersenneTwister(1)
rng2 = MersenneTwister(2)

# 不同种子产生不同序列
println(rand(rng1, 3))
println(rand(rng2, 3))

# 但相同种子产生相同序列
rng3 = MersenneTwister(42)
rng4 = MersenneTwister(42)
println(rand(rng3, 3) == rand(rng4, 3))  # true

# Xoshiro256++ (Julia 1.7+ 默认)
rng = Random.Xoshiro(1234)
rand(rng, 5)

# 指定默认生成器类型
Random.default_rng()  # 获取默认 RNG

using Random

# 完整的可重复实验框架
function run_experiment(; seed=42, n_samples=10000)
    Random.seed!(seed)
    
    # 生成数据
    X = randn(n_samples, 10)
    β = [1.0, -2.0, 3.0, 0, 0, 0.5, -0.5, 0, 0, 0]
    y = X * β .+ 0.1 .* randn(n_samples)
    
    # 线性回归
    β_hat = X \ y
    
    # 评估
    mse = mean((β - β_hat).^2)
    return β_hat, mse
end

# 多次运行相同结果
for i in 1:3
    β_hat, mse = run_experiment(seed=42)
    println("MSE: $(round(mse, digits=6))")
end

# 不同种子得到不同但统计一致的结果
mses = [run_experiment(seed=s)[2] for s in 1:100]
println("平均 MSE: $(round(mean(mses), digits=6))")
println("MSE 标准差: $(round(std(mses), digits=6))")

using Statistics

data = randn(1000) .* 10 .+ 50  # N(50, 100)

println("均值: ", round(mean(data), digits=2))    # ≈50
println("中位数: ", round(median(data), digits=2)) # ≈50
println("标准差: ", round(std(data), digits=2))    # ≈10
println("方差: ", round(var(data), digits=2))      # ≈100
println("最小: ", round(minimum(data), digits=2))
println("最大: ", round(maximum(data), digits=2))
println("偏度: ", round(skewness(data), digits=4))
println("峰度: ", round(kurtosis(data), digits=4))

# 分位数
println("25%分位: ", round(quantile(data, 0.25), digits=2))
println("75%分位: ", round(quantile(data, 0.75), digits=2))

# using HypothesisTests

# 单样本 t 检验
# data = randn(100) .+ 0.5
# OneSampleTTest(data, 0.0)  # H0: μ = 0

# 双样本 t 检验
# TwoSampleTTest(group1, group2)

# 正态性检验
# ExactOneSampleKSTest(data, Normal(0, 1))

using Statistics

function confidence_interval(data; level=0.95)
    n = length(data)
    μ = mean(data)
    σ = std(data)
    
    # z 分位数（大样本近似）
    z = quantile(Normal(), 1 - (1-level)/2)
    margin = z * σ / sqrt(n)
    
    return (μ - margin, μ + margin)
end

data = randn(100) .+ 5.0
ci = confidence_interval(data)
println("95% 置信区间: ($(round(ci[1], digits=3)), $(round(ci[2], digits=3)))")

using Random, Statistics

function monte_carlo_call(S0, K, r, σ, T, n_paths; seed=nothing)
    """
    S0: 初始股价
    K:  执行价
    r:  无风险利率
    σ:  波动率
    T:  到期时间
    n_paths: 蒙特卡洛路径数
    """
    seed !== nothing && Random.seed!(seed)
    
    # 生成终值
    Z = randn(n_paths)
    ST = S0 .* exp.((r - σ^2/2) .* T .+ σ .* sqrt(T) .* Z)
    
    # 计算收益
    payoff = max.(ST .- K, 0)
    
    # 折现
    price = exp(-r * T) * mean(payoff)
    se = exp(-r * T) * std(payoff) / sqrt(n_paths)
    
    return price, se
end

# 参数
S0 = 100.0    # 初始价格
K = 105.0     # 执行价
r = 0.05      # 无风险利率
σ = 0.2       # 波动率
T = 1.0       # 1年到期

price, se = monte_carlo_call(S0, K, r, σ, T, 100_000; seed=42)
println("期权价格: $(round(price, digits=4))")
println("标准误差: $(round(se, digits=4))")

# 与 Black-Scholes 解析解对比
using Distributions
d1 = (log(S0/K) + (r + σ^2/2)*T) / (σ*sqrt(T))
d2 = d1 - σ*sqrt(T)
bs_price = S0 * cdf(Normal(), d1) - K * exp(-r*T) * cdf(Normal(), d2)
println("BS 公式价格: $(round(bs_price, digits=4))")

function convergence_analysis(S0, K, r, σ, T)
    ns = [100, 1000, 10000, 100000, 1000000]
    
    println("样本数    价格估计     标准误差     95%置信区间")
    println("-" ^ 60)
    
    for n in ns
        price, se = monte_carlo_call(S0, K, r, σ, T, n; seed=42)
        ci_lower = price - 1.96 * se
        ci_upper = price + 1.96 * se
        println("$(lpad(n, 8))  $(round(price, digits=4))" *
                "       $(round(se, digits=4))" *
                "       [$(round(ci_lower, digits=4)), $(round(ci_upper, digits=4))]")
    end
end

convergence_analysis(S0, K, r, σ, T)

function compute_delta(S0, K, r, σ, T, n_paths; ΔS=1.0, seed=42)
    # 中心差分法计算 Delta
    price_up, _ = monte_carlo_call(S0 + ΔS, K, r, σ, T, n_paths; seed=seed)
    price_dn, _ = monte_carlo_call(S0 - ΔS, K, r, σ, T, n_paths; seed=seed)
    return (price_up - price_dn) / (2 * ΔS)
end

delta = compute_delta(S0, K, r, σ, T, 100_000)
println("Delta: $(round(delta, digits=4))")

function portfolio_var(returns, weights; confidence=0.95, n_sim=10000)
    # 历史模拟法
    portfolio_returns = returns * weights
    
    # 排序取分位数
    sorted = sort(portfolio_returns)
    var_index = Int(floor((1 - confidence) * n_sim)) + 1
    var = -sorted[var_index]
    
    # Expected Shortfall (CVaR)
    cvar = -mean(sorted[1:var_index])
    
    return var, cvar
end

# 模拟1000天的5个资产收益
Random.seed!(42)
n_assets = 5
n_days = 1000
returns = randn(n_days, n_assets) .* 0.02 .+ 0.0005  # 日收益
weights = [0.3, 0.25, 0.2, 0.15, 0.1]

var, cvar = portfolio_var(returns, weights)
println("VaR (95%): $(round(var*100, digits=2))%")
println("CVaR (95%): $(round(cvar*100, digits=2))%")