嘿，针对你遇到的两个问题，我整理了具体的解决思路和代码调整方案，帮你一步步排查和验证：

如果绘制的直方图和对应的密度曲线（不管是核密度还是理论正态密度）不匹配，可以从这几个角度检查：

• 刻度一致性问题：直方图默认展示的是频数，而密度函数是概率密度，两者纵坐标尺度不统一。解决方法是绘制直方图时加上freq=FALSE参数，让直方图的纵坐标变为密度，这样才能和密度曲线对齐，比如：hist(est, freq=FALSE)。
• 参数准确性问题：如果是叠加理论正态密度，要确保用样本的均值和方差来计算密度值，比如dnorm(x, mean=mean(est), sd=sd(est))，不能用默认的标准正态参数（均值0，方差1），否则肯定不匹配。
• 样本量不足：你代码里设置的M=1000，样本量不算特别大，直方图的波动会比较明显。可以尝试增大M（比如调到5000或10000），样本量越大，直方图越接近真实分布，和密度曲线的贴合度会更好。

首先先修正你代码里的一个关键问题：循环内的变量M和外层定义的M<-1000重名了，这会导致循环提前终止，结果完全错误。我把循环内的变量改成了M_hat，避免命名冲突。

接下来，正态性检验可以从可视化和统计检验两个维度进行：

可视化方法（直观判断） #

• Q-Q图：这是最常用的正态性可视化工具，代码如下：

qqnorm(iest)
qqline(iest, col="red", lwd=2)

如果图里的点基本沿着红色直线分布，说明样本近似正态；如果点明显偏离直线，尤其是两端，说明不符合正态性。

• 直方图+双密度曲线：同时绘制核密度曲线和理论正态密度曲线，对比贴合度：

hist(iest, freq=FALSE, main="Histogram of Estimates with Density Fits", xlab="Estimated Values")
# 核密度曲线
lines(density(iest), col="blue", lwd=2)
# 理论正态密度曲线（用样本均值和方差）
x_seq <- seq(min(iest), max(iest), length.out=100)
lines(x_seq, dnorm(x_seq, mean=mean(iest), sd=sd(iest)), col="red", lwd=2, lty=2)
legend("topright", legend=c("Kernel Density", "Normal Density"), col=c("blue", "red"), lwd=2, lty=c(1,2))

统计检验方法（量化判断） #

• Shapiro-Wilk检验：适合样本量≤5000的场景，代码：

shapiro.test(iest)

如果输出的p值大于你的显著性水平（通常是0.05），则不能拒绝“样本来自正态分布”的假设；反之则说明不符合正态性。

• Kolmogorov-Smirnov检验：检验样本分布和理论正态分布的拟合程度，代码：

ks.test(iest, "pnorm", mean=mean(iest), sd=sd(iest))

同样，p值大于0.05时支持正态性假设。

修正后的完整R代码 #

# 初始化参数
y0 <- 50
a <- 10
b <- 0.4
theta <- 4
T <- 30
h <- 0.01
t <- seq(0, T, h)
n <- length(t)
M <- 1000  # 外层实验次数，避免和循环内变量重名
iest <- rep(0, M)

for(j in 1:M) {
  # 生成布朗运动路径
  W <- c(0, cumsum(rnorm(n-1, 0, 1)*sqrt(h)))
  y <- rep(y0, n)
  # 求解SDE
  for(i in 1:(n-1)) {
    y[i+1] <- y[i] + a*t[i]*h - b*y[i]*h + theta*(W[i+1]-W[i])
  }
  delta <- 1/n
  X <- y
  # 重命名变量，避免和外层M冲突
  M_hat <- (a/b) + (X - (a/b))*exp(-b*delta)
  # 计算估计值
  iest[j] <- sqrt( (2*b)/(1-exp(-2*b*delta)) * (1/n)*sum(X[2:n]-M_hat[1:(n-1)]) )
}

# 正态性检验流程
# 1. Q-Q图
qqnorm(iest)
qqline(iest, col="red", lwd=2)

# 2. 直方图+密度曲线
hist(iest, freq=FALSE, main="Histogram of Estimates with Density Fits", xlab="Estimated Values")
lines(density(iest), col="blue", lwd=2)
x_seq <- seq(min(iest), max(iest), length.out=100)
lines(x_seq, dnorm(x_seq, mean=mean(iest), sd=sd(iest)), col="red", lwd=2, lty=2)
legend("topright", legend=c("Kernel Density", "Normal Density"), col=c("blue", "red"), lwd=2, lty=c(1,2))

# 3. Shapiro-Wilk检验
cat("Shapiro-Wilk正态性检验结果：\n")
print(shapiro.test(iest))

# 4. Kolmogorov-Smirnov检验
cat("\nKolmogorov-Smirnov正态性检验结果：\n")
print(ks.test(iest, "pnorm", mean=mean(iest), sd=sd(iest)))

内容的提问来源于stack exchange，提问作者user8810618