拒绝采样以从柯西样本生成正态样本

Question

我尝试编写拒绝抽样方法以生成遵循正态分布的样本。这些样本乍一看像是正态分布，但 Shapiro-Wilk 检验的 p 值始终 <0.05。我真的不知道我哪里错了，我只从我的老师那里得到了伪代码（它不是家庭作业）。任何帮助表示赞赏。在我的代码下面：

f <- function(x,m,v) {    #target distribution, m=mean,v=variance
  dnorm(x,m,sqrt(v))
}

g <- function(x,x0,lambda) {  #cauchy distribution for sampling
  dcauchy(x,x0,lambda)
}

genSamp <- function(n,m,v) {  #I want the user to be able to choose mean and sd
                              #and size of the sample
  stProbe <- rep(0,n)         #the sample vector
  interval = c(m-10*sqrt(v),m+10*sqrt(v)) #wanted to go sure that everything
                                          #is covered, so I took a range
                                          #that depends on the mean
  M = max(f(interval,m,v)/g(interval,m,v))  #rescaling coefficient, so the cauchy distribution
                              #is never under the normal distribution
  #I chose x0 = m and lambda = v, so the cauchy distribution is close to a
  #the target normal distribution

  for (i in 1:n) {
    repeat{
      x <- rcauchy(1,m,v)
      u <- runif(1,0,max(f(interval,m,v)))
      if(u < (f(x,m,v)/(M*g(x,m,v)))) {
        break
      }
    }
    stProbe[i] <- x
  }

  return(stProbe)
}

然后我试了一下：

test <- genSamp(100,2,0.5)
hist(test,prob=T,breaks=30)#looked not bad
shapiro.test(test) #p-value way below 0.05

预先感谢您的帮助。

Answer 1

你有

f <- function(x,m,v) {    #target distribution, m=mean,v=variance
  dnorm(x,e,sqrt(v))
}

哪些样本的平均值为 e，但从未定义。

Answer 2

其实我首先检查的是样本均值和样本方差。当我用你的 genSamp 抽取 1000 个样本时，我得到的样本均值为 2，但样本方差约为 2.64，与目标 0.5 相去甚远。

第一个问题是您对 M 的计算。注意：

interval = c(m - 10 * sqrt(v), m + 10 * sqrt(v))

只给你2个值，而不是区间上等距点的网格。在偏离均值 10 个标准差时，正态密度几乎为 0，因此 M 几乎为 0。您需要执行类似

的操作

interval <- seq(m - 10 * sqrt(v), m + 10 * sqrt(v), by = 0.01)

第二个问题是repeat中均匀随机变量的生成。你为什么这样做

u <- runif(1,0,max(f(interval,m,v)))

你想要

u <- runif(1, 0, 1)

通过这些修复，我测试了 genSamp 获得了正确的样本均值和样本方差。样本通过了 Shapiro–Wilk 检验和 Kolmogorov-Smirnov 检验 (?ks.test).

完整的工作代码

f <- function(x,m,v) dnorm(x,m,sqrt(v))

g <- function(x,x0,lambda) dcauchy(x,x0,lambda)

genSamp <- function(n,m,v) {

  stProbe <- rep(0,n)
  interval <- seq(m - 10 * sqrt(v), m + 10 * sqrt(v), by = 0.01)
  M = max(f(interval,m,v)/g(interval,m,v))

  for (i in 1:n) {
    repeat{
      x <- rcauchy(1,m,v)
      u <- runif(1,0,1)
      if(u < (f(x,m,v)/(M*g(x,m,v)))) break
      }
    stProbe[i] <- x
    }

  return(stProbe)
  }

set.seed(0)
test <- genSamp(1000, 2, 0.5)
shapiro.test(test)$p.value
#[1] 0.1563038

ks.test(test, rnorm(1000, 2, sqrt(0.5)))$p.value
#[1] 0.7590978

拒绝采样以从柯西样本生成正态样本

Rejection Sampling to generate Normal samples from Cauchy samples

r

statistics

normal-distribution