2016-10-22 4 views
4

[폭설 량 데이터]장난감 R 부호 I는 눈 관측의 개수가

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) 

및 I는 알려진 표준 편차 정규 분포를 따른다 들었다 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") 

답변

15

나의 첫 번째 제안은이 배경의 통계를 이해했는지 확인하는 것입니다. 나는 당신의

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

을보고 난 당신이 몇 가지 개념을 엉망 유망주. 이것은 베이 즈 공식 인 것으로 보이지만 가능성에 대한 잘못된 코드가 있습니다. 정확한 우도 함수는이 함수의 변수되어야하므로 mu은 미지이다

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

주이고; 또한, 확률은 관측치에서 모든 개별 확률 밀도의 산물이다. 이제, 우리는 성공적인 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 패키지 LearnBayeshistprior을 사용하여 불연속 근사를 원한 것처럼 보입니다.

prior distribution

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

Baye의 공식을 적용하기 전에, 우리는 먼저 분모에 정규화 상수 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) 너무, 이산화된다.

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

posterior density

와우, 이것은 아름답다! : mu

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

이 하하,하자 스케치 후방은 추론 결과를 볼 수 있습니다