На самом деле, первое, что я проверил, это среднее значение выборки и дисперсия выборки. Когда я рисую 1000 выборок с вашим genSamp
, я получаю среднее значение выборки на уровне 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
получает правильное среднее значение выборки и дисперсию выборки. Образцы проходят как тест Шапиро – Вилка, так и тест Колмогорова-Смирнова (?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