玩具RコードIは、降雪観測の数を有する

Question

[降雪量のデータが】：

x <- c(98.044, 107.696, 146.050, 102.870, 131.318, 170.434, 84.836, 154.686, 
     162.814, 101.854, 103.378, 16.256) 

と私はそれに知られている標準偏差の正規分布に従うと言われました25.4しかし未知の平均mu。ベイジアンフォーミュラを使ってmuの推論をしなければなりません。

これは、次のmu

の前

mean of snow | 50.8 | 76.2 | 101.6 | 127.0 | 152.4 | 177.8 
--------------------------------------------------------------- 
probability | 0.1 | 0.15 | 0.25 |0.25 | 0.15 | 0.1 
--------------------------------------------------------------- 

の情報では、私がこれまで試してみましたものですが、postについての最終ラインが正常に動作しません。結果のプロットは水平線を与えます。

library(LearnBayes) 
midpts <- c(seq(50.8, 177.8, 30)) 
prob <- c(0.1, 0.15, 0.25, 0.25, 0.15, 0.1) 
p <- seq(50, 180, length = 40000) 
histp <- histprior(p, midpts, prob) 
plot(p, histp, type = "l") 

# posterior density 
post <- round(histp * dnorm(x, 115, 42)/sum(histp * dnorm(x, 115, 42)), 3) 
plot(p, post, type = "l") 

Answer 1

私の最初の提案は、この背後にある統計情報を理解していることを確認してください。私があなたを見たとき

post <- round(histp * dnorm(x, 115, 42)/sum(histp * dnorm(x, 115, 42)), 3) 

私はあなたがいくつかのコンセプトを混乱させたと思った。これはベイズフォーミュラのようですが、あなたはその可能性について誤ったコードを持っています。正しい尤度関数は

## likelihood function: `L(obs | mu)` 
## standard error is known (to make problem easy) at 25.4 
Lik <- function (obs, mu) prod(dnorm(obs, mu, 25.4)) 

注、muが未知であるので、この関数の変数であるべきです。また、尤度は、観測におけるすべての個々の確率密度の積である。今、我々は成功したRの実装について

Lik(x, 100) 
# [1] 6.884842e-30 

によってmu = 100で、例えば可能性を評価することができ、我々は機能Likのためのベクトル化バージョンが必要です。つまり、スカラー入力ではなく、muのベクトル入力で評価できる関数です。私はベクトル化のためのsapplyを使用します。

vecLik <- function (obs, mu) sapply(mu, Lik, obs = obs) 

のは

vecLik(x, c(80, 90, 100)) 
# [1] 6.248416e-34 1.662366e-31 6.884842e-30 

今ではmuための事前分布を得るための時間ですを試してみましょう。原理的にはこれは連続関数ですが、RパッケージLearnBayesのhistpriorを使用して、離散近似が必要なように見えます。ベイズの式を適用する前に

## prior distribution for `mu`: `prior(mu)` 
midpts <- c(seq(50.8, 177.8, 30)) 
prob <- c(0.1, 0.15, 0.25, 0.25, 0.15, 0.1) 
mu_grid <- seq(50, 180, length = 40000) ## a grid of `mu` for discretization 
library(LearnBayes) 
prior_mu_grid <- histprior(mu_grid, midpts, prob) ## discrete prior density 
plot(mu_grid, prior_mu_grid, type = "l") 

、我々は最初の分母の正規化定数NCをうまく。これは、Lik(obs | mu) * prior(mu)の整数になります。しかし、prior(mu)の離散近似があるので、この積分を近似するためにリーマン和を使用します。

delta <- mu_grid[2] - mu_grid[1] ## division size 
NC <- sum(vecLik(x, mu_grid) * prior_mu_grid * delta) ## Riemann sum 
# [1] 2.573673e-28 

グレート、すべての準備ができている、我々はベイズ式を使用することができますprior(mu)が離散化されると、再び

posterior(mu | obs) = Lik(obs | mu) * prior(mu)/NC 

を、posterior(mu)はあまりにも、離散化されます。

post_mu <- vecLik(x, mu_grid) * prior_mu_grid/NC 

ハハ、推論結果を参照するのmuのスケッチ後部を聞かせて：

plot(mu_grid, post_mu, type = "l") 

うわー、これは美しいです！