5

我想找到给定区间内标准正态分布的平均值。

例如，如果我将标准正态分布分成两个 ([-Inf:0] [0:Inf]) 我想得到每一半的平均值。

以下代码几乎完全符合我的要求：

divide <- 2
boundaries <- qnorm(seq(0,1,length.out=divide+1))
t <- sort(rnorm(100000))
means.1 <- rep(NA,divide)
for (i in 1:divide) {
    means.1[i] <- mean(t[(t>boundaries[i])&(t<boundaries[i+1])])
  }

但我需要一种更精确（和优雅）的方法来计算这些数字（means.1）。

我尝试了以下代码，但它不起作用（可能是因为我缺乏概率知识）。

divide <- 2
boundaries <- qnorm(seq(0,1,length.out=divide+1))
means.2 <- rep(NA,divide)
f <- function(x) {x*dnorm(x)}
for (i in 1:divide) {
  means.2[i] <- integrate(f,lower=boundaries[i],upper=boundaries[i+1])$value
}

有任何想法吗？提前致谢。

4

5 回答 5

3

问题是 dnorm(x) 在区间（-Inf 到 0）中的积分不是 1，这就是你得到错误答案的原因。要更正，您必须将得到的结果除以 0.5（积分结果）。像：

func <- function(x, ...) x * dnorm(x, ...)
integrate(func, -Inf, 0, mean=0, sd=1)$value / (pnorm(0, mean=0, sd=1) - pnorm(-Inf, mean=0, sd=1))

适应不同的间隔应该很容易。

于 2013-04-12T18:03:33.980 回答
3

谢谢回答我的问题。

据我了解，我结合了所有答案：

    divide <- 5
    boundaries <- qnorm(seq(0,1,length.out=divide+1))
# My original thinking        
    t <- sort(rnorm(1e6))
    means.1 <- rep(NA,divide)
    for (i in 1:divide) {
        means.1[i] <- mean(t[((t>boundaries[i])&(t<boundaries[i+1]))])
      }    

# Based on @DWin
    t <- sort(rnorm(1e6))
    means.2 <- tapply(t, findInterval(t, boundaries), mean)

# Based on @Rcoster
    means.3 <- rep(NA,divide)
    f <- function(x, ...) x * dnorm(x, ...)
    for (i in 1:divide) {
      means.3[i] <- integrate(f, boundaries[i], boundaries[i+1])$value / (pnorm(boundaries[i+1]) - pnorm(boundaries[i]))
    }   

# Based on @Kith
    t <- sort(rnorm(1e6))
    means.4 <- rep(NA,divide)    
    for (i in 1:divide) {
      means.4[i] <- fitdistr(t[t > boundaries[i] & t < boundaries[i+1]], densfun="normal")$estimate[1]
    }

结果

>   means.1
[1] -1.4004895486 -0.5323784986 -0.0002590746  0.5313539906  1.3978177100
>   means.2   
[1] -1.3993590768 -0.5329465789 -0.0002875593  0.5321381745  1.3990997391 
>   means.3
[1] -1.399810e+00 -5.319031e-01  1.389222e-16  5.319031e-01  1.399810e+00
>   means.4
[1] -1.399057073 -0.531946615 -0.000250952  0.531615180  1.400086731

我相信@Rcoster 是我想要的。与我的相比，休息是创新的方法，但仍然是近似的。谢谢。

于 2013-04-12T20:29:04.253 回答
2

您可以结合使用 fitdistr 和向量索引。

这是一个如何获取正值的均值和标准差的示例：

library("MASS")
x = rnorm(10000)
fitdistr(x[x > 0], densfun="normal")

或者只是区间 (0,2) 中的值：

fitdistr(x[x > 0 & x < 2], densfun="normal")
于 2013-04-12T17:56:25.047 回答
2

假设您的切点是 -1、0、1 和 2，并且您对模拟标准法线的截面的平均值感兴趣。

 samp <-   rnorm(1e5)
 (res <- tapply(samp, findInterval(samp, c( -1, 0, 1, 2)), mean) )
#         0          1          2          3          4 
#-1.5164151 -0.4585519  0.4608587  1.3836470  2.3824633

请注意，标签可以改进。一项改进可能是：

names(res) <-  paste("[", c(-Inf, -1, 0, 1, 2, Inf)[-6],  " , ", 
                      c(-Inf, -1, 0, 1, 2, Inf)[-1], ")", sep="")
> res
[-Inf , -1)    [-1 , 0)     [0 , 1)     [1 , 2)   [2 , Inf) 
 -1.5278185  -0.4623743   0.4621885   1.3834442   2.3835116
于 2013-04-12T18:00:18.550 回答
1

使用distrEx和distr包：

library(distrEx)
E(Truncate(Norm(mean=0, sd=1), lower=0, upper=Inf))
# [1] 0.797884

（有关vignette(distr)distr套件和相关包的出色概述，请参阅distrDoc包。）

lb或者，仅使用基数 R，这是一个替代方案，它在和之间的区间内构造期望的离散近似值ub。调整近似矩形的底边，使它们都具有相等的面积（即，一个点落在每个矩形中的概率相同）。

intervalMean <- function(lb, ub, n=1e5, ...) {
    ## Get x-values at n evenly-spaced quantiles between lower and upper bounds
    xx <- qnorm(seq(pnorm(lb, ...), pnorm(ub, ...), length = n), ...)
    ## Calculate expectation
    mean(xx[is.finite(xx)])
}

## Your example
intervalMean(lb=0, ub=1)
# [1] 0.4598626

## The mean of the complete normal distribution
intervalMean(-Inf, Inf)
## [1] -6.141351e-17

## Right half of standard normal distribution
intervalMean(lb=0, ub=Inf)
# [1] 0.7978606

## Right half of normal distribution with mean 0 and standard deviation 100
intervalMean(lb=0, ub=Inf, mean=0, sd=100)
# [1] 79.78606
于 2013-04-12T18:05:48.470 回答