import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 设置样式
sns.set_theme(style="whitegrid")

def plot_discrete_rv(values, pmf, cdf, samples, rv_name):
    # 准备经验 CDF（ECDF）
    n_samples = len(samples)
    sorted_samples = np.sort(samples)
    ecdf_x = np.unique(sorted_samples)
    ecdf_y = [np.sum(sorted_samples <= x) / n_samples for x in ecdf_x]
    plt.figure(figsize=(15, 8))

    # 理论 PMF
    plt.subplot(2, 2, 1)
    plt.stem(values, pmf, basefmt=" ", linefmt='-.')
    plt.title(f"Theoretical PMF: {rv_name}")
    plt.xlabel("x")
    plt.ylabel("f(X=x) = P(X = x)")
    plt.ylim(0, 1.1)

    # 理论 CDF
    plt.subplot(2, 2, 2)
    plt.step(values, cdf, where='post', color='green')
    plt.title(f"Theoretical CDF: {rv_name}")
    plt.xlabel("x")
    plt.ylabel("F(x) = P(X ≤ x)")
    plt.ylim(0, 1.1)
    plt.grid(True)

    # 采样直方图
    
    plt.subplot(2, 2, 3)
    sns.countplot(x=samples, hue=samples, legend=False, palette='pastel', stat='proportion', order=values)
    plt.title(f"Empirical Distribution ({n_samples} samples)")
    plt.xlabel("x")
    plt.ylabel("Relative Frequency")

    # 经验 CDF（ECDF）
    plt.subplot(2, 2, 4)
    plt.step(ecdf_x, ecdf_y, where='post', color='orange')
    plt.title("Empirical CDF")
    plt.xlabel("x")
    plt.ylabel("ECDF")
    plt.ylim(0, 1.1)
    plt.grid(True)
    plt.tight_layout()
    plt.show()

import numpy as np

# 设置随机种子以保证结果可重复
np.random.seed(42)

# 1. 定义参数
a, b = 1, 6  # 均匀分布的范围
values = np.arange(a, b+1)  # 离散取值：1~6
n = len(values)
pmf = np.ones(n) / n  # 每个值的概率均等
cdf = np.cumsum(pmf)  # 累积分布函数

# 2. 抽样
n_samples = 1000
samples = np.random.choice(values, size=n_samples, p=pmf)


# 3. 可视化：理论分布 + 采样频率对比
plot_discrete_rv(values, pmf, cdf, samples, f"DiscreteUniform({a},{b})")

import random

def discrete_uniform_sample(a, b, n):
    N = b - a + 1
    # U ~ Uniform(0, 1)
    U = [random.random() for _ in range(n)] # random.random() 返回区间 [0.0, 1.0) 的浮点数
    # X ~ Discrete Uniform(a, b)
    X = [a+int(u*N) for u in U]  # When u = 0.999, u*N = 0.999 * N, which is close to N-1, so a + int(u*N) will be b.
    return U, X


discrete_uniform_sample(0, 1, 10)

[random.randint(0, 1) for _ in range(10)]  # 使用内置函数直接采样，验证结果是否正确

import random
import matplotlib.pyplot as plt
from math import comb

# 采样
a, b, N = 1, 10, 30000
origin_samples, samples = discrete_uniform_sample(a, b, N)
print(f"Empirical mean: {sum(samples)/len(samples):.3f}")
# 统计频率
counts = [samples.count(k) / N for k in range(a, b+1)]

# 计算理论概率（PMF）
theoretical = [1/(b-a+1) for _ in range(a, b+1)]
print(f"PMF = {theoretical}")

# Step 5: 可视化
plt.figure(figsize=(10, 6))
plt.bar(range(a, b+1), counts, width=0.4, label='Sampled Frequency', color='skyblue', align='center')
plt.bar(range(a, b+1), theoretical, width=0.4, label='Theoretical PMF', color='orange', align='edge')
plt.xlabel("Value")
plt.ylabel("Probability")
plt.title(f"Discrete Uniform Sample ({a}~{b})")
plt.legend()
plt.grid(True)
plt.show()

Step 1：生成一个 U ~ Uniform(0, 1)
Step 2：通过线性变换 X = a + (b - a) * U
Step 3：X 就是你要的 sample from Uniform(a, b)

import random

def sample_uniform(a, b):
    U = random.random()        # U ~ Uniform(0,1)
    X = a + (b - a) * U        # X ~ Uniform(a, b)
    return X


def sample_uniform_list(a, b, n):
    return [sample_uniform(a, b) for _ in range(n)]

sample_uniform_list(0, 1, 10)

random.uniform(0, 1)

import matplotlib.pyplot as plt

a, b, n = 0, 1, 10000
pdf = 1 / (b - a)  # 均匀分布的概率密度函数
samples = sample_uniform_list(a, b, n)

plt.figure(figsize=(15, 10))
# 采样结果
plt.subplot(2, 2, 1)
plt.scatter(range(n), samples, alpha=0.5, color='blue')
plt.title(f"Sample results from Uniform({a}, {b})")
plt.xlabel("n")
plt.ylabel("Sample Value (X)")
plt.ylim(0, 1.1)

# 采样直方图（频率）
N = 10
plt.subplot(2, 2, 3)
plt.hist(samples, bins=N, density=False, edgecolor='black')
plt.title(f"Histogram of Sample from Uniform({a}, {b})")
plt.xlabel("Sample Value (X)")
plt.ylabel("Absolute Frequency")

# 采样直方图（相对频率）
plt.subplot(2, 2, 4)
plt.hist(samples, bins=N, density=True, edgecolor='black')
plt.hlines(pdf, a, b, colors='red', linestyles='solid', label='PDF')
plt.title(f"Histogram of Sample from Uniform({a}, {b})")
plt.xlabel("Sample Value (X)")
plt.ylabel(f"h = Relative Frequency * N (N={N})")

plt.show()

import numpy as np

# 设置随机种子以保证结果可重复
np.random.seed(42)

# 1. 定义参数
p = 0.7
values = np.array([0, 1])  # 离散取值：1~6
n = len(values)
pmf = np.array([1-p, p])
cdf = np.cumsum(pmf)  # 累积分布函数

# 2. 抽样
n_samples = 1000
samples = np.array([np.random.binomial(n=1, p=p) for _ in range(n_samples)])


# 3. 可视化：理论分布 + 采样频率对比
plot_discrete_rv(values, pmf, cdf, samples, f"Bernoulli(p={p})")

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# -----------------------------
# 参数设置
# -----------------------------
p = 0.3                # 成功概率
n_trials = 10000        # 总实验次数
interval = 50          # 动画间隔时间（ms）

# -----------------------------
# 生成伯努利实验数据
# -----------------------------
np.random.seed(0)
samples = np.random.binomial(n=1, p=p, size=n_trials)
cumulative_freq = np.cumsum(samples) / np.arange(1, n_trials + 1)

# -----------------------------
# 创建图形
# -----------------------------
fig, ax = plt.subplots(figsize=(8, 4))
ax.set_ylim(0, 1)
ax.set_xlim(1, n_trials)
ax.axhline(y=p, color='red', linestyle='--', label=f'True Probability p = {p}')
line, = ax.plot([], [], lw=2, label='Empirical Frequency')
text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
ax.set_xlabel('Number of Trials')
ax.set_ylabel('Success Frequency')
ax.set_title('Bernoulli Theorem Animation')
ax.legend()

# -----------------------------
# 动画更新函数
# -----------------------------
def update(frame):
    x = np.arange(1, frame + 1)
    y = cumulative_freq[:frame]
    line.set_data(x, y)
    text.set_text(f'n = {frame}, freq = {y[-1]:.3f}')
    return line, text

# -----------------------------
# 创建动画
# -----------------------------
ani = FuncAnimation(fig, update, frames=np.arange(1, n_trials + 1, 10),
                    interval=interval, blit=True)
# Save as GIF
from matplotlib.animation import PillowWriter
ani.save("bernoulli_theorem_animation.gif", writer=PillowWriter(fps=10))

plt.close(fig)

from IPython.display import Image
Image(filename="bernoulli_theorem_animation.gif")

#plt.close()
# 在 notebook 中显示动画
#HTML(ani.to_jshtml())

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
from IPython.display import HTML

# Parameters
p = 0.3
n_frames = 100

# Pre-generate uniform samples
uniform_samples = np.random.uniform(0, 1, n_frames)
bernoulli_samples = (uniform_samples < p).astype(int)

# Set up the figure and axes
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
uniform_ax = axs[0]
bernoulli_ax = axs[1]

# Initialize plots
uniform_scatter = uniform_ax.scatter([], [], c='blue', alpha=0.6)
bernoulli_bar = bernoulli_ax.bar([0, 1], [0, 0], color='orange', edgecolor='black')

# Configure uniform axis
uniform_ax.axvline(p, color='red', linestyle='--', label=f'p = {p}')
uniform_ax.set_xlim(0, 1)
uniform_ax.set_ylim(0, 1)
uniform_ax.set_title('Sampling from Uniform(0,1)')
uniform_ax.set_xlabel('Value')
uniform_ax.set_ylabel('Random Height')
uniform_ax.legend()
uniform_ax.grid(True)

# Configure Bernoulli axis
bernoulli_ax.set_xlim(-0.5, 1.5)
bernoulli_ax.set_ylim(0, n_frames)
bernoulli_ax.set_title('Bernoulli Sample Counts')
bernoulli_ax.set_xlabel('Value')
bernoulli_ax.set_ylabel('Count')
bernoulli_ax.set_xticks([0, 1])
bernoulli_ax.grid(True)

# Store sample counts
count_0 = 0
count_1 = 0
x_vals = []
y_vals = []

# Animation update function
def update(frame):
    global count_0, count_1, x_vals, y_vals

    u = uniform_samples[frame]
    x_vals.append(u)
    y_vals.append(np.random.rand())  # random y position for scatter

    bern_sample = bernoulli_samples[frame]
    if bern_sample == 0:
        count_0 += 1
    else:
        count_1 += 1

    # Update scatter
    uniform_scatter.set_offsets(np.column_stack((x_vals, y_vals)))

    # Update bar chart
    bernoulli_bar[0].set_height(count_0)
    bernoulli_bar[1].set_height(count_1)

    return uniform_scatter, bernoulli_bar

# Create animation
ani = animation.FuncAnimation(fig, update, frames=n_frames, interval=100, blit=False)
# 保存为 GIF
from matplotlib.animation import PillowWriter
ani.save("sample_uniform_to_bernoulli.gif", writer=PillowWriter(fps=10))

plt.close(fig)

from IPython.display import Image
Image(filename="sample_uniform_to_bernoulli.gif")
#plt.close()
# 在 notebook 中显示动画
#HTML(ani.to_jshtml())

import random

def sample_bernoulli(p):
    U = random.random() # U ~ Uniform(0, 1)
    return 1 if U < p else 0

sample_bernoulli(0.7)

p, N = 0.7, 10000
samples = [sample_bernoulli(p) for _ in range(N)]
print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 0.7

import numpy as np

np.random.binomial(n=1, p=0.7)

p, N = 0.7, 10000
samples = [np.random.binomial(n=1, p=p) for _ in range(N)]
print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 0.7

from scipy.stats import binom
import numpy as np

# 设置随机种子以保证结果可重复
np.random.seed(42)

# 1. 定义参数
n, p = 5, 0.3 # X ～ Binomial(n, p)
values = range(n + 1)  # 离散取值：0, 1, ..., n
pmf = binom.pmf(values, n, p)
cdf = np.cumsum(pmf)  # 累积分布函数

# 2. 抽样
n_samples = 1000
samples = np.random.binomial(n=n, p=p, size=n_samples)

# 3. 可视化：理论分布 + 采样频率对比
plot_discrete_rv(values, pmf, cdf, samples, f'Binomial(n={n}, p={p})')

Binomial(n, p)
│
├── 小 n: 直接伯努利重复
├── 小 p: 用 Poisson(np)
├── 大 n, np(1-p)>10: 用正态(np, np(1-p))
├── 任意 n,p: 反函数 or 查表
└── 高效工业级: BTPE 算法（NumPy）

import matplotlib.pyplot as plt
from math import comb
def verify_binomial_sample(n, p, samples):
    """
    验证二项分布采样的正确性
    :param n: 二项分布的试验次数
    :param p: 成功概率
    :param samples: 采样结果列表
    :return: None
    """
    N = len(samples)
    empirical_mean = sum(samples) / N
    theoretical_mean = n * p
    print(f"Empirical mean: {empirical_mean:.3f}, Theoretical mean: {theoretical_mean:.3f}")
    
    empirical_variance = sum((x - empirical_mean) ** 2 for x in samples) / len(samples)
    theoretical_variance = n * p * (1 - p)
    print(f"Empirical variance: {empirical_variance:.3f}, Theoretical variance: {theoretical_variance:.3f}")

    # 统计频率
    counts = [samples.count(k) / N for k in range(n + 1)]
    # 计算理论概率（PMF）
    theoretical = [comb(n, k) * (p ** k) * ((1 - p) ** (n - k)) for k in range(n + 1)]
    # 可视化
    plt.figure(figsize=(10, 6))
    plt.bar(range(n + 1), counts, width=0.4, label='Sampled Frequency', color='skyblue', align='edge')
    plt.bar([k - 0.4 for k in range(n + 1)], theoretical, width=0.4, label='Theoretical PMF', color='orange', align='edge')
    plt.xlabel("Number of Successes")
    plt.ylabel("Probability")
    plt.title(f"Binomial Distribution: n={n}, p={p}")
    plt.legend()
    plt.grid(True)
    plt.show()

import numpy as np
np.random.binomial(n=10, p=0.5, size=1000)

import random

def sample_binomial_mimic(n, p):
    count = 0
    for _ in range(n): # repeat N times Bernouli sampling
        u = random.random()  # 从 [0,1) 中采样一个均匀变量
        if 1-p <= u < 1:
            count += 1       # 成功就加一
    return count

def sample_binomial_mimic_list(n, p, num_samples):
    return [sample_binomial_mimic(n, p) for _ in range(num_samples)]

n, p, num_samples = 10, 0.5, 10
sample_binomial_mimic_list(n, p, num_samples)

import matplotlib.pyplot as plt

n, p, num_samples = 10, 0.3, 1000
samples = sample_binomial_mimic_list(n, p, num_samples)
plt.hist(samples, bins=range(12), align='left', rwidth=0.8, color='skyblue')
plt.xlabel("Number of Successes")
plt.ylabel("Frequency")
plt.title("Sampling from Binomial(n=10, p=0.3)")
plt.show()

# 采样
n = 10
p = 0.3
N = 10000  # 采样次数
samples = sample_binomial_mimic_list(n, p, N)
#print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 np = 10x0.4 = 4.0

verify_binomial_sample(n, p, samples)

import matplotlib.pyplot as plt
from math import comb
import numpy as np
import random

# 参数
n = 10
p = 0.4

# 计算 PMF 和 CDF
k_vals = list(range(n + 1))
pmf = [comb(n, k) * p**k * (1 - p)**(n - k) for k in k_vals]
cdf = np.cumsum(pmf)

####### 采样 ####### 
# 生成一个随机数 u
u = 0.52  # 也可以改成 random.random()

# 找到第一个 cdf[k] >= u
k_selected = None
for k, value in enumerate(cdf):
    if u <= value:
        k_selected = k
        break

####### 绘制验证 ####### 
# 绘图
plt.figure(figsize=(10, 6))

# PMF 图（柱状图）
plt.bar(k_vals, pmf, alpha=0.6, label='PMF: P(X = k)', color='skyblue')

# CDF 图（阶梯线）
plt.step(k_vals, cdf, where='mid', color='orange', label='CDF: P(X ≤ k)', linewidth=2)

# 画出 u 和对应的 k_selected
plt.axhline(y=u, color='red', linestyle='--', label=f'u = {u:.2f}')
plt.axvline(x=k_selected, color='green', linestyle='--', label=f'sampled k = {k_selected}')

# 标注
plt.title(f"Inverse Transform Sampling for Binomial(n={n}, p={p})")
plt.xlabel("k")
plt.ylabel("Probability")
plt.legend()
plt.grid(True)
plt.xticks(k_vals)
plt.ylim(0, 1.05)

plt.show()

import random
from math import comb

def sample_binomial_inverse(n, p):
    # Step 1: generate a uniform random number
    u = random.random()
    
    # Step 2: initialize cumulative probability
    cumulative = 0.0

    for k in range(n + 1):
        prob = comb(n, k) * (p ** k) * ((1 - p) ** (n - k))
        cumulative += prob
        if u <= cumulative:
            return k

    return n  # fallback

def sample_binomial_inverse_list(n, p, num_samples):
    return [sample_binomial_inverse(n, p) for _ in range(num_samples)]

n, p, num_samples = 10, 0.5, 10
sample_binomial_inverse_list(n, p, num_samples)

import matplotlib.pyplot as plt

n, p, num_samples = 10, 0.3, 1000
samples = sample_binomial_inverse_list(n, p, num_samples)

plt.hist(samples, bins=range(12), align='left', rwidth=0.8, color='skyblue')
plt.xlabel("Number of Successes")
plt.ylabel("Frequency")
plt.title("Sampling from Binomial(n=10, p=0.3)")
plt.show()

# 采样
n = 10
p = 0.3
N = 10000  # 采样次数
samples = sample_binomial_inverse_list(n, p, N)
#print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 np = 10x0.4 = 4.0

verify_binomial_sample(n, p, samples)

import random
from math import comb

binomial_table, binomial_table_n, binomial_table_p = None, None, None
def generate_binomial_table(n, p):
    """
    生成二项分布的 CDF 查找表
    :param n: 二项分布的试验次数
    :param p: 成功概率
    :return: None
    """
    global binomial_table
    binomial_table = {}
    cumulative = 0.0
    for k in range(n + 1):
        prob = comb(n, k) * (p ** k) * ((1 - p) ** (n - k))
        cumulative += prob
        binomial_table[k] = cumulative

def sample_binomial_table(n, p):
    # Step 1: generate a uniform random number
    u = random.random()
    
    # Step 2: look up in the CDF table
    if binomial_table is None or binomial_table_n != n or binomial_table_p != p:
        generate_binomial_table(n, p)
    # Step 3: find the first k such that CDF[k] >= u
    for k in range(n + 1):
        if u <= binomial_table[k]:
            return k
    # If no k found, return n
    return n  # fallback

def sample_binomial_table_list(n, p, num_samples):
    return [sample_binomial_table(n, p) for _ in range(num_samples)]

# 采样
n = 10
p = 0.3
N = 10000  # 采样次数
samples = sample_binomial_table_list(n, p, N)
#print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 np = 10x0.4 = 4.0

verify_binomial_sample(n, p, samples)

import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt

def rejection_sample_binomial(n, p, N_samples=1000):
    samples = []

    # 目标分布：Bin(n, p)
    bin_pmf = [binom.pmf(k, n, p) for k in range(n+1)]

    # Proposal 分布：Uniform(0, n)
    g = 1.0 / (n + 1)

    # 找最大值用于计算 M
    M = max(bin_pmf) / g

    while len(samples) < N_samples:
        k = np.random.randint(0, n+1)
        u = np.random.uniform(0, 1)

        accept_prob = bin_pmf[k] / (M * g)
        if u < accept_prob:
            samples.append(k)

    return samples

# 验证
# 采样
n = 10
p = 0.3
N = 10000  # 采样次数
samples = rejection_sample_binomial(n, p, N)
#print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 np = 10x0.4 = 4.0

verify_binomial_sample(n, p, samples)

import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt

def normal_approx_binomial(n, p, N_samples=10000):
    mu = n * p
    sigma = np.sqrt(n * p * (1 - p))
    
    # 采样
    samples = np.random.normal(loc=mu, scale=sigma, size=N_samples)
    
    # 四舍五入并裁剪
    samples = np.round(samples).astype(int)
    samples = np.clip(samples, 0, n)
    
    return samples.tolist()

# 验证
# 采样
n = 10
p = 0.3
N = 10000  # 采样次数
samples = normal_approx_binomial(n, p, N)
#print(f"Empirical mean: {sum(samples)/len(samples):.3f}")  # 应该接近 np = 10x0.4 = 4.0

verify_binomial_sample(n, p, samples)

import numpy as np
import matplotlib.pyplot as plt

# 设置随机种子以保证可重复性
np.random.seed(42)

# 指数分布参数 λ（速率参数）
lambda_val = 12

# 采样数量
n_samples = 10000

# 使用 numpy 采样
samples = np.random.exponential(scale=1/lambda_val, size=n_samples)

# 创建图形
plt.figure(figsize=(8, 5))

# 绘制采样的直方图（归一化为密度）
plt.hist(samples, bins=50, density=True, alpha=0.6, color='skyblue', label='Sampled Histogram')

# 生成理论曲线
x_vals = np.linspace(0, np.max(samples), 200)
pdf = lambda_val * np.exp(-lambda_val * x_vals)
plt.plot(x_vals, pdf, 'r-', lw=2, label=f'Theoretical PDF (λ={lambda_val})')
a, b = 2, 1
y_vals = a+b*x_vals
ypdf = a*pdf+b
plt.plot(y_vals, pdf, 'b-', lw=2, label=f'Theoretical PDF (λ={lambda_val};y={a}*x+{b})')

# 图形设置
plt.title('Exponential Distribution Visualization')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import random
import math

def sample_exponential_inverse(lambda_val):
    """
    采样指数分布
    :param lambda_val: 指数分布的速率参数 λ
    :param n_samples: 采样数量
    :return: 采样结果
    """
    u = random.random()  # 从 [0, 1) 中采样一个均匀变量
    return -math.log(1-u) / lambda_val  # 使用逆变换采样公式 X = -ln(1-U)/λ
    
def sample_exponential_inverse_list(lambda_val, n_samples):
    return [sample_exponential_inverse(lambda_val) for _ in range(n_samples)]

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
lambda_val = 1.5
n = 10000

# 使用反函数法采样
X = sample_exponential_inverse_list(lambda_val, n)

# 验证：绘图
plt.hist(X, bins=100, density=True, alpha=0.6, label='Sampled Histogram')
x_vals = np.linspace(0, 5, 200)
plt.plot(x_vals, lambda_val * np.exp(-lambda_val * x_vals), 'r-', label='True PDF')
plt.xlabel('X')
plt.ylabel('PDF')
plt.title('Exponential Distribution (λ = {:.1f}, n = {})'.format(lambda_val, n))
plt.legend()
plt.show()

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

# 设置随机种子以保证结果可重复
np.random.seed(42)
# 1. 定义参数
p = 0.3
values = np.arange(1, 15)
n = len(values)
pmf = stats.geom(p).pmf(values)
cdf = np.cumsum(pmf)  # 累积分布函数

# 2. 抽样
n_samples = 1000
samples = np.random.geometric(p=p, size=n_samples)

# 3. 可视化：理论分布 + 采样频率对比
plot_discrete_rv(values, pmf, cdf, samples, f"Geometric(p={p})")

import matplotlib.pyplot as plt
from scipy import stats

def verify_geometric_sample(p, samples):
    """
    验证几何分布采样的正确性
    :param p: 成功概率
    :param samples: 采样结果列表（值域应为 1, 2, 3, ...）
    :return: None
    """
    N = len(samples)

    # 均值与方差验证
    empirical_mean = sum(samples) / N
    theoretical_mean = 1 / p
    print(f"Empirical mean: {empirical_mean:.3f}, Theoretical mean: {theoretical_mean:.3f}")
    
    empirical_variance = sum((x - empirical_mean) ** 2 for x in samples) / N
    theoretical_variance = (1 - p) / (p ** 2)
    print(f"Empirical variance: {empirical_variance:.3f}, Theoretical variance: {theoretical_variance:.3f}")

    # 统计最大值作为绘图范围
    k_max = max(samples)
    ks = list(range(1, k_max + 1))

    # 统计频率（经验 PMF）
    from collections import Counter
    counter = Counter(samples)
    counts = [counter.get(k, 0) / N for k in ks]

    # 理论 PMF
    theoretical = stats.geom(p).pmf(ks)

    # 可视化
    plt.figure(figsize=(10, 6))
    plt.bar([k - 0.2 for k in ks], counts, width=0.4, label='Sampled Frequency', color='skyblue', align='center')
    plt.bar([k + 0.2 for k in ks], theoretical, width=0.4, label='Theoretical PMF', color='orange', align='center')
    plt.xlabel("Number of Trials Until First Success (k)")
    plt.ylabel("Probability")
    plt.title(f"Geometric Distribution Verification (p={p})")
    plt.legend()
    plt.grid(True)
    plt.show()

import numpy as np

p = 0.3
n = 10000  # 采样次数
samples = np.random.geometric(p=p, size=n).tolist()
verify_geometric_sample(p, samples)

def sample_geometric_inverse(p):
    """
    使用逆变换法采样几何分布
    :param p: 成功概率
    :return: 采样结果
    """
    u = np.random.uniform(0, 1)  # 从 [0, 1) 中采样一个均匀变量
    return int(np.ceil(np.log(1 - u) / np.log(1 - p)))  # 使用逆变换公式

def sample_geometric_inverse_list(p, n_samples):
    """
    生成几何分布的采样列表
    :param p: 成功概率
    :param n_samples: 采样数量
    :return: 采样结果列表
    """
    return [sample_geometric_inverse(p) for _ in range(n_samples)]

p = 0.3
n = 10000  # 采样次数
samples = sample_geometric_inverse_list(p, n)
verify_geometric_sample(p, samples)

import random

def build_cdf_table(p, max_k=100):
    cdf = []
    total = 0.0
    for k in range(1, max_k + 1):
        prob = (1 - p) ** (k - 1) * p
        total += prob
        cdf.append(total)
    return cdf

def sample_geometric_lookup(cdf_table):
    u = random.random()
    for i, value in enumerate(cdf_table):
        if u <= value:
            return i + 1
    return len(cdf_table)  # fallback

def sample_geometric_lookup_list(p, n_samples, max_k=100):
    cdf_table = build_cdf_table(p, max_k)
    return [sample_geometric_lookup(cdf_table) for _ in range(n_samples)]

p = 0.3
n = 10000  # 采样次数
samples = sample_geometric_lookup_list(p, n)
verify_geometric_sample(p, samples)

import random

def sample_geometric_mimic(p):
    """
    模拟几何分布采样
    :param p: 成功概率
    :return: 采样结果
    """
    k = 1
    while True:
        u = random.random()  # 从 [0, 1) 中采样一个均匀变量
        if 1-p <=u <=1:  # 成功
            return k
        k += 1  # 失败，增加计数

def sample_geometric_mimic_list(p, n_samples):
    """
    生成几何分布的采样列表
    :param p: 成功概率
    :param n_samples: 采样数量
    :return: 采样结果列表
    """
    return [sample_geometric_mimic(p) for _ in range(n_samples)]

# 采样
p = 0.3
N = 10000  # 采样次数
samples = sample_geometric_mimic_list(p, N)

verify_geometric_sample(p, samples)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

x = np.linspace(-4, 4, 500)
pdf = norm.pdf(x, loc=0, scale=1)

plt.plot(x, pdf, label='Standard Normal PDF', color='blue')
plt.fill_between(x, pdf, alpha=0.3, color='skyblue')
plt.title("Standard Normal Distribution")
plt.xlabel("x")
plt.ylabel("Probability Density")
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def verify_standard_normal_sample(samples):
    """
    验证标准正态分布采样的正确性
    :param samples: 从 N(0,1) 中采样的结果 (list or array)
    """
    N = len(samples)
    
    # 经验统计量
    empirical_mean = np.mean(samples)
    empirical_var = np.var(samples)
    
    # 理论统计量
    theoretical_mean = 0
    theoretical_var = 1

    print(f"Empirical mean: {empirical_mean:.4f} | Theoretical mean: {theoretical_mean}")
    print(f"Empirical var: {empirical_var:.4f} | Theoretical var: {theoretical_var}")
    
    # 直方图 + 理论 PDF
    plt.figure(figsize=(12, 5))

    plt.subplot(1, 2, 1)
    x = np.linspace(-4, 4, 500)
    plt.hist(samples, bins=50, density=True, alpha=0.6, color='skyblue', label='Sampled Histogram')
    plt.plot(x, norm.pdf(x, loc=theoretical_mean, scale=np.sqrt(theoretical_var)), 'r-', lw=2, label='Theoretical PDF')
    plt.title('PDF Comparison')
    plt.xlabel('x')
    plt.ylabel('Density')
    plt.legend()
    plt.grid(True)

    # 经验 CDF vs 理论 CDF
    plt.subplot(1, 2, 2)
    sorted_samples = np.sort(samples)
    empirical_cdf = np.arange(1, N+1) / N
    plt.plot(sorted_samples, empirical_cdf, label='Empirical CDF', lw=2)
    plt.plot(x, norm.cdf(x, loc=theoretical_mean, scale=np.sqrt(theoretical_var)), 'r--', label='Theoretical CDF')
    plt.title('CDF Comparison')
    plt.xlabel('x')
    plt.ylabel('CDF')
    plt.legend()
    plt.grid(True)

    plt.tight_layout()
    plt.show()

import numpy as np

samples = np.random.normal(loc=0, scale=1, size=10000)
# verify the sampling results
verify_standard_normal_sample(samples)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def box_muller_sample(n_samples=10000, visualize=True):
    # Step 1: 生成两个均匀分布变量
    U1 = np.random.uniform(0, 1, n_samples)
    U2 = np.random.uniform(0, 1, n_samples)

    # Step 2: Box-Muller变换
    R = np.sqrt(-2 * np.log(U1))
    theta = 2 * np.pi * U2
    Z1 = R * np.cos(theta)
    Z2 = R * np.sin(theta)

    if visualize:
        num_bins = 100
        plt.figure(figsize=(12, 6))
        plt.subplot(1, 2, 1)
        plt.hist(Z1, bins=num_bins, density=True, alpha=0.6, label='Z1')
        x = np.linspace(-4, 4, 500)
        plt.plot(x, norm.pdf(x), 'r--', label='Standard Normal PDF')
        plt.title("Z1 ~ N(0,1)")
        plt.legend()

        plt.subplot(1, 2, 2)
        plt.hist(Z2, bins=num_bins, density=True, alpha=0.6, label='Z2', color='orange')
        plt.plot(x, norm.pdf(x), 'r--', label='Standard Normal PDF')
        plt.title("Z2 ~ N(0,1)")
        plt.legend()

        plt.suptitle("Box-Muller Sampling of Standard Normal Distribution")
        plt.grid(True)
        plt.show()

    return Z1, Z2

box_muller_sample(n_samples=10000, visualize=True)

import numpy as np

def sample_normal_via_clt(a=0, b=1, n=12, num_samples=10000):
    """
    利用中心极限定理从 Uniform(a, b) 采样生成近似标准正态分布的样本
    :param a: Uniform 分布下界
    :param b: Uniform 分布上界
    :param n: 每次采样的 Uniform 个数
    :param num_samples: 总共生成多少个样本
    :return: 一个近似标准正态的样本数组
    """
    # 计算均值和标准差
    mu = (a + b) / 2
    sigma = (b - a) / np.sqrt(12) # 均匀分布的标准差
    # Step 1: 从 Uniform(a, b) 中生成 n 个样本，重复 num_samples 次
    uniform_samples = np.random.uniform(a, b, size=(num_samples, n))
    #print(uniform_samples.shape) #(num_samples, n)
    # Step 2: 对每一组样本求和并标准化
    sample_sums = np.sum(uniform_samples, axis=1) # (num_samples,)
    z_samples = (sample_sums - n * mu) / (np.sqrt(n) * sigma)
    return z_samples

samples = sample_normal_via_clt(a=0, b=1, n=10, num_samples=100000)
verify_standard_normal_sample(samples)

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, laplace

# 目标分布: 标准正态
target_pdf = lambda x: norm.pdf(x)
# proposal: Laplace 分布
proposal_pdf = lambda x: laplace.pdf(x)
proposal_sampler = lambda size: laplace.rvs(size=size)

# M 是 f(x)/q(x) 的上界，估计为最大值（略大于实际最大）
x_vals = np.linspace(-10, 10, 1000)
M = np.max(target_pdf(x_vals) / proposal_pdf(x_vals)) * 1.1

# 采样函数
def rejection_sample(n):
    samples = []
    while len(samples) < n:
        x = proposal_sampler(1)[0]
        u = np.random.uniform()
        if u < target_pdf(x) / (M * proposal_pdf(x)):
            samples.append(x)
    return np.array(samples)

# 生成样本
samples = rejection_sample(10000)

# 可视化
x = np.linspace(-4, 4, 1000)
plt.figure(figsize=(10,6))
plt.hist(samples, bins=50, density=True, alpha=0.5, label="Sampled (AR)")
plt.plot(x, norm.pdf(x), label="Target N(0,1)", lw=2)
plt.plot(x, M * proposal_pdf(x), '--', label=f"M * Proposal (M={M:.2f})", color="red")
plt.title("Rejection Sampling from N(0,1)")
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
lambda_ = 1  # parameter for Exponential distribution
a = 2
b = 3

# Define the original PDF of X ~ Exp(lambda)
def f_X(x):
    return lambda_ * np.exp(-lambda_ * x) * (x >= 0)

# Define the transformed PDF of Y = aX + b
def f_Y(y):
    x = (y - b) / a
    return f_X(x) * (1 / abs(a)) * (y >= b)

# Create x and y ranges
x_vals = np.linspace(0, 10, 400)
y_vals = np.linspace(3, 20, 400)

# Evaluate PDFs
fx_vals = f_X(x_vals)
fy_vals = f_Y(y_vals)

# Plot
plt.figure(figsize=(10, 5))
plt.plot(x_vals, fx_vals, label=r'$f_X(x)$ (Exponential)', color='blue')
plt.plot(y_vals, fy_vals, label=r'$f_Y(y)$ (Transformed)', color='orange')
plt.title('PDF of X and Y = aX + b (a=2, b=3)')
plt.xlabel('x or y')
plt.ylabel('Density')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the original PDF: X ~ Uniform(0, 1)
def f_X(x):
    return np.ones_like(x) * ((x >= 0) & (x <= 1))

# Define the transformation: Y = sqrt(X) => X = Y^2
def f_Y(y):
    return 2 * y * ((y >= 0) & (y <= 1))

# Create value ranges
x_vals = np.linspace(-0.2, 1.2, 400)
y_vals = np.linspace(-0.2, 1.2, 400)

# Evaluate PDFs
fx_vals = f_X(x_vals)
fy_vals = f_Y(y_vals)

# Plot the PDFs
plt.figure(figsize=(10, 5))
plt.plot(x_vals, fx_vals, label=r'$f_X(x)$: Uniform(0, 1)', color='blue')
plt.plot(y_vals, fy_vals, label=r'$f_Y(y)$: $Y = \sqrt{X}$', color='orange')
plt.title("""PDF Transformation: $Y = \sqrt{X}$ with $X \sim U(0, 1)$""")
plt.xlabel('x or y')
plt.ylabel('Density')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定义非线性函数和它的导数
def g(x):
    return np.log(x)

def g_prime(x):
    return 1 / x

# 取均值点 mu
mu = 1.0
x = np.linspace(0.5, 2.0, 400)

# 原函数值
y = g(x)

# 一阶线性近似：在 mu 附近
y_approx = g(mu) + g_prime(mu) * (x - mu)

# 绘图
plt.figure(figsize=(8, 6))
plt.plot(x, y, label=r'$g(x) = \log(x)$', color='blue')
plt.plot(x, y_approx, '--', label='Linear approximation at $\mu=1$', color='red')
plt.axvline(mu, color='gray', linestyle=':', label=r'$\mu_X$')

# 标记点
plt.scatter([mu], [g(mu)], color='black', zorder=5)
plt.text(mu+0.02, g(mu)+0.1, r'$g(\mu_X)$', fontsize=12)

plt.title('Function $g(x)$ and Its Linear Approximation at $\mu_X$')
plt.xlabel('x')
plt.ylabel('g(x)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# 假设原始分布 X 是指数分布（不是正态）
X = np.random.exponential(scale=1.0, size=(10000,))  # 原始样本

# 不同样本数量下的 sample mean
sample_sizes = [1, 5, 10, 30, 100]
means = []

for n in sample_sizes:
    sample_means = [np.mean(np.random.choice(X, n, replace=False)) for _ in range(1000)]
    means.append(sample_means)

# 可视化不同 n 下的 sample mean 分布
fig, axs = plt.subplots(1, len(sample_sizes), figsize=(18, 3))

for i, n in enumerate(sample_sizes):
    axs[i].hist(means[i], bins=30, color='skyblue', edgecolor='black', density=True)
    axs[i].set_title(f'n={n}')
    axs[i].axvline(np.mean(X), color='red', linestyle='--', label='True Mean')
    axs[i].legend()

plt.suptitle('Sample Mean Distribution as Sample Size Increases')
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 被积函数
def f(x):
    return np.exp(-x**2)

# 积分区间
a, b = 0, 1

# Monte Carlo 采样数量
n = 10000
x_samples = np.random.uniform(a, b, n) # 1️⃣ 均匀采样 [a, b] 区间

# 样本均值估计
estimate = (b - a) * np.mean(f(x_samples)) # 2️⃣ 计算样本均值并乘以区间长度

print(f"[Estimate Value] Monte Carlo estimate of ∫₀¹ e^(-x²) dx ≈ {estimate:.6f}")

# 计算真实值
from scipy.special import erf
true_val = np.sqrt(np.pi)/2 * erf(1)
print(f"[True Value] The true value of ∫₀¹ e^(-x²) dx ≈ {true_val:.6f}")

import numpy as np
import plotly.graph_objs as go
from plotly.subplots import make_subplots

# 目标分布：指数分布 X ~ Exp(λ=1)
x_vals = np.linspace(0, 6, 500)
cdf_vals = 1 - np.exp(-x_vals)  # CDF of Exp(1)
inv_cdf = lambda u: -np.log(1 - u)  # x = 反函数 F⁻¹(u)

# 随机采样 5 个均匀数
np.random.seed(42)
u_samples = np.sort(np.random.uniform(0, 1, 5))
x_samples = inv_cdf(u_samples)

# 创建子图：左侧显示 CDF 映射过程，右侧显示 Histogram
fig = make_subplots(rows=1, cols=2,
                    subplot_titles=("Inverse Transform Sampling", "Histogram of Transformed Samples"),
                    column_widths=[0.6, 0.4])

# 左图：CDF 曲线
fig.add_trace(go.Scatter(x=x_vals, y=cdf_vals, mode='lines', name='CDF F(x)', line=dict(color='blue')),
              row=1, col=1)

# 添加每个 u 值的水平线和对应的 x 映射
for u, x in zip(u_samples, x_samples):
    fig.add_trace(go.Scatter(x=[0, x], y=[u, u], mode='lines',
                             line=dict(dash='dot', color='gray'), showlegend=False), row=1, col=1)
    fig.add_trace(go.Scatter(x=[x], y=[u], mode='markers+text',
                             marker=dict(color='red', size=8),
                             text=[f"x={x:.2f}"], textposition='top right', showlegend=False), row=1, col=1)

# 添加 u 样本点
fig.add_trace(go.Scatter(x=[0]*len(u_samples), y=u_samples, mode='markers',
                         marker=dict(symbol='line-ns-open', color='green', size=10),
                         name='u ~ Uniform(0,1)'), row=1, col=1)

# 右图：x 样本的直方图
fig.add_trace(go.Histogram(x=x_samples, nbinsx=10, name='Sampled X', marker_color='orange'), row=1, col=2)

# 布局设置
fig.update_layout(height=500, width=900, title_text="Inverse Transform Sampling Visualization (Exponential RV)",
                  showlegend=True)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_yaxes(title_text="F(x) or u", row=1, col=1)
fig.update_xaxes(title_text="Sampled x", row=1, col=2)
fig.update_yaxes(title_text="Frequency", row=1, col=2)

fig.show()

import numpy as np
import matplotlib.pyplot as plt

# Target distribution: Exponential(lambda=1)
from scipy.stats import expon

# Generate Uniform samples
n = 1000
U = np.random.uniform(0, 1, n)
U_sorted = np.sort(U)

# Compute inverse CDF (quantile function) for exponential
Y = expon.ppf(U_sorted)

# Prepare plot
fig, axs = plt.subplots(1, 2, figsize=(14, 5))

# Plot CDF of Exponential
x = np.linspace(0, 5, 500)
cdf = expon.cdf(x)

axs[0].plot(x, cdf, label='CDF of Exponential(λ=1)', color='blue')
axs[0].scatter(Y, U_sorted, color='red', alpha=0.7, label='(Y, U)')
axs[0].set_title('Inverse Transform Sampling Visualization')
axs[0].set_xlabel('y')
axs[0].set_ylabel('u = F(y)')
axs[0].legend()
axs[0].grid(True)

# Plot histogram of sampled Y
axs[1].hist(Y, bins=50, density=True, alpha=0.7, color='orange', label='Sampled Y')
axs[1].plot(x, expon.pdf(x), color='blue', lw=2, label='True PDF')
axs[1].set_title('Sampled Distribution via Inverse Transform')
axs[1].set_xlabel('y')
axs[1].set_ylabel('Density')
axs[1].legend()
axs[1].grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import random

def flatten_joint_distribution(x_vals, y_vals, joint_probs):
    """
    展平二维联合分布矩阵
    :param joint_probs: 二维联合分布矩阵
    :return: 展平后的概率分布
    """
    # 展平联合分布
    flattened = joint_probs.flatten()
    #print("Flatten joint probabilities:", flattened)

    # 所有可能的(x, y)组合（顺序与flattened一致）
    xy_pairs = [(x, y) for x in x_vals for y in y_vals]
    #print("XY pairs:", xy_pairs)

    return flattened, xy_pairs


def sample_2d_discrete_inverse_transform(x_vals, y_vals, joint_probs):
    # flatten the joint distribution
    flattened, xy_pairs = flatten_joint_distribution(x_vals, y_vals, joint_probs)
    # apply inverse transform sampling
    cdf = np.cumsum(flattened)
    u = random.random()
    for i, threshold in enumerate(cdf):
        if u <= threshold:
            return xy_pairs[i]
    return xy_pairs[-1]  # fallback for u ~ 1.0

import numpy as np
import matplotlib.pyplot as plt

# 定义二维离散联合分布
x_vals = [0, 1, 2]
y_vals = [0, 1, 2]

joint_probs = np.array([
    [0.1, 0.2, 0.1],
    [0.1, 0.2, 0.1],
    [0.05, 0.05, 0.1]
])  # shape (3, 3)

assert np.isclose(joint_probs.sum(), 1.0), "Joint probabilities must sum to 1"

# 采样
N = 1000  # 采样次数
samples = [sample_2d_discrete_inverse_transform(x_vals, y_vals, joint_probs) for _ in range(N)]
x_sample, y_sample = zip(*samples)

# 可视化联合频率
import seaborn as sns
import pandas as pd

df = pd.DataFrame({'x': x_sample, 'y': y_sample})
pivot_table = pd.crosstab(df['x'], df['y'], normalize='all')

plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
sns.heatmap(joint_probs, annot=True, cmap="YlGnBu", fmt=".3f", xticklabels=y_vals, yticklabels=x_vals)
plt.title("Theoretical Joint Distribution from Samples")
plt.xlabel("Y")
plt.ylabel("X")

plt.subplot(1, 2, 2)
sns.heatmap(pivot_table, annot=True, cmap="YlGnBu", fmt=".3f", xticklabels=y_vals, yticklabels=x_vals)
plt.title("Empirical Joint Distribution from Samples")
plt.xlabel("Y")
plt.ylabel("X")
plt.show()

import numpy as np
import random

def sample_2d_marginal_conditional(x_vals, y_vals, joint_probs):
    # 计算边缘分布 P(X)
    P_X = joint_probs.sum(axis=1)  # shape (3,)

    # 计算条件分布 P(Y | X)
    P_Y_given_X = joint_probs / P_X[:, np.newaxis]  # 每行除以对应 P(X)
    # 1. 采样 X
    x = random.choices(x_vals, weights=P_X)[0]
    x_idx = x_vals.index(x)

    # 2. 根据 P(Y|X=x) 采样 Y
    y = random.choices(y_vals, weights=P_Y_given_X[x_idx])[0]
    return (x, y)

import numpy as np
import matplotlib.pyplot as plt
import random

# 值域
x_vals = [0, 1, 2]
y_vals = [0, 1, 2]

# 联合概率矩阵
joint_probs = np.array([
    [0.1, 0.2, 0.1],
    [0.1, 0.2, 0.1],
    [0.05, 0.05, 0.1]
])

N = 1000  # 采样次数
# 采样
samples = [sample_2d_marginal_conditional(x_vals, y_vals, joint_probs) for _ in range(N)]
x_sample, y_sample = zip(*samples)

# 可视化
import seaborn as sns
import pandas as pd

df = pd.DataFrame({'x': x_sample, 'y': y_sample})
pivot_table = pd.crosstab(df['x'], df['y'], normalize='all')

plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
sns.heatmap(joint_probs, annot=True, cmap="YlGnBu", fmt=".3f", xticklabels=y_vals, yticklabels=x_vals)
plt.title("Theoretical Joint Distribution from Marginal + Conditional Sampling")
plt.xlabel("Y")
plt.ylabel("X")

plt.subplot(1, 2, 2)
sns.heatmap(pivot_table, annot=True, cmap="YlGnBu", fmt=".3f", xticklabels=y_vals, yticklabels=x_vals)
plt.title("Empirical Joint Distribution from Marginal + Conditional Sampling")
plt.xlabel("Y")
plt.ylabel("X")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 设置参数
N = 6 # 树的深度（也即将 [0,1] 划分为 2^6 = 64 个区间）

def draw_path(ax, level, x0, x1, path=[]):
    if level > N:
        x_mid = (x0 + x1) / 2
        y = -level
        ax.plot(x_mid, y, 'ko')  # 画叶子节点
        bin_label = ''.join(str(b) for b in path)
        ax.text(x_mid, y - 0.2, f"{bin_label}", ha='center', fontsize=8)
        return
    
    x_mid = (x0 + x1) / 2
    y = -level
    ax.plot(x_mid, y, 'ko')  # 当前节点
    
    # 左边
    draw_path(ax, level + 1, x0, x_mid, path + [0])
    ax.plot([x_mid, (x0 + x_mid) / 2], [y, y - 1], 'gray', lw=1)

    # 右边
    draw_path(ax, level + 1, x_mid, x1, path + [1])
    ax.plot([x_mid, (x1 + x_mid) / 2], [y, y - 1], 'gray', lw=1)

# 画图
fig, ax = plt.subplots(figsize=(12, 6))
draw_path(ax, 0, 0.0, 1.0)
ax.set_ylim(-N - 1.5, 1)
ax.set_xlim(-0.05, 1.05)
ax.set_yticks([])
ax.set_xlabel("Interval on [0,1]")
ax.set_title("Grid Sampling: Binary Tree Partition (N=6)")
plt.grid(True, axis='x', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.show()