机器算法验证 - 正态截断随机变量的总和 - 吾爱随笔录

正态截断随机变量的总和

机器算法验证 r 正态分布密度函数截断鞍点近似

2022-03-15 19:24:42

假设我有 $n$ 独立正态随机变量

X_{1} \sim N (μ_{1}, σ_{1}^{2}) X_{2} \sim N (μ_{2}, σ_{2}^{2}) ⋮ X_{n} \sim N (μ_{n}, σ_{n}^{2})

$X_1 \sim \mathrm{N}(\mu_1, \sigma_1^2)\\X_2 \sim \mathrm{N}(\mu_2, \sigma_2^2)\\\vdots\\X_n \sim \mathrm{N}(\mu_n, \sigma_n^2)$

和 $Y=X_1+X_2+\dotsm+X_n$ . 我将如何表征密度 $Y$ 如果每个分布 $X_i$ 每个都被截断到内 $(\mu_i - 2\sigma_i, \mu_i + 2\sigma_i)$ ? 换句话说，我从 $n$ 独立正态分布，丢弃不在范围内的样本 $2\sigma_i$ 每个平均值，并将它们相加。

现在，我正在使用下面的 R 代码执行此操作：

x_mu <- c(12, 18, 7)
x_sd <- c(1.5, 2, 0.8)
a <- x_mu - 2 * x_sd
b <- x_mu + 2 * x_sd

samples <- sapply(1:3, function(i) {
  return(rtruncnorm(100000, a[i], b[i], x_mu[i], x_sd[i]))
})

y <- rowSums(samples)

有什么方法可以生成密度 $Y$ 直接地？

2个回答

对于截断法线的总和，您可以使用鞍点方法的近似值。我现在不会给出详细信息，您可以查看我对 Gamma 分布的一般总和的回答以获取提示。我们需要的是找到截断法线的矩生成函数，这很容易。我将在这里为标准法线截断 $\pm 2$ , 有密度

f (x) = {\begin{cases} \frac{1}{C} ϕ (x), & | x | \leq 2 \\ 0, & | x | > 2 \end{cases}

$f(x) =\begin{cases} \frac1{C} \phi(x), & |x| \le 2 \\ 0, & |x| > 2 \end{cases}$ 在哪里

C = Φ (2) - Φ (- 2)

$C=\Phi(2) - \Phi(-2)$ 这里

ϕ (x), Φ (x)

$\phi(x), \Phi(x)$ 分别是标准法线的密度和 cdf。

矩生成函数可以计算为

M (t) = E e^{t X} = \frac{1}{C} \int_{- 2}^{2} e^{t x} ϕ (x) d x = \frac{1}{C} e^{\frac{1}{2} t^{2}} [Φ (2 - t) - Φ (- 2 - t)]

$\DeclareMathOperator{\E}{\mathbb{E}} M(t) = \E e^{tX}=\frac1{C}\int_{-2}^2 e^{tx} \phi(x)\; dx=\frac1{C}e^{\frac12 t^2} [\Phi(2-t)-\Phi(-2-t) ]$ 然后我们可以使用鞍点近似。

我很好奇为什么，但是是的，有一种简单的方法可以生成这个分布总和的 pdf：

## install.packages("truncnorm")
## install.packages("caTools")
library(truncnorm)

x.mu <- c(12, 18, 7)
x.sd <- c(1.5, 2, 0.8)
x.a <- x.mu - 2*x.sd
x.b <- x.mu + 2*x.sd

dmulti <- function(x, a, b, mu, sd)
  rowSums(
    sapply(1:length(mu),
           function(idx)
             dtruncnorm(x, a=a[idx], b=b[idx], mean=mu[idx], sd=sd[idx])))/length(mu)
pmulti <- function(q, a, b, mu, sd)
  rowSums(
    sapply(1:length(mu),
           function(idx)
             ptruncnorm(q, a=a[idx], b=b[idx], mean=mu[idx], sd=sd[idx])))/length(mu)

pointrange <- range(c(x.a, x.b))
pointseq <- seq(pointrange[1], pointrange[2], length.out=100)
## Plot the probability density function
plot(pointseq, dmulti(pointseq, x.a, x.b, x.mu, x.sd),
     type="l")

## Plot the cumulative distribution function
plot(pointseq, pmulti(pointseq, x.a, x.b, x.mu, x.sd),
     type="l")

其它你可能感兴趣的问题

上一篇SVM 模型的概率是如何计算的？下一篇用 3 类结果解释优势比小于 1