与えられたディストリビューションの区間確率を見つける方法は？

Question

私はいくつかのデータを持っており、gammaディストリビューションに合うと仮定します。Pr(1 < x <= 1.5)の区間確率を求めるにはどうすればいいですか？xはサンプル外のデータポイントですか？

require(fitdistrplus) 

a <- c(2.44121289,1.70292449,0.30550832,0.04332383,1.0553436,0.26912546,0.43590885,0.84514809, 
0.36762336,0.94935435,1.30887437,1.08761895,0.66581035,0.83108270,1.7567334,1.00241339, 
0.96263021,1.67488277,0.87400413,0.34639636,1.16804671,1.4182144,1.7378907,1.7462686, 
1.7427784,0.8377457,0.1428738,0.71473956,0.8458882,0.2140742,0.9663167,0.7933085, 
0.0475603,1.8657773,0.18307362,1.13519144) 

fit <- fitdist(a, "gamma",lower = c(0, 0)) 

Answer 1

誰かが私の上記のアプローチを好まない。これはMLEの条件付きである。無条件に何かを見てみましょう。直接統合する場合は、shapeの場合は1つ、rateの場合は1つ、xの場合は3つの統合が必要です。これは魅力的ではありません。代わりにMonte Carloの見積もりを作成します。

セントラルリミット定理では、MLEは通常分布しています。 fitdistrplus::fitdistは標準エラーを出さないが、ここでは正確な推論を行うMASS::fitdistrを使うことができる。

fit <- fitdistr(a, "gamma", lower = c(0,0)) 

b <- fit$estimate 
# shape  rate 
#1.739737 1.816134 

V <- fit$vcov ## covariance 
      shape  rate 
shape 0.1423679 0.1486193 
rate 0.1486193 0.2078086 

ここでは、パラメータ分布からサンプルし、ターゲット確率のサンプルを取得したいと考えています。

set.seed(0) 
## sample from bivariate normal with mean `b` and covariance `V` 
## Cholesky method is used here 
X <- matrix(rnorm(1000 * 2), 1000) ## 1000 `N(0, 1)` normal samples 
R <- chol(V) ## upper triangular Cholesky factor of `V` 
X <- X %*% R ## transform X under desired covariance 
X <- X + b ## shift to desired mean 
## you can use `cov(X)` to check it is very close to `V` 

## now samples for `Pr(1 < x < 1.5)` 
p <- pgamma(1.5, X[,1], X[,2]) - pgamma(1, X[,1], X[,2]) 

我々はpのヒストグラムを作る（そして、あなたがしたい場合は、多分密度推定を行う）ことができます。

hist(p, prob = TRUE) 

今、私たちは多くの場合、サンプルは、予測のために意味します：

mean(p) 
# [1] 0.1906975 

Answer 2

あなただけのfitに推定されたパラメータでpgammaを使用することができます。

b <- fit$estimate 
# shape  rate 
#1.739679 1.815995 

pgamma(1.5, b[1], b[2]) - pgamma(1, b[1], b[2]) 
# [1] 0.1896032 

---

感謝。しかし、どうすればP(x > 2)？デフォルトでは

pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, log.p = FALSE) 

、pgamma(q)はPr(x <= q)を評価：

はlower.tail引数をチェックしてください。 lower.tail = FALSEを設定すると、Pr(x > q)となります。だから、行うことができます。

pgamma(2, b[1], b[2], lower.tail = FALSE) 
# [1] 0.08935687 

それとも、またここで

1 - pgamma(2, b[1], b[2]) 
# [1] 0.08935687 

Answer 3

を使用することができ、新たな観測がに落ちることを事後確率を推定するために、MCMC技術を使用する例と推論のベイズモードを行くをインターバル（1：1.5）。これは、ガンマ分布を最尤パラメータ推定値と積分することによって得られる条件付き推定値とは対照的に、無条件推定値である。

このコードでは、JAGSをコンピュータにインストールする必要があります（無料で簡単にインストールできます）。

library(rjags) 

a <- c(2.44121289,1.70292449,0.30550832,0.04332383,1.0553436,0.26912546,0.43590885,0.84514809, 
     0.36762336,0.94935435,1.30887437,1.08761895,0.66581035,0.83108270,1.7567334,1.00241339, 
     0.96263021,1.67488277,0.87400413,0.34639636,1.16804671,1.4182144,1.7378907,1.7462686, 
     1.7427784,0.8377457,0.1428738,0.71473956,0.8458882,0.2140742,0.9663167,0.7933085, 
     0.0475603,1.8657773,0.18307362,1.13519144) 

# Specify the model in JAGS language using diffuse priors for shape and scale 
sink("GammaModel.txt") 
cat("model{ 

    # Priors 
    shape ~ dgamma(.001,.001) 
    rate ~ dgamma(.001,.001) 

    # Model structure 
    for(i in 1:n){ 
    a[i] ~ dgamma(shape, rate) 
    } 
    } 
    ", fill=TRUE) 
sink() 

jags.data <- list(a=a, n=length(a)) 

# Give overdispersed initial values (not important for this simple model, but very important if running complicated models where you need to check convergence by monitoring multiple chains) 
inits <- function(){list(shape=runif(1,0,10), rate=runif(1,0,10))} 

# Specify which parameters to monitor 
params <- c("shape", "rate") 

# Set-up for MCMC run 
nc <- 1 # number of chains 
n.adapt <-1000 # number of adaptation steps 
n.burn <- 1000 # number of burn-in steps 
n.iter <- 500000 # number of posterior samples 
thin <- 10 # thinning of posterior samples 

# Running the model 
gamma_mod <- jags.model('GammaModel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt) 
update(gamma_mod, n.burn) 
gamma_samples <- coda.samples(gamma_mod,params,n.iter=n.iter, thin=thin) 
# Summarize the result 
summary(gamma_samples) 

# Compute improper (non-normalized) probability distribution for x 
x <- rep(NA, 50000) 
for(i in 1:50000){ 
    x[i] <- rgamma(1, gamma_samples[[1]][i,1], rate = gamma_samples[[1]][i,2]) 
} 

# Find which values of x fall in the desired range and normalize. 
length(which(x>1 & x < 1.5))/length(x) 

回答： Pr(1 < x <= 1.5) = 0.194 とてもきれい条件付き推定値に近いが、これは一般的ケースであることが保証されていません。