X∼N(0,1),f(X)= 1+X 2 1−X 2 ​ ，用采样的方法求E(f(X))连续随机变量的采样python

在Python中，我们可以使用numpy库来进行高斯分布（正态分布）的随机生成，并通过蒙特卡洛模拟的方法估计函数的期望值。对于给定的连续随机变量 \( X \sim N(0,1) \)，其密度函数 \( f(x) = (1 + x^2) / (1 - x^2) \)，由于x的取值范围受限于单位圆内(-1到1)，我们通常需要对这个函数做一些处理。

首先，让我们导入必要的库并定义函数：

import numpy as np

def pdf_x(x):
    return (1 + x**2) / (1 - x**2)

# 定义一个辅助函数，确保x在(-1, 1)范围内
def within_unit_circle(x):
    return -1 < x < 1

def estimate_expectation(n_samples=1000000):
    # 生成n_samples个标准正态分布的随机数
    x_samples = np.random.normal(size=n_samples)
    
    # 过滤掉超出单位圆的样本
    valid_samples = x_samples[within_unit_circle(x_samples)]
    
    # 计算并返回期望值的估算
    expectation_estimate = np.mean(pdf_x(valid_samples))
    return expectation_estimate

# 调用函数并打印结果
estimated_expectation = estimate_expectation()
print("估计的期望值为:", estimated_expectation)

这里的`estimate_expectation`函数会生成大量标准正态分布的点，然后只计算那些落在单位圆内的点的函数值的平均，作为函数\( f(X) \)的期望值的近似。由于采样的数量越多，估计的结果越接近真实值。你可以调整`n_samples`来获得更精确的估计。