在数据科学和机器学习领域，概率论和统计学扮演着至关重要的角色。Python作为一种强大而灵活的编程语言，提供了丰富的库和工具来实现这些概念。本文将通过20个Python实例，展示如何在实际应用中运用概率论和统计学知识。
1. 基本概率计算
让我们从一个简单的硬币投掷实验开始：
import random
def coin_flip(n):
return [random.choice(['H', 'T']) for _ in range(n)]
flips = coin_flip(1000)
probability_head = flips.count('H') / len(flips)
print(f"Probability of getting heads: {probability_head:.2f}")
这个例子模拟了1000次硬币投掷，并计算出现正面的概率。
2. 描述性统计
使用NumPy和Pandas来计算一些基本的描述性统计量：
import numpy as np
import pandas as pd
data = np.random.normal(0, 1, 1000)
df = pd.DataFrame(data, columns=['values'])
print(df.describe())
这个例子生成了1000个服从标准正态分布的随机数，并计算了均值、标准差等统计量。
3. 概率分布
使用SciPy绘制正态分布的概率密度函数：
import scipy.stats as stats
import matplotlib.pyplot as plt
x = np.linspace(-5, 5, 100)
plt.plot(x, stats.norm.pdf(x, 0, 1))
plt.title("Standard Normal Distribution")
plt.xlabel("x")
plt.ylabel("Probability Density")
plt.show()
4. 中心极限定理
演示中心极限定理：
sample_means = [np.mean(np.random.exponential(1, 100)) for _ in range(1000)]
plt.hist(sample_means, bins=30, edgecolor='black')
plt.title("Distribution of Sample Means")
plt.xlabel("Sample Mean")
plt.ylabel("Frequency")
plt.show()
这个例子展示了指数分布的样本均值趋向于正态分布。
5. 假设检验
进行t检验：
from scipy import stats
group1 = np.random.normal(0, 1, 100)
group2 = np.random.normal(0.5, 1, 100)
t_statistic, p_value = stats.ttest_ind(group1, group2)
print(f"T-statistic: {t_statistic:.4f}")
print(f"P-value: {p_value:.4f}")
这个例子比较两组数据，检验它们的均值是否有显著差异。
6. 置信区间
计算均值的置信区间：
data = np.random.normal(0, 1, 100)
mean = np.mean(data)
se = stats.sem(data)
ci = stats.t.interval(0.95, len(data)-1, loc=mean, scale=se)
print(f"95% Confidence Interval: {ci}")
7. 线性回归
使用sklearn进行简单线性回归：
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
X = np.random.rand(100, 1)
y = 2 * X + 1 + np.random.randn(100, 1) * 0.1
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
model = LinearRegression()
model.fit(X_train, y_train)
print(f"Coefficient: {model.coef_[0][0]:.2f}")
print(f"Intercept: {model.intercept_[0]:.2f}")
8. 多项式回归
使用numpy的polyfit函数进行多项式回归：
x = np.linspace(0, 1, 100)
y = x**2 + np.random.randn(100) * 0.1
coeffs = np.polyfit(x, y, 2)
p = np.poly1d(coeffs)
plt.scatter(x, y)
plt.plot(x, p(x), color='red')
plt.title("Polynomial Regression")
plt.show()
9. 贝叶斯推断
使用PyMC3进行简单的贝叶斯推断：
import pymc3 as pm
with pm.Model() as model:
mu = pm.Normal('mu', mu=0, sd=1)
obs = pm.Normal('obs', mu=mu, sd=1, observed=np.random.randn(100))
trace = pm.sample(1000)
pm.plot_posterior(trace)
plt.show()
这个例子展示了如何对正态分布的均值进行贝叶斯推断。
10. 蒙特卡罗模拟
使用蒙特卡罗方法估算π：
def estimate_pi(n):
inside_circle = 0
total_points = n
for _ in range(total_points):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
if x**2 + y**2 <= 1:
inside_circle += 1
return 4 * inside_circle / total_points
print(f"Estimated value of π: {estimate_pi(1000000):.6f}")
这个例子通过随机点的方法估算π的值。
11. 马尔可夫链
实现简单的马尔可夫链：
states = ['A', 'B', 'C']
transition_matrix = {
'A': {'A': 0.3, 'B': 0.6, 'C': 0.1},
'B': {'A': 0.4, 'B': 0.2, 'C': 0.4},
'C': {'A': 0.1, 'B': 0.3, 'C': 0.6}
}
def next_state(current):
return random.choices(states, weights=list(transition_matrix[current].values()))[0]
current = 'A'
for _ in range(10):
print(current, end=' -> ')
current = next_state(current)
print(current)
12. 主成分分析 (PCA)
使用sklearn进行PCA：
from sklearn.decomposition import PCA
data = np.random.randn(100, 5)
pca = PCA(n_components=2)
reduced_data = pca.fit_transform(data)
plt.scatter(reduced_data[:, 0], reduced_data[:, 1])
plt.title("PCA Reduced Data")
plt.xlabel("First Principal Component")
plt.ylabel("Second Principal Component")
plt.show()
13. 时间序列分析
使用statsmodels进行ARIMA模型拟合：
from statsmodels.tsa.arima.model import ARIMA
np.random.seed(1)
ts = pd.Series(np.cumsum(np.random.randn(100)))
model = ARIMA(ts, order=(1,1,1))
results = model.fit()
print(results.summary())
14. 核密度估计
使用seaborn进行核密度估计：
import seaborn as sns
data = np.concatenate([np.random.normal(-2, 1, 1000), np.random.normal(2, 1, 1000)])
sns.kdeplot(data)
plt.title("Kernel Density Estimation")
plt.show()
15. Bootstrap方法
使用Bootstrap方法估计均值的置信区间：
def bootstrap_mean(data, num_samples, size):
means = [np.mean(np.random.choice(data, size=size)) for _ in range(num_samples)]
return np.percentile(means, [2.5, 97.5])
data = np.random.normal(0, 1, 1000)
ci = bootstrap_mean(data, 10000, len(data))
print(f"95% CI for the mean: {ci}")
16. 假设检验的功效分析
进行t检验的功效分析：
from statsmodels.stats.power import TTestIndPower
effect = 0.5
alpha = 0.05
power = 0.8
analysis = TTestIndPower()
sample_size = analysis.solve_power(effect, power=power, nobs1=None, ratio=1.0, alpha=alpha)
print(f"Required sample size: {sample_size:.0f}")
17. 贝叶斯信息准则 (BIC)
使用BIC进行模型选择：
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
X = np.random.rand(100, 3)
y = X[:, 0] + 2*X[:, 1] + np.random.randn(100) * 0.1
def bic(y, y_pred, n_params):
mse = mean_squared_error(y, y_pred)
return len(y) * np.log(mse) + n_params * np.log(len(y))
models = [
LinearRegression().fit(X[:, :1], y),
LinearRegression().fit(X[:, :2], y),
LinearRegression().fit(X, y)
]
bic_scores = [bic(y, model.predict(X[:, :i+1]), i+1) for i, model in enumerate(models)]
best_model = np.argmin(bic_scores)
print(f"Best model (lowest BIC): {best_model + 1} features")
18. 非参数检验
使用Mann-Whitney U检验：
group1 = np.random.normal(0, 1, 100)
group2 = np.random.normal(0.5, 1, 100)
statistic, p_value = stats.mannwhitneyu(group1, group2)
print(f"Mann-Whitney U statistic: {statistic}")
print(f"P-value: {p_value:.4f}")
19. 生存分析
使用lifelines进行Kaplan-Meier生存分析：
from lifelines import KaplanMeierFitter
T = np.random.exponential(10, size=100)
E = np.random.binomial(1, 0.7, size=100)
kmf = KaplanMeierFitter()
kmf.fit(T, E, label="KM Estimate")
kmf.plot()
plt.title("Kaplan-Meier Survival Curve")
plt.show()
20. 聚类分析
使用K-means聚类：
from sklearn.cluster import KMeans
X = np.random.randn(300, 2)
kmeans = KMeans(n_clusters=3)
kmeans.fit(X)
plt.scatter(X[:, 0], X[:, 1], c=kmeans.labels_)
plt.title("K-means Clustering")
plt.show()