Интеграция Метрополис-Гастингс - почему не работает моя стратегия?

Предположим, у меня есть функция которую я хочу интегрировать Конечно, предполагая, что стремится к нулю в конечных точках, нет сбоев, хорошая функция. Один из способов , который я теребил является использование алгоритма Метрополиса-Гастингса , чтобы создать список образцов от распределения пропорционально к , который отсутствует константа нормировки который я назову , а затем вычислю некоторую статистику по этим : ∫ ∞ - ∞ g ( x ) d x . г ( х ) х 1 , х 2 , … , х $g(x)$

$$\int_{-\infty}^\infty g(x) dx.$$ $g(x)$$x_1, x_2, \dots, x_n$N = ∫ ∞ - ∞ g ( x ) d x p ( x ) f ( x ) x $g(x)$ $$N = \int_{-\infty}^{\infty} g(x)dx$$ $p(x)$$f(x)$$x$ $$\frac{1}{n} \sum_{i=0}^n f(x_i) \approx \int_{-\infty}^\infty f(x)p(x)dx.$$

Поскольку , я могу заменить в f (x) = U (x) / g (x), чтобы отменить g из интеграла, что приведет к выражению вида \ frac {1} {N} \ int _ {- \ infty} ^ {\ infty} \ frac {U (x)} {g (x)} g (x) dx = \ frac {1} {N} \ int _ {- \ infty} ^ \ infty U (x) dx. Таким образом, при условии, что U (x) интегрируется в 1 вдоль этой области, я должен получить результат 1 / N , который я мог бы просто взять обратно, чтобы получить желаемый ответ. Поэтому я мог бы взять диапазон моей выборки (чтобы наиболее эффективно использовать точки) r = x_ \ max - x_ \ min и позволить U (x) = 1 / r для каждой выборки, которую я нарисовал. Таким образом, U (х)$p(x) = g(x)/N$g $f(x) = U(x)/g(x)$$g$U(x)11/Nr=xmax-xminU(x)=1/rU(x

$$\frac{1}{N}\int_{-\infty}^{\infty}\frac{U(x)}{g(x)} g(x) dx = \frac{1}{N}\int_{-\infty}^\infty U(x) dx.$$ $U(x)$$1$$1/N$$r = x_\max - x_\min$$U(x) = 1/r$$U(x)$оценивает в ноль за пределами региона, где мои выборки не, но интегрируется в $1$ в этом регионе. Поэтому, если я теперь возьму ожидаемое значение, я должен получить: $$E\left [\frac{U(x)}{g(x)}\right ] = \frac{1}{N} \approx \frac{1}{n} \sum_{i=0}^n \frac{U(x)}{g(x)}.$$

Я попытался проверить это в R для примера функции . В этом случае я не использую Metropolis-Hastings для генерации выборок, но использую фактические вероятности для генерации выборок (просто для тестирования). Я не совсем получаю результаты, которые я ищу. По сути, полное выражение того, что я бы вычислял, таково: В моей теории это должно быть равно . Это близко, но это, конечно, не сходится ожидаемым образом, я делаю что-то не так?$g(x) = e^{-x^2}$rnorm1/

$$\frac{1}{n(x_{\max} - x_\min)} \sum_{i=0}^n \frac{1}{ e^{-x_i^2}}.$$ $1/\sqrt{\pi}$

ys = rnorm(1000000, 0, 1/sqrt(2))
r = max(ys) - min(ys)
sum(sapply(ys, function(x) 1/( r * exp(-x^2))))/length(ys)
## evaluates to 0.6019741. 1/sqrt(pi) = 0.5641896

Редактировать для CliffAB

Причина, по которой я использую диапазон, заключается в простом определении функции, которая не равна нулю в области, где находятся мои точки, но которая интегрируется в в диапазоне . Полная спецификация функции: Мне не нужно было использовать качестве этой равномерной плотности. Я мог бы использовать некоторую другую плотность, интегрированную в , например, плотность вероятности Однако это сделало бы суммирование отдельных образцов тривиальным, т.е. [ - ∞ , ∞ ] U ( x ) = { $1$$[-\infty, \infty]$U(х

$$U(x) = \begin{cases} \frac{1}{x_\max - x_\min} & x_\max > x > x_\min \\ 0 & \text{otherwise.} \end{cases}$$ $U(x)$P ( x ) = $1$ $$P(x) = \frac{1}{\sqrt{\pi}} e^{-x^2}.$$ $$\frac{1}{n} \sum_{i=0}^n \frac{P(x)}{g(x)} = \frac{1}{n} \sum_{i=0}^n \frac{e^{-x_i^2}/\sqrt{\pi}}{e^{-x_i^2} } = \frac{1}{n} \sum_{i=0}^n \frac{1}{\sqrt{\pi}} = \frac{1}{\sqrt{\pi}}.$$

Я мог бы попробовать эту технику для других дистрибутивов, которые интегрируются в . Тем не менее, я все еще хотел бы знать, почему это не работает для равномерного распределения.$1$

$g$$g$$g$

$$\int_\mathcal{X} g(x) \,\text{d}x$$ $p(x)\propto g(x)$$\alpha(x)$$U(x)$ $$\{x;\alpha(x)>0\}\subset\{x;g(x)>0\}$$ $$\int_\mathcal{X} \dfrac{\alpha(x)}{g(x)}p(x) \,\text{d}x=\int_\mathcal{X} \dfrac{\alpha(x)}{N} \,\text{d}x=\dfrac{1}{N}$$ $p$$1/N$ $$\hat\eta=\frac{1}{n}\sum_{i=1}^n \dfrac{\alpha(x_i)}{g(x_i)}\qquad x_i\stackrel{\text{iid}}{\sim}p(x)$$ $\hat\eta$$\alpha$$\alpha=\pi$ $$\dfrac{\alpha(x)}{g(x)} = \dfrac{1}{\ell(x)}$$ $\ell(x)$$g(x)=\pi(x)\ell(x)$ $$\hat{N}=\dfrac{n}{\sum_{i=1}^n1\big/\ell(x_i)}$$

$(\min(x_i),\max(x_i))$$\exp\{x^2\}$

$\alpha$$(q_{.25}(x_i),q_{.75}(x_i))$$g$

$1/\sqrt{\pi}$

ys = rnorm(1e6, 0, 1/sqrt(2))
r = quantile(ys,.75) - quantile(ys,.25)
yc=ys[(ys>quantile(ys,.25))&(ys<quantile(ys,.75))]
sum(sapply(yc, function(x) 1/( r * exp(-x^2))))/length(ys)
## evaluates to 0.5649015. 1/sqrt(pi) = 0.5641896

Мы подробно обсуждаем этот метод в двух статьях с Дарреном Рейтом и с Жаном-Мишелем Марином .