简单随机抽样：统计学的基石与实践

E(x̄) = μ

lim(n→∞) P(|x̄ - μ| > ε) = 0

Var(x̄) = σ²/n × (N-n)/(N-1)

# 使用操作系统熵源
import os
import struct
 
def get_trng_bytes(num_bytes):
    """获取真随机字节"""
    return os.urandom(num_bytes)
 
# 生成真随机整数
trng_int = struct.unpack('I', get_trng_bytes(4))[0]
print(f"真随机整数: {trng_int}")

// Python的Mersenne Twister实现
class MersenneTwister {
  constructor(seed = Date.now()) {
    this.mt = new Uint32Array(624);
    this.index = 624;
    this.initialize(seed);
  }
 
  initialize(seed) {
    this.mt[0] = seed >>> 0;
    for (let i = 1; i < 624; i++) {
      this.mt[i] = (1812433253 * (this.mt[i-1] ^ (this.mt[i-1] >>> 30)) + i) >>> 0;
    }
  }
 
  extractNumber() {
    if (this.index >= 624) {
      this.generateNumbers();
    }
 
    let y = this.mt[this.index++];
    y ^= y >>> 11;
    y ^= (y << 7) & 0x9d2c5680;
    y ^= (y << 15) & 0xefc60000;
    y ^= y >>> 18;
 
    return y >>> 0;
  }
 
  generateNumbers() {
    for (let i = 0; i < 624; i++) {
      let y = this.mt[i] & 0x80000000 + this.mt[(i + 1) % 624] & 0x7fffffff;
      this.mt[i] = this.mt[(i + 397) % 624] ^ (y >>> 1);
      if (y % 2) {
        this.mt[i] ^= 0x9908b0df;
      }
    }
    this.index = 0;
  }
 
  // 简单随机抽样实现
  simpleRandomSample(population, sampleSize) {
    const shuffled = [...population];
    // Fisher-Yates洗牌
    for (let i = shuffled.length - 1; i > 0; i--) {
      const j = Math.floor(this.extractNumber() / 0xFFFFFFFF * (i + 1));
      [shuffled[i], shuffled[j]] = [shuffled[j], shuffled[i]];
    }
    return shuffled.slice(0, sampleSize);
  }
}
 
// 使用示例
const mt = new MersenneTwister(12345);
const population = Array.from({ length: 1000 }, (_, i) => `单位${i + 1}`);
const sample = mt.simpleRandomSample(population, 50);
console.log("抽中的样本:", sample);

随机数表片段：
57385 84234 82937 48392 48234 93847 23482 84732 39482 48293
84729 39482 48293 84729 39482 48293 84729 39482 48293 84729
...

import random
import numpy as np
 
def simple_random_sampling(population, sample_size, seed=None):
    """
    简单随机抽样
 
    参数:
    - population: 总体列表
    - sample_size: 样本容量
    - seed: 随机种子（可选）
 
    返回:
    - 样本列表
    """
    if seed is not None:
        random.seed(seed)
 
    return random.sample(population, sample_size)
 
# 使用示例
population = list(range(1, 10001))  # 1-10000的总体
sample = simple_random_sampling(population, 100, seed=42)
print(f"样本容量: {len(sample)}")
print(f"前10个样本单位: {sample[:10]}")

def numpy_simple_sampling(population_size, sample_size, seed=42):
    """使用NumPy进行简单随机抽样"""
    np.random.seed(seed)
 
    # 无放回抽样
    sample_indices = np.random.choice(
        population_size,
        size=sample_size,
        replace=False
    )
 
    return sample_indices + 1  # 转为1-based编号
 
# 使用示例
sample_indices = numpy_simple_sampling(10000, 100)
print("样本索引:", sample_indices)

def fisher_yates_sample(population, sample_size):
    """使用Fisher-Yates算法进行抽样"""
    population_copy = population.copy()
 
    # 洗牌过程
    for i in range(len(population_copy) - 1, 0, -1):
        j = random.randint(0, i)
        population_copy[i], population_copy[j] = population_copy[j], population_copy[i]
 
    return population_copy[:sample_size]
 
# 性能对比
import time
 
population = list(range(1, 100001))  # 10万总体
sample_size = 1000
 
# 方法1：random.sample
start = time.time()
sample1 = simple_random_sampling(population, sample_size)
time1 = time.time() - start
 
# 方法2：Fisher-Yates
start = time.time()
sample2 = fisher_yates_sample(population, sample_size)
time2 = time.time() - start
 
print(f"random.sample耗时: {time1:.4f}秒")
print(f"Fisher-Yates耗时: {time2:.4f}秒")

x̄ ~ N(μ, σ²/n)

SE = σ/√n

(x̄ - μ) / (s/√n) ~ t(n-1)

(n-1)s²/σ² ~ χ²(n-1)

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
 
# 模拟总体数据（全国消费者满意度，1-10分）
np.random.seed(42)
population_satisfaction = np.random.normal(7.2, 1.5, 100000000)
 
print(f"总体均值: {np.mean(population_satisfaction):.2f}")
print(f"总体标准差: {np.std(population_satisfaction):.2f}")
 
# 计算所需样本量
z_score = 1.96  # 95%置信水平的z值
margin_error = 0.015
sigma = 1.5
 
required_sample_size = (z_score * sigma / margin_error) ** 2
print(f"所需样本量: {int(required_sample_size)}")
 
# 进行简单随机抽样
sample_size = 5000
sample_indices = np.random.choice(
    len(population_satisfaction),
    size=sample_size,
    replace=False
)
sample_data = population_satisfaction[sample_indices]
 
print(f"\n样本统计:")
print(f"样本均值: {np.mean(sample_data):.3f}")
print(f"样本标准差: {np.std(sample_data, ddof=1):.3f}")
 
# 计算置信区间
standard_error = np.std(sample_data, ddof=1) / np.sqrt(sample_size)
margin = z_score * standard_error
ci_lower = np.mean(sample_data) - margin
ci_upper = np.mean(sample_data) + margin
 
print(f"95%置信区间: [{ci_lower:.3f}, {ci_upper:.3f}]")
 
# 验证：是否包含真实总体均值
contains_true_mean = ci_lower <= np.mean(population_satisfaction) <= ci_upper
print(f"置信区间是否包含真实均值: {contains_true_mean}")
 
# 绘制分布图
plt.figure(figsize=(12, 6))
 
plt.subplot(1, 2, 1)
plt.hist(population_satisfaction[:10000], bins=50, alpha=0.7, label='总体')
plt.hist(sample_data, bins=50, alpha=0.7, label='样本')
plt.xlabel('满意度评分')
plt.ylabel('频次')
plt.title('总体 vs 样本分布对比')
plt.legend()
 
plt.subplot(1, 2, 2)
# 模拟多次抽样，验证置信区间的有效性
n_simulations = 1000
coverage_count = 0
sample_means = []
 
for i in range(n_simulations):
    sim_sample = np.random.choice(
        population_satisfaction,
        size=sample_size,
        replace=False
    )
    sample_means.append(np.mean(sim_sample))
 
    # 计算置信区间
    sim_se = np.std(sim_sample, ddof=1) / np.sqrt(sample_size)
    sim_ci = [np.mean(sim_sample) - z_score * sim_se,
              np.mean(sim_sample) + z_score * sim_se]
 
    if sim_ci[0] <= np.mean(population_satisfaction) <= sim_ci[1]:
        coverage_count += 1
 
coverage_rate = coverage_count / n_simulations
print(f"1000次模拟中置信区间覆盖率: {coverage_rate:.3f}")
 
plt.hist(sample_means, bins=30, alpha=0.7)
plt.axvline(np.mean(population_satisfaction), color='red',
           linestyle='--', label='真实总体均值')
plt.axvline(np.mean(sample_data), color='green',
           linestyle='--', label='当前样本均值')
plt.xlabel('样本均值')
plt.ylabel('频次')
plt.title(f'样本均值分布 (n={sample_size})')
plt.legend()
 
plt.tight_layout()
plt.show()

import pandas as pd
from datetime import datetime, timedelta
 
class ClinicalTrialRandomizer:
    """临床试验随机分组系统"""
 
    def __init__(self, participants, treatment_ratio=1.0, seed=None):
        self.participants = participants
        self.treatment_ratio = treatment_ratio
        self.assignments = {}
        if seed:
            np.random.seed(seed)
 
    def stratified_randomization(self, strata_columns):
        """
        分层随机分组
        """
        # 按层分组
        strata = {}
        for idx, participant in self.participants.iterrows():
            # 创建层标识
            stratum_key = tuple(participant[col] for col in strata_columns)
            if stratum_key not in strata:
                strata[stratum_key] = []
            strata[stratum_key].append(idx)
 
        # 每层内随机分配
        for stratum_key, indices in strata.items():
            # 随机打乱
            np.random.shuffle(indices)
 
            # 计算分组数量
            total_in_stratum = len(indices)
            treatment_count = int(total_in_stratum * self.treatment_ratio / (1 + self.treatment_ratio))
 
            # 分配治疗组和对照组
            for i, idx in enumerate(indices):
                if i < treatment_count:
                    self.assignments[idx] = 'treatment'
                else:
                    self.assignments[idx] = 'control'
 
    def simple_randomization(self):
        """
        简单随机分组
        """
        participant_indices = list(self.participants.index)
        np.random.shuffle(participant_indices)
 
        n_treatment = int(len(participant_indices) * self.treatment_ratio / (1 + self.treatment_ratio))
 
        for i, idx in enumerate(participant_indices):
            if i < n_treatment:
                self.assignments[idx] = 'treatment'
            else:
                self.assignments[idx] = 'control'
 
    def get_balance_statistics(self):
        """获取分组平衡性统计"""
        assignment_df = pd.DataFrame.from_dict(self.assignments, orient='index', columns=['group'])
 
        balance_stats = {}
        for col in self.participants.columns:
            if self.participants[col].dtype in ['object', 'category']:
                # 分类变量
                crosstab = pd.crosstab(assignment_df['group'], self.participants[col])
                balance_stats[col] = crosstab
 
                # 卡方检验
                chi2, p_value, dof, expected = stats.chi2_contingency(crosstab)
                print(f"{col} - 卡方检验 p值: {p_value:.4f}")
 
        return balance_stats
 
# 模拟临床试验数据
np.random.seed(42)
n_participants = 1000
 
participants_data = {
    'participant_id': [f'P{i:04d}' for i in range(1, n_participants + 1)],
    'age': np.random.normal(65, 12, n_participants),
    'gender': np.random.choice(['M', 'F'], n_participants),
    'disease_stage': np.random.choice(['I', 'II', 'III'], n_participants, p=[0.3, 0.5, 0.2]),
    'baseline_score': np.random.normal(50, 15, n_participants)
}
 
participants_df = pd.DataFrame(participants_data)
 
# 简单随机分组
randomizer = ClinicalTrialRandomizer(participants_df, seed=42)
randomizer.simple_randomization()
 
# 分析分组结果
assignments_df = pd.DataFrame.from_dict(randomizer.assignments, orient='index', columns=['group'])
participants_df['group'] = assignments_df['group']
 
print("=== 分组统计 ===")
print(participants_df['group'].value_counts())
 
print("\n=== 年龄分布平衡性 ===")
age_by_group = participants_df.groupby('group')['age'].describe()
print(age_by_group)
 
# t检验验证年龄分布
treatment_ages = participants_df[participants_df['group'] == 'treatment']['age']
control_ages = participants_df[participants_df['group'] == 'control']['age']
 
t_stat, p_value = stats.ttest_ind(treatment_ages, control_ages)
print(f"年龄t检验 p值: {p_value:.4f}")
 
print("\n=== 性别分布平衡性 ===")
gender_crosstab = pd.crosstab(participants_df['group'], participants_df['gender'])
print(gender_crosstab)
 
# 卡方检验性别分布
chi2, p_value, dof, expected = stats.chi2_contingency(gender_crosstab)
print(f"性别卡方检验 p值: {p_value:.4f}")
 
print("\n=== 疾病分期分布平衡性 ===")
stage_crosstab = pd.crosstab(participants_df['group'], participants_df['disease_stage'])
print(stage_crosstab)
 
chi2, p_value, dof, expected = stats.chi2_contingency(stage_crosstab)
print(f"疾病分期卡方检验 p值: {p_value:.4f}")
 
# 绘制平衡性可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
 
# 年龄分布
axes[0, 0].hist(treatment_ages, alpha=0.7, label='治疗组', bins=20)
axes[0, 0].hist(control_ages, alpha=0.7, label='对照组', bins=20)
axes[0, 0].set_xlabel('年龄')
axes[0, 0].set_ylabel('频次')
axes[0, 0].set_title('年龄分布对比')
axes[0, 0].legend()
 
# 性别分布
gender_crosstab.plot(kind='bar', ax=axes[0, 1])
axes[0, 1].set_title('性别分布对比')
axes[0, 1].set_xlabel('组别')
axes[0, 1].set_ylabel('人数')
axes[0, 1].legend(title='性别')
 
# 疾病分期分布
stage_crosstab.plot(kind='bar', ax=axes[1, 0])
axes[1, 0].set_title('疾病分期分布对比')
axes[1, 0].set_xlabel('组别')
axes[1, 0].set_ylabel('人数')
axes[1, 0].legend(title='分期')
 
# 基线评分分布
treatment_scores = participants_df[participants_df['group'] == 'treatment']['baseline_score']
control_scores = participants_df[participants_df['group'] == 'control']['baseline_score']
 
axes[1, 1].hist(treatment_scores, alpha=0.7, label='治疗组', bins=20)
axes[1, 1].hist(control_scores, alpha=0.7, label='对照组', bins=20)
axes[1, 1].set_xlabel('基线评分')
axes[1, 1].set_ylabel('频次')
axes[1, 1].set_title('基线评分分布对比')
axes[1, 1].legend()
 
plt.tight_layout()
plt.show()

import pandas as pd
import numpy as np
from scipy import stats
 
class QualityControlSampler:
    """质量控制抽样系统"""
 
    def __init__(self, batch_size, defect_rate=0.05, seed=None):
        self.batch_size = batch_size
        self.defect_rate = defect_rate
        if seed:
            np.random.seed(seed)
 
        # 生成模拟批次数据
        self.batch_data = self.generate_batch_data()
 
    def generate_batch_data(self):
        """生成批次数据"""
        # 生成产品ID
        product_ids = [f'P{i:06d}' for i in range(1, self.batch_size + 1)]
 
        # 生成质量指标（正态分布）
        length = np.random.normal(100.0, 0.5, self.batch_size)
        weight = np.random.normal(50.0, 2.0, self.batch_size)
 
        # 生成缺陷标识（基于缺陷率）
        is_defect = np.random.choice([0, 1], self.batch_size,
                                    p=[1-self.defect_rate, self.defect_rate])
 
        # 添加一些超规产品（非缺陷但质量接近边界）
        near_limit = np.random.choice([0, 1], self.batch_size, p=[0.95, 0.05])
        length[near_limit == 1] += np.random.normal(0, 0.3,
                                                   np.sum(near_limit))
 
        return pd.DataFrame({
            'product_id': product_ids,
            'length': length,
            'weight': weight,
            'is_defect': is_defect
        })
 
    def simple_random_sample(self, sample_size):
        """简单随机抽样"""
        sample_indices = np.random.choice(
            self.batch_size,
            size=sample_size,
            replace=False
        )
        return self.batch_data.iloc[sample_indices]
 
    def estimate_quality_parameters(self, sample_data):
        """估计质量参数"""
        # 估计缺陷率
        defect_rate_est = sample_data['is_defect'].mean()
 
        # 估计长度均值和标准差
        length_mean = sample_data['length'].mean()
        length_std = sample_data['length'].std(ddof=1)
 
        # 估计重量均值和标准差
        weight_mean = sample_data['weight'].mean()
        weight_std = sample_data['weight'].std(ddof=1)
 
        return {
            'defect_rate': defect_rate_est,
            'length_mean': length_mean,
            'length_std': length_std,
            'weight_mean': weight_mean,
            'weight_std': weight_std
        }
 
    def calculate_confidence_intervals(self, sample_data, confidence=0.95):
        """计算置信区间"""
        alpha = 1 - confidence
        z_score = stats.norm.ppf(1 - alpha/2)
 
        intervals = {}
 
        # 缺陷率的置信区间
        p_hat = sample_data['is_defect'].mean()
        n = len(sample_data)
        se_p = np.sqrt(p_hat * (1 - p_hat) / n)
        intervals['defect_rate'] = [
            p_hat - z_score * se_p,
            p_hat + z_score * se_p
        ]
 
        # 长度均值的置信区间
        length_mean = sample_data['length'].mean()
        length_se = sample_data['length'].std(ddof=1) / np.sqrt(n)
        intervals['length_mean'] = [
            length_mean - z_score * length_se,
            length_mean + z_score * length_se
        ]
 
        # 重量均值的置信区间
        weight_mean = sample_data['weight'].mean()
        weight_se = sample_data['weight'].std(ddof=1) / np.sqrt(n)
        intervals['weight_mean'] = [
            weight_mean - z_score * weight_se,
            weight_mean + z_score * weight_se
        ]
 
        return intervals
 
# 实施质量控制
batch_size = 100000
qc = QualityControlSampler(batch_size, defect_rate=0.03, seed=42)
 
# 进行抽样
sample_size = 2000  # 2%抽检
sample_data = qc.simple_random_sample(sample_size)
 
print("=== 批次质量评估报告 ===")
print(f"批次规模: {qc.batch_data.shape[0]:,}件")
print(f"抽检数量: {sample_size:,}件 ({sample_size/batch_size*100:.1f}%)")
 
# 估计质量参数
quality_params = qc.estimate_quality_parameters(sample_data)
 
print(f"\n估计缺陷率: {quality_params['defect_rate']:.4f} "
      f"(真实值: {qc.defect_rate:.4f})")
 
print(f"\n长度指标:")
print(f"  估计均值: {quality_params['length_mean']:.3f}mm")
print(f"  标准差: {quality_params['length_std']:.3f}mm")
 
print(f"\n重量指标:")
print(f"  估计均值: {quality_params['weight_mean']:.3f}g")
print(f"  标准差: {quality_params['weight_std']:.3f}g")
 
# 计算置信区间
confidence_intervals = qc.calculate_confidence_intervals(sample_data)
 
print(f"\n95%置信区间:")
print(f"缺陷率: [{confidence_intervals['defect_rate'][0]:.4f}, "
      f"{confidence_intervals['defect_rate'][1]:.4f}]")
print(f"长度均值: [{confidence_intervals['length_mean'][0]:.3f}, "
      f"{confidence_intervals['length_mean'][1]:.3f}]mm")
print(f"重量均值: [{confidence_intervals['weight_mean'][0]:.3f}, "
      f"{confidence_intervals['weight_mean'][1]:.3f}]g")
 
# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
 
# 缺陷率分布
axes[0, 0].hist(sample_data['is_defect'], bins=[-0.5, 0.5, 1.5],
               alpha=0.7, edgecolor='black')
axes[0, 0].set_xticks([0, 1])
axes[0, 0].set_xticklabels(['合格', '缺陷'])
axes[0, 0].set_ylabel('频次')
axes[0, 0].set_title(f'样本缺陷分布 (n={sample_size})')
 
# 长度分布
axes[0, 1].hist(sample_data['length'], bins=30, alpha=0.7, edgecolor='black')
axes[0, 1].axvline(100, color='red', linestyle='--', label='目标值')
axes[0, 1].axvline(confidence_intervals['length_mean'][0],
                  color='green', linestyle=':', label='95% CI下限')
axes[0, 1].axvline(confidence_intervals['length_mean'][1],
                  color='green', linestyle=':', label='95% CI上限')
axes[0, 1].set_xlabel('长度 (mm)')
axes[0, 1].set_ylabel('频次')
axes[0, 1].set_title('长度分布')
axes[0, 1].legend()
 
# 重量分布
axes[1, 0].hist(sample_data['weight'], bins=30, alpha=0.7, edgecolor='black')
axes[1, 0].axvline(50, color='red', linestyle='--', label='目标值')
axes[1, 0].axvline(confidence_intervals['weight_mean'][0],
                  color='green', linestyle=':', label='95% CI下限')
axes[1, 0].axvline(confidence_intervals['weight_mean'][1],
                  color='green', linestyle=':', label='95% CI上限')
axes[1, 0].set_xlabel('重量 (g)')
axes[1, 0].set_ylabel('频次')
axes[1, 0].set_title('重量分布')
axes[1, 0].legend()
 
# 散点图：长度vs重量
scatter = axes[1, 1].scatter(sample_data['length'], sample_data['weight'],
                           c=sample_data['is_defect'],
                           cmap='RdYlBu_r', alpha=0.6)
axes[1, 1].set_xlabel('长度 (mm)')
axes[1, 1].set_ylabel('重量 (g)')
axes[1, 1].set_title('长度-重量关系')
plt.colorbar(scatter, ax=axes[1, 1], label='缺陷状态')
 
plt.tight_layout()
plt.show()
 
# 模拟多次抽样的稳定性
print(f"\n=== 抽样稳定性分析 ===")
n_simulations = 100
defect_rates = []
 
for i in range(n_simulations):
    sim_sample = qc.simple_random_sample(sample_size)
    sim_defect_rate = sim_sample['is_defect'].mean()
    defect_rates.append(sim_defect_rate)
 
defect_rates = np.array(defect_rates)
print(f"100次抽样中缺陷率的标准差: {np.std(defect_rates):.4f}")
print(f"95%置信区间宽度: {confidence_intervals['defect_rate'][1] - confidence_intervals['defect_rate'][0]:.4f}")
 
plt.figure(figsize=(10, 6))
plt.hist(defect_rates, bins=20, alpha=0.7, edgecolor='black')
plt.axvline(qc.defect_rate, color='red', linestyle='--',
           label=f'真实缺陷率 ({qc.defect_rate:.4f})')
plt.axvline(np.mean(defect_rates), color='green', linestyle='--',
           label=f'平均估计值 ({np.mean(defect_rates):.4f})')
plt.xlabel('估计缺陷率')
plt.ylabel('频次')
plt.title(f'100次简单随机抽样的缺陷率分布')
plt.legend()
plt.show()

抽样方法
├── 概率抽样
│   ├── 简单随机抽样（SRS）
│   ├── 系统抽样
│   ├── 分层抽样
│   └── 整群抽样
└── 非概率抽样
    ├── 便利抽样
    ├── 判断抽样
    ├── 配额抽样
    └── 雪球抽样

# 简单随机抽样
def simple_random_sampling(population, n):
    return random.sample(population, n)
 
# 系统抽样
def systematic_sampling(population, n):
    N = len(population)
    k = N // n  # 抽样间隔
    start = random.randint(0, k-1)
    return population[start::k]
 
# 分层抽样
def stratified_sampling(population, strata, n):
    sample = []
    for stratum in set(strata):
        stratum_indices = [i for i, s in enumerate(strata) if s == stratum]
        stratum_size = len(stratum_indices)
        stratum_sample_size = int(stratum_size / len(population) * n)
        stratum_sample = random.sample(
            [population[i] for i in stratum_indices],
            stratum_sample_size
        )
        sample.extend(stratum_sample)
    return sample
 
# 整群抽样
def cluster_sampling(clusters, n_clusters):
    return random.sample(clusters, n_clusters)
 
# 性能对比测试
import time
import matplotlib.pyplot as plt
 
population = list(range(1, 100001))  # 10万总体
n = 1000  # 样本量
strata = [random.choice(['A', 'B', 'C']) for _ in range(len(population))]
 
methods = {
    '简单随机抽样': simple_random_sampling,
    '系统抽样': systematic_sampling,
    '分层抽样': stratified_sampling
}
 
times = {}
for name, method in methods.items():
    start = time.time()
    if name == '分层抽样':
        result = method(population, strata, n)
    else:
        result = method(population, n)
    end = time.time()
    times[name] = end - start
    print(f"{name}: {times[name]:.4f}秒, 样本量: {len(result)}")
 
# 绘制对比图
plt.figure(figsize=(10, 6))
methods_list = list(times.keys())
times_list = list(times.values())
plt.bar(methods_list, times_list, alpha=0.7)
plt.ylabel('执行时间 (秒)')
plt.title('不同抽样方法性能对比')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

def calculate_sample_size(population_size, margin_error, confidence_level, estimated_proportion=0.5):
    """
    计算简单随机抽样所需样本量
 
    参数:
    - population_size: 总体规模
    - margin_error: 误差界限
    - confidence_level: 置信水平
    - estimated_proportion: 估计比例
 
    返回:
    - 所需样本量
    """
    from scipy.stats import norm
 
    # Z值
    z_score = norm.ppf((1 + confidence_level) / 2)
 
    # 无限总体样本量
    n_infinite = (z_score**2 * estimated_proportion * (1 - estimated_proportion)) / (margin_error**2)
 
    # 有限总体修正
    n_finite = n_infinite / (1 + (n_infinite - 1) / population_size)
 
    return int(np.ceil(n_finite))
 
# 使用示例
population_sizes = [1000, 10000, 100000, 1000000]
margin_errors = [0.05, 0.03, 0.01]
 
plt.figure(figsize=(12, 8))
 
for i, margin_error in enumerate(margin_errors):
    sample_sizes = [calculate_sample_size(N, margin_error, 0.95) for N in population_sizes]
    plt.subplot(2, 2, i+1)
    plt.plot(population_sizes, sample_sizes, 'o-')
    plt.xlabel('总体规模')
    plt.ylabel('所需样本量')
    plt.title(f'误差界限: ±{margin_error*100}% (95%置信水平)')
    plt.grid(True, alpha=0.3)
    plt.xscale('log')
    plt.yscale('log')
 
plt.tight_layout()
plt.show()

def compare_sampling_methods(population, strata, sample_size=1000, n_simulations=1000):
    """
    比较简单随机抽样与分层抽样的效果
    """
    stratified_means = []
    simple_means = []
 
    for _ in range(n_simulations):
        # 简单随机抽样
        simple_sample = random.sample(population, sample_size)
        simple_means.append(np.mean(simple_sample))
 
        # 分层抽样
        stratified_sample = stratified_sampling(population, strata, sample_size)
        stratified_means.append(np.mean(stratified_sample))
 
    # 计算方差
    simple_var = np.var(simple_means)
    stratified_var = np.var(stratified_means)
    efficiency_gain = (simple_var - stratified_var) / simple_var * 100
 
    # 可视化
    plt.figure(figsize=(12, 5))
 
    plt.subplot(1, 2, 1)
    plt.hist(simple_means, bins=30, alpha=0.7, label='简单随机抽样')
    plt.hist(stratified_means, bins=30, alpha=0.7, label='分层抽样')
    plt.axvline(np.mean(population), color='red', linestyle='--', label='总体均值')
    plt.xlabel('样本均值')
    plt.ylabel('频次')
    plt.title('抽样分布对比')
    plt.legend()
 
    plt.subplot(1, 2, 2)
    methods = ['简单随机抽样', '分层抽样']
    variances = [simple_var, stratified_var]
    plt.bar(methods, variances, alpha=0.7)
    plt.ylabel('样本均值的方差')
    plt.title(f'效率比较 (分层抽样效率提升: {efficiency_gain:.1f}%)')
 
    plt.tight_layout()
    plt.show()
 
    return simple_var, stratified_var, efficiency_gain
 
# 模拟数据
np.random.seed(42)
population = np.concatenate([
    np.random.normal(50, 10, 300),   # 第一层
    np.random.normal(70, 15, 400),   # 第二层
    np.random.normal(90, 20, 300)    # 第三层
])
strata = (['A'] * 300) + (['B'] * 400) + (['C'] * 300)
 
print(f"总体均值: {np.mean(population):.2f}")
print(f"总体标准差: {np.std(population):.2f}")
 
simple_var, stratified_var, efficiency = compare_sampling_methods(
    population, strata, sample_size=300, n_simulations=500
)

def bootstrap_simple_random_sample(population, original_sample, n_bootstrap=1000, confidence_level=0.95):
    """
    使用自助法增强简单随机抽样
    """
    bootstrap_means = []
 
    for _ in range(n_bootstrap):
        # 有放回抽样
        bootstrap_sample = np.random.choice(original_sample,
                                          size=len(original_sample),
                                          replace=True)
        bootstrap_means.append(np.mean(bootstrap_sample))
 
    # 计算自助法置信区间
    alpha = 1 - confidence_level
    lower_percentile = (alpha/2) * 100
    upper_percentile = (1 - alpha/2) * 100
 
    ci_lower = np.percentile(bootstrap_means, lower_percentile)
    ci_upper = np.percentile(bootstrap_means, upper_percentile)
 
    # 可视化
    plt.figure(figsize=(12, 5))
 
    plt.subplot(1, 2, 1)
    plt.hist(bootstrap_means, bins=50, alpha=0.7, edgecolor='black')
    plt.axvline(ci_lower, color='green', linestyle='--', label=f'{confidence_level*100}% CI')
    plt.axvline(ci_upper, color='green', linestyle='--')
    plt.axvline(np.mean(bootstrap_means), color='red', linestyle='--', label='Bootstrap均值')
    plt.xlabel('Bootstrap样本均值')
    plt.ylabel('频次')
    plt.title(f'Bootstrap抽样分布 (n={n_bootstrap})')
    plt.legend()
 
    plt.subplot(1, 2, 2)
    plt.hist(original_sample, bins=30, alpha=0.7, edgecolor='black')
    plt.axvline(np.mean(original_sample), color='red', linestyle='--', label='原始样本均值')
    plt.xlabel('原始样本值')
    plt.ylabel('频次')
    plt.title('原始样本分布')
    plt.legend()
 
    plt.tight_layout()
    plt.show()
 
    return {
        'bootstrap_mean': np.mean(bootstrap_means),
        'bootstrap_std': np.std(bootstrap_means),
        'confidence_interval': (ci_lower, ci_upper),
        'bias': np.mean(bootstrap_means) - np.mean(original_sample)
    }
 
# 使用示例
np.random.seed(42)
population = np.random.normal(100, 15, 10000)
original_sample = np.random.choice(population, size=200, replace=False)
 
bootstrap_results = bootstrap_simple_random_sample(population, original_sample)
 
print("=== Bootstrap分析结果 ===")
print(f"Bootstrap均值: {bootstrap_results['bootstrap_mean']:.3f}")
print(f"Bootstrap标准误: {bootstrap_results['bootstrap_std']:.3f}")
print(f"95%置信区间: [{bootstrap_results['confidence_interval'][0]:.3f}, "
      f"{bootstrap_results['confidence_interval'][1]:.3f}]")
print(f"Bootstrap偏差: {bootstrap_results['bias']:.6f}")

# 展示样本量与总体规模的关系
population_sizes = np.logspace(3, 8, 100)  # 1千到1亿
sample_sizes_005 = [calculate_sample_size(N, 0.05, 0.95) for N in population_sizes]
sample_sizes_001 = [calculate_sample_size(N, 0.01, 0.95) for N in population_sizes]
 
plt.figure(figsize=(10, 6))
plt.plot(population_sizes, sample_sizes_005, label='误差界限 ±5%', linewidth=2)
plt.plot(population_sizes, sample_sizes_001, label='误差界限 ±1%', linewidth=2)
plt.axhline(y=1000, color='red', linestyle='--', label='固定1000样本')
plt.xlabel('总体规模')
plt.ylabel('所需样本量')
plt.title('样本量与总体规模的关系')
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 模拟抽样框偏差
np.random.seed(42)
 
# 真实总体（假设）：收入分布
true_population = np.random.lognormal(10, 1, 100000)
 
# 有偏差的抽样框（排除低收入群体）
biased_frame = true_population[true_population > np.percentile(true_population, 20)]
 
print(f"真实总体均值: {np.mean(true_population):.2f}")
print(f"有偏差抽样框均值: {np.mean(biased_frame):.2f}")
print(f"偏差幅度: {np.mean(biased_frame) - np.mean(true_population):.2f}")
 
# 简单随机抽样结果对比
unbiased_sample = np.random.choice(true_population, 1000, replace=False)
biased_sample = np.random.choice(biased_frame, 1000, replace=False)
 
print(f"\n无偏样本均值: {np.mean(unbiased_sample):.2f}")
print(f"有偏样本均值: {np.mean(biased_sample):.2f}")
print(f"无偏样本误差: {abs(np.mean(unbiased_sample) - np.mean(true_population)):.2f}")
print(f"有偏样本误差: {abs(np.mean(biased_sample) - np.mean(true_population)):.2f}")

# 展示小样本问题
np.random.seed(42)
 
# 非正态分布总体（指数分布）
population = np.random.exponential(2, 10000)
 
# 不同样本量下的抽样分布
sample_sizes = [5, 10, 30, 100]
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()
 
for i, n in enumerate(sample_sizes):
    sample_means = []
    for _ in range(1000):
        sample = np.random.choice(population, n, replace=False)
        sample_means.append(np.mean(sample))
 
    axes[i].hist(sample_means, bins=30, alpha=0.7, edgecolor='black')
    axes[i].set_title(f'样本量 n={n}')
    axes[i].set_xlabel('样本均值')
    axes[i].set_ylabel('频次')
 
    # 绘制理论正态分布曲线
    theoretical_mean = np.mean(population)
    theoretical_std = np.std(population) / np.sqrt(n)
    x = np.linspace(min(sample_means), max(sample_means), 100)
    y = stats.norm.pdf(x, theoretical_mean, theoretical_std)
    axes[i].plot(x, y, 'r-', linewidth=2, label='理论正态分布')
    axes[i].legend()
 
plt.tight_layout()
plt.show()
 
# Shapiro-Wilk正态性检验
from scipy.stats import shapiro
 
print("=== 正态性检验 (Shapiro-Wilk) ===")
for n in sample_sizes:
    sample_means = []
    for _ in range(1000):
        sample = np.random.choice(population, n, replace=False)
        sample_means.append(np.mean(sample))
 
    # 对前100个样本进行检验
    test_sample = sample_means[:100]
    statistic, p_value = shapiro(test_sample)
    print(f"样本量 n={n}: 统计量={statistic:.4f}, p值={p_value:.4f}")

def calculate_design_effect(srs_var, complex_var):
    """计算设计效应（Deff）"""
    return complex_var / srs_var
 
# 模拟复杂抽样设计
np.random.seed(42)
population = np.random.normal(100, 15, 10000)
 
# 简单随机抽样
srs_samples = []
for _ in range(1000):
    sample = np.random.choice(population, 100, replace=False)
    srs_samples.append(np.mean(sample))
 
srs_var = np.var(srs_samples)
 
# 整群抽样（效率较低）
cluster_samples = []
for _ in range(1000):
    # 假设有100个群，每个群100个单位
    clusters = np.array_split(population, 100)
    selected_clusters = np.random.choice(len(clusters), 10, replace=False)
    cluster_sample = np.concatenate([clusters[i] for i in selected_clusters])
    cluster_samples.append(np.mean(cluster_sample))
 
complex_var = np.var(cluster_samples)
deff = calculate_design_effect(srs_var, complex_var)
 
print(f"简单随机抽样方差: {srs_var:.4f}")
print(f"整群抽样方差: {complex_var:.4f}")
print(f"设计效应 (Deff): {deff:.4f}")
print(f"效率损失: {(deff-1)*100:.1f}%")
 
# 可视化
plt.figure(figsize=(12, 5))
 
plt.subplot(1, 2, 1)
plt.hist(srs_samples, bins=30, alpha=0.7, label='简单随机抽样')
plt.hist(cluster_samples, bins=30, alpha=0.7, label='整群抽样')
plt.axvline(np.mean(population), color='red', linestyle='--', label='总体均值')
plt.xlabel('样本均值')
plt.ylabel('频次')
plt.title('抽样方法对比')
plt.legend()
 
plt.subplot(1, 2, 2)
methods = ['简单随机抽样', '整群抽样']
variances = [srs_var, complex_var]
plt.bar(methods, variances, alpha=0.7)
plt.ylabel('样本均值方差')
plt.title('方差对比')
plt.axhline(y=srs_var * deff, color='red', linestyle='--', label=f'Deff = {deff:.2f}')
plt.legend()
 
plt.tight_layout()
plt.show()

class SamplingPlan:
    """抽样计划制定"""
 
    def __init__(self):
        self.objectives = None
        self.target_population = None
        self.sampling_frame = None
        self.sample_size = None
        self.confidence_level = None
        self.margin_error = None
 
    def define_objectives(self, objectives):
        """定义调查目标"""
        self.objectives = objectives
        return self
 
    def define_population(self, population_size, population_characteristics):
        """定义目标总体"""
        self.target_population = {
            'size': population_size,
            'characteristics': population_characteristics
        }
        return self
 
    def assess_sampling_frame(self, frame_completeness, frame_accuracy, frame_currency):
        """评估抽样框质量"""
        # 计算抽样框质量得分
        frame_score = (frame_completeness + frame_accuracy + frame_currency) / 3
 
        recommendations = []
        if frame_score < 0.8:
            recommendations.append("需要改进抽样框质量")
        if frame_completeness < 0.9:
            recommendations.append("补充缺失单位")
        if frame_accuracy < 0.9:
            recommendations.append("更新单位信息")
 
        return {
            'score': frame_score,
            'recommendations': recommendations
        }
 
    def calculate_sample_size(self, margin_error, confidence_level, estimated_proportion=0.5):
        """计算样本量"""
        self.margin_error = margin_error
        self.confidence_level = confidence_level
 
        self.sample_size = calculate_sample_size(
            self.target_population['size'],
            margin_error,
            confidence_level,
            estimated_proportion
        )
        return self
 
    def select_sampling_method(self, population_heterogeneity, cost_constraints, precision_requirements):
        """选择抽样方法"""
        methods = []
 
        # 基于总体同质性
        if population_heterogeneity == 'high':
            methods.append('分层抽样')
 
        # 基于成本约束
        if cost_constraints == 'tight':
            methods.append('整群抽样')
 
        # 基于精度要求
        if precision_requirements == 'high':
            methods.append('简单随机抽样')
 
        return methods
 
    def generate_sampling_protocol(self):
        """生成抽样协议"""
        protocol = {
            'research_objectives': self.objectives,
            'target_population': self.target_population,
            'sample_size': self.sample_size,
            'confidence_level': self.confidence_level,
            'margin_error': self.margin_error,
            'estimated_proportion': 0.5,
            'sampling_method': '简单随机抽样',
            'data_collection_plan': '结构化问卷',
            'quality_control': '双重验证',
            'ethical_considerations': '知情同意、隐私保护'
        }
 
        return protocol
 
# 使用示例
plan = SamplingPlan()
sampling_plan = (plan
    .define_objectives('评估产品市场接受度')
    .define_population(5000000, {'age': '18-65', 'geographic': '全国'})
    .calculate_sample_size(0.03, 0.95)
)
 
protocol = sampling_plan.generate_sampling_protocol()
print("=== 抽样协议 ===")
for key, value in protocol.items():
    print(f"{key}: {value}")

class QualityControlSampler:
    """带质量控制的抽样器"""
 
    def __init__(self, population, sample_size, seed=None):
        self.population = population
        self.sample_size = sample_size
        if seed:
            np.random.seed(seed)
        self.sampling_log = []
 
    def random_sample_with_logging(self):
        """记录抽样过程的随机抽样"""
        # 记录种子
        seed = np.random.randint(0, 2**31)
        np.random.seed(seed)
 
        # 执行抽样
        sample_indices = np.random.choice(
            len(self.population),
            size=self.sample_size,
            replace=False
        )
 
        sample = self.population[sample_indices]
 
        # 记录日志
        log_entry = {
            'timestamp': datetime.now(),
            'seed': seed,
            'sample_size': self.sample_size,
            'population_size': len(self.population),
            'sample_indices': sample_indices.tolist(),
            'sampling_method': 'simple_random_sampling'
        }
        self.sampling_log.append(log_entry)
 
        return sample, log_entry
 
    def validate_sample_representativeness(self, sample):
        """验证样本代表性"""
        # 这里可以添加各种验证指标
        validations = {
            'sample_size_check': len(sample) == self.sample_size,
            'no_duplicate_check': len(np.unique(sample)) == len(sample),
            'within_range_check': all(x in self.population for x in sample)
        }
 
        return validations
 
    def audit_trail(self):
        """生成审计追踪"""
        return pd.DataFrame(self.sampling_log)
 
# 实施质量控制
np.random.seed(42)
population = list(range(1, 10001))
 
qc_sampler = QualityControlSampler(population, 500, seed=42)
sample, log = qc_sampler.random_sample_with_logging()
 
validations = qc_sampler.validate_sample_representativeness(sample)
print("=== 质量控制验证 ===")
for check, result in validations.items():
    print(f"{check}: {'✓ 通过' if result else '✗ 失败'}")
 
audit_trail = qc_sampler.audit_trail()
print(f"\n=== 审计追踪 ===")
print(audit_trail[['timestamp', 'seed', 'sample_size', 'sampling_method']])

class PostStratificationWeighter:
    """后分层加权"""
 
    def __init__(self, sample_data, target_distributions):
        self.sample_data = sample_data
        self.target_distributions = target_distributions
        self.weights = None
 
    def calculate_post_stratification_weights(self, stratify_columns):
        """计算后分层权重"""
        # 计算样本分布
        sample_dist = {}
        for col in stratify_columns:
            sample_dist[col] = self.sample_data[col].value_counts(normalize=True)
 
        # 计算权重
        weights = pd.DataFrame(index=self.sample_data.index)
 
        for col in stratify_columns:
            col_weights = []
            for value in self.sample_data[col]:
                target_prop = self.target_distributions[col].get(value, 0)
                sample_prop = sample_dist[col].get(value, 0)
 
                if sample_prop > 0:
                    weight = target_prop / sample_prop
                else:
                    weight = 1.0
 
                col_weights.append(weight)
 
            weights[col] = col_weights
 
        # 计算综合权重（简单平均）
        self.weights = weights.mean(axis=1)
 
        return self.weights
 
    def apply_weights(self, variable):
        """应用权重到变量"""
        return variable * self.weights
 
    def weighted_analysis(self, variables):
        """加权分析"""
        results = {}
        for var in variables:
            weighted_mean = np.average(self.sample_data[var], weights=self.weights)
            results[var] = weighted_mean
 
        return results
 
# 示例：市场调研后分层
np.random.seed(42)
sample_data = pd.DataFrame({
    'respondent_id': range(1, 1001),
    'age_group': np.random.choice(['18-29', '30-44', '45-59', '60+'], 1000,
                                 p=[0.2, 0.3, 0.3, 0.2]),
    'region': np.random.choice(['东部', '中部', '西部'], 1000,
                              p=[0.4, 0.35, 0.25]),
    'satisfaction': np.random.normal(7.5, 1.5, 1000)
})
 
# 目标总体分布（来自人口普查）
target_distributions = {
    'age_group': {
        '18-29': 0.22,
        '30-44': 0.28,
        '45-59': 0.31,
        '60+': 0.19
    },
    'region': {
        '东部': 0.38,
        '中部': 0.36,
        '西部': 0.26
    }
}
 
# 应用后分层加权
weighter = PostStratificationWeighter(sample_data, target_distributions)
weights = weighter.calculate_post_stratification_weights(['age_group', 'region'])
 
print("=== 后分层加权结果 ===")
print(f"权重统计:")
print(f"  最小值: {weights.min():.3f}")
print(f"  最大值: {weights.max():.3f}")
print(f"  平均值: {weights.mean():.3f}")
print(f"  标准差: {weights.std():.3f}")
 
# 比较加权前后的估计
unweighted_mean = sample_data['satisfaction'].mean()
weighted_mean = weighter.weighted_analysis(['satisfaction'])['satisfaction']
 
print(f"\n满意度估计:")
print(f"  加权前: {unweighted_mean:.3f}")
print(f"  加权后: {weighted_mean:.3f}")
print(f"  差异: {abs(weighted_mean - unweighted_mean):.3f}")
 
# 可视化权重分布
plt.figure(figsize=(12, 5))
 
plt.subplot(1, 2, 1)
plt.hist(weights, bins=30, alpha=0.7, edgecolor='black')
plt.xlabel('权重')
plt.ylabel('频次')
plt.title('后分层权重分布')
 
plt.subplot(1, 2, 2)
plt.scatter(range(len(weights)), weights, alpha=0.6)
plt.xlabel('样本序号')
plt.ylabel('权重')
plt.title('权重分布散点图')
 
plt.tight_layout()
plt.show()

def stratified_train_test_split(X, y, test_size=0.2, random_state=42):
    """分层划分（保持类别比例）"""
    from sklearn.model_selection import train_test_split
 
    return train_test_split(X, y, test_size=test_size,
                           random_state=random_state, stratify=y)
 
def simple_random_cv_split(X, n_splits=5, test_size=0.2, random_state=42):
    """简单随机交叉验证"""
    from sklearn.model_selection import KFold
 
    kfold = KFold(n_splits=n_splits, shuffle=True, random_state=random_state)
    splits = list(kfold.split(X))
 
    return splits
 
# 演示
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
 
# 生成数据
X, y = make_classification(n_samples=1000, n_features=20, n_informative=15,
                          n_redundant=5, n_classes=2, random_state=42)
 
# 简单随机划分
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
                                                   random_state=42)
 
# 训练模型
model = LogisticRegression(random_state=42)
model.fit(X_train, y_train)
 
# 评估
y_pred = model.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
 
print(f"测试集准确率: {accuracy:.3f}")
 
# 交叉验证
from sklearn.model_selection import cross_val_score
 
cv_scores = cross_val_score(model, X, y, cv=5, scoring='accuracy')
print(f"5折交叉验证准确率: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

from collections import Counter
 
def handle_imbalanced_data(X, y, method='undersample', random_state=42):
    """处理不平衡数据"""
    np.random.seed(random_state)
 
    counter = Counter(y)
    print(f"原始类别分布: {counter}")
 
    if method == 'undersample':
        # 简单随机欠采样
        minority_class = min(counter, key=counter.get)
        majority_class = max(counter, key=counter.get)
 
        minority_indices = np.where(y == minority_class)[0]
        majority_indices = np.where(y == majority_class)[0]
 
        np.random.shuffle(majority_indices)
        selected_majority = majority_indices[:len(minority_indices)]
 
        indices = np.concatenate([minority_indices, selected_majority])
        np.random.shuffle(indices)
 
        return X[indices], y[indices]
 
    elif method == 'oversample':
        # 简单随机过采样
        minority_class = min(counter, key=counter.get)
        minority_indices = np.where(y == minority_class)[0]
 
        oversample_indices = np.random.choice(minority_indices,
                                            size=counter[majority_class],
                                            replace=True)
 
        indices = np.concatenate([np.where(y != minority_class)[0], oversample_indices])
        np.random.shuffle(indices)
 
        return X[indices], y[indices]
 
# 生成不平衡数据
X_imbalanced, y_imbalanced = make_classification(n_samples=2000, n_features=10,
                                               n_informative=8, n_redundant=2,
                                               n_classes=2, weights=[0.9, 0.1],
                                               random_state=42)
 
print("=== 不平衡数据处理 ===")
print(f"原始分布: {Counter(y_imbalanced)}")
 
# 欠采样
X_under, y_under = handle_imbalanced_data(X_imbalanced, y_imbalanced, 'undersample')
print(f"欠采样后分布: {Counter(y_under)}")
 
# 过采样
X_over, y_over = handle_imbalanced_data(X_imbalanced, y_imbalanced, 'oversample')
print(f"过采样后分布: {Counter(y_over)}")

# 理想场景示例
ideal_scenarios = [
    {
        'scenario': '大学班级满意度调查',
        'reason': '总体小且同质',
        'sample_size': '全班调查或随机抽取30%',
        'advantage': '实施简单，代表性强'
    },
    {
        'scenario': '工厂产品质量抽检',
        'reason': '产品同质，检测破坏性',
        'sample_size': '2-5%抽检',
        'advantage': '成本低，客观公正'
    },
    {
        'scenario': '临床试验基线特征比较',
        'reason': '需要确保组间可比性',
        'sample_size': '根据统计功效计算',
        'advantage': '保证随机性和无偏性'
    }
]
 
print("=== 理想应用场景 ===")
for scenario in ideal_scenarios:
    print(f"场景: {scenario['scenario']}")
    print(f"原因: {scenario['reason']}")
    print(f"样本量: {scenario['sample_size']}")
    print(f"优势: {scenario['advantage']}\n")

# 展示大总体情况下的抽样效果
population_sizes = [1000, 10000, 100000, 1000000]
sample_size = 1000
true_mean = 50
true_std = 10
 
results = []
 
for N in population_sizes:
    # 生成总体
    population = np.random.normal(true_mean, true_std, N)
 
    # 多次抽样
    sample_means = []
    for _ in range(1000):
        sample = np.random.choice(population, sample_size, replace=False)
        sample_means.append(np.mean(sample))
 
    bias = abs(np.mean(sample_means) - true_mean)
    std_error = np.std(sample_means)
 
    results.append({
        'population_size': N,
        'bias': bias,
        'standard_error': std_error
    })
 
print("=== 大总体抽样效果 ===")
for result in results:
    print(f"总体规模: {result['population_size']:,}")
    print(f"抽样偏差: {result['bias']:.6f}")
    print(f"标准误: {result['standard_error']:.4f}\n")
 
# 可视化
plt.figure(figsize=(12, 5))
 
plt.subplot(1, 2, 1)
pop_sizes = [r['population_size'] for r in results]
biases = [r['bias'] for r in results]
plt.plot(pop_sizes, biases, 'o-')
plt.xlabel('总体规模')
plt.ylabel('抽样偏差')
plt.title('抽样偏差与总体规模的关系')
plt.xscale('log')
plt.grid(True, alpha=0.3)
 
plt.subplot(1, 2, 2)
std_errors = [r['standard_error'] for r in results]
plt.plot(pop_sizes, std_errors, 'o-', label='实际标准误')
plt.axhline(y=true_std/np.sqrt(sample_size), color='red', linestyle='--',
           label=f'理论值 ({true_std/np.sqrt(sample_size):.4f})')
plt.xlabel('总体规模')
plt.ylabel('标准误')
plt.title('标准误与总体规模的关系')
plt.xscale('log')
plt.legend()
plt.grid(True, alpha=0.3)
 
plt.tight_layout()
plt.show()

# 模拟不同质量的抽样框
np.random.seed(42)
 
# 真实总体
true_population = np.random.normal(100, 15, 100000)
 
# 情况1：完美的抽样框
perfect_frame = true_population.copy()
 
# 情况2：不完整的抽样框（缺失10%）
missing_indices = np.random.choice(len(true_population),
                                  size=int(0.1 * len(true_population)),
                                  replace=False)
incomplete_frame = np.delete(true_population, missing_indices)
 
# 情况3：有错误的抽样框（包含错误单位10%）
error_indices = np.random.choice(len(true_population),
                                size=int(0.1 * len(true_population)),
                                replace=False)
error_frame = true_population.copy()
error_frame[error_indices] = np.random.normal(200, 15, len(error_indices))
 
# 情况4：过时的抽样框（过时信息20%）
outdated_indices = np.random.choice(len(true_population),
                                   size=int(0.2 * len(true_population)),
                                   replace=False)
outdated_frame = true_population.copy()
outdated_frame[outdated_indices] = np.random.normal(80, 20, len(outdated_indices))
 
# 比较抽样结果
frames = {
    '完美抽样框': perfect_frame,
    '不完整抽样框': incomplete_frame,
    '错误抽样框': error_frame,
    '过时抽样框': outdated_frame
}
 
print("=== 抽样框质量影响分析 ===")
for name, frame in frames.items():
    samples = []
    for _ in range(100):
        sample = np.random.choice(frame, 1000, replace=False)
        samples.append(np.mean(sample))
 
    mean_error = abs(np.mean(samples) - 100)
    std_error = np.std(samples)
 
    print(f"{name}:")
    print(f"  平均误差: {mean_error:.3f}")
    print(f"  标准误: {std_error:.3f}")
    print(f"  总体均值: {np.mean(frame):.3f}\n")
 
# 可视化
plt.figure(figsize=(12, 10))
 
for i, (name, frame) in enumerate(frames.items()):
    plt.subplot(2, 2, i+1)
    plt.hist(frame, bins=50, alpha=0.7, edgecolor='black')
    plt.axvline(100, color='red', linestyle='--', label='真实均值')
    plt.axvline(np.mean(frame), color='green', linestyle='--', label='框架均值')
    plt.xlabel('数值')
    plt.ylabel('频次')
    plt.title(f'{name}')
    plt.legend()
 
plt.tight_layout()
plt.show()