Trouver la moyenne d'une distribution normale standard dans un intervalle donné

Question

Je veux trouver la moyenne de la distribution normale standard dans un intervalle donné.

Par exemple, si je divise la distribution normale standard en deux ([-Inf:0] [0:Inf]), je veux obtenir la moyenne de chaque moitié.

Le code suivant fait presque exactement ce que je veux :

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])])
  }    

Mais j'ai besoin d'une méthode plus précise (et élégante) pour calculer ces chiffres (moyens.1).

J'ai essayé le code suivant mais cela n'a pas fonctionné (peut-être à cause du manque de mes connaissances en probabilité).

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
}    

Avez-vous des idées ? Merci d'avance.

Answer 1

Le problème est que l'intégrale de dnorm(x) dans l'intervalle (-Inf à 0) n'est pas 1, c'est pourquoi vous avez obtenu la mauvaise réponse. Pour corriger, vous devez diviser le résultat que vous avez obtenu par 0,5 (le résultat de l'intégrale). Comme :

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)) 

Il devrait être facile de l'adapter à différents intervalles.

Answer 2

Merci d'avoir répondu à ma question.

J'ai combiné toutes les réponses comme je le comprends :

    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]
    }    

Résultats

>   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

Je crois que @Rcoster est celui que je voulais. Le reste est des approches innovantes par rapport à la mienne mais toujours approximatives. Merci.

Answer 3

Vous pouvez utiliser une combinaison de fitdistr et d'indexation vectorielle.

Voici un exemple de comment obtenir la moyenne et l'écart-type des valeurs positives uniquement :

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

ou seulement les valeurs dans l'intervalle (0,2) :

fitdistr(x[x > 0 & x < 2], densfun="normal")

Answer 4

Disons que vos points de coupure sont -1, 0, 1 et 2 et que vous vous intéressez à la moyenne des sections simulant une normale standard.

 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 

Veuillez noter que l'étiquetage pourrait être amélioré. Une amélioration serait possible :

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

Answer 5

Utilisation de la distrEx y distr distr paquets :

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

(Voir vignette(distr) dans le distrDoc pour un excellent aperçu de l'ensemble des distr distr et les paquets connexes).

---

Ou, en utilisant seulement le R de base, voici une alternative qui construit une approximation discrète de l'espérance dans l'intervalle entre lb y ub . Les bases des rectangles d'approximation sont ajustées de manière à ce qu'ils aient tous la même surface (c'est-à-dire que la probabilité qu'un point tombe dans chacun d'eux soit identique).

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