Statistics/Bayesian Statistics (23) μΈλ€μΌν 리μ€νΈν νΌν© λͺ¨λΈ(Mixture model) λ€μμ λ°μ΄ν°μ νμ€ν κ·Έλ¨μ νμΈνκ³ , μ΄λ₯Ό μ ν©ν λΆν¬μ λνμ¬ μκ°ν΄ 보μ. In: data = read.csv('../input/mixture.csv', header = F) y = data$V1 n = length(y) hist(y, breaks = 20) Out: β· μΌλ°μ μΈ λΆν¬μλ λ€λ₯΄κ² -2μ 1 κ·Όμ²μ λ κ°μ λ΄μ°λ¦¬λ₯Ό νμ±νκ³ μλ κ²μ νμΈν μ μλ€. μ°λ¦¬κ° μκ³ μλ μ κ·λΆν¬, μ§μλΆν¬, κ°λ§λΆν¬λ λ¨λ΄ννλ₯Ό λκ³ μλ€. λ°λΌμ μ΄λ₯Ό μμ λ°μ΄ν°μ μ ν©ν κ²½μ°, λ κ°μ λ΄μ°λ¦¬μ λν΄ μ ν©ν λΆν¬λ₯Ό μ»μ μ μλ€. βΆ νΌν© λͺ¨λΈ(Mixture model)μ λ κ°μ΄μμ λΆν¬λ₯Ό ν©μ³μ λ§λ λͺ¨λΈλ‘ κΈ°μ‘΄ λΆν¬μ λΉνμ¬ μμ λ λμ μ ν©μ΄ κ°λ₯νλ€λ μ₯μ μ κ°μ§κ³ μλ€. λν 볡μ‘ν λΆ.. μμμ νΈ λͺ¨λΈ(Random intercept model) Rμ car ν¨ν€μ§μ Leinhardt λ°μ΄ν°λ₯Ό μ΄μ©νμ¬ μμμ νΈ λͺ¨λΈ(Random intercept model)μ λͺ¨λΈμ λ§λ€μ΄λ³΄μ. μνκ³Όμ μ λ€μκ³Ό κ°λ€. 1. λ°μ΄ν° νμΈ 2. λͺ¨λΈλ§ 3. λͺ¨λΈνμΈ 1. λ°μ΄ν° νμΈ In: library(car) data('Leinhardt') pairs(Leinhardt) head(Leinhardt) Out: income infant region oil Australia 3426 26.7 Asia no Austria 3350 23.7 Europe no Belgium 3346 17.0 Europe no Canada 4751 16.8 Americas no Denmark 5029 13.5 Europe no Finland 3312 10.1 Europe no β· μ°μν .. Bayesian linear model for New York air quality measurements μΊλ¦¬ν¬λμ λνκ΅μ "Bayesian Statistics: Techniques and Models"μ μ΄μνκΈ° μν νλ‘μ νΈ κ²°κ³Όλ¬Όμ΄λ€. ########################## # setting & loading data # ########################## set.seed(777) library(dplyr) library(tidyr) library(ggplot2) library(GGally) library(rjags) theme_set(theme_light() + theme(plot.title = element_text(face = 'bold', colour = 'grey10'), plot.subtitle = element_text(colour = 'grey25'), panel... κ³μΈ΅μ λͺ¨λΈ(Hierarchical model) λ€μμ μμ λ₯Ό ν΅ν΄ κ³μΈ΅μ λͺ¨λΈ(Hierarchical model)μ νΉμ§μ μμλ³΄κ³ , λͺ¨λΈλ§ κ²°κ³Όμ λν΄ λΆμν΄ λ³΄μ. λ¬Έμ ) μΉμ΄μ μμ°νλ 5κ°μ 곡μ₯μ΄ μλ€. κ° κ³΅μ₯μμ μμ°λ μΉμ΄ κ³Όμ 1κ°μ λ°ν μλ μ΄μ½μΉ© κ°μκ° ν¬μμ‘ λΆν¬λ₯Ό λ°λ₯΄κ³ , ν¬μμ‘ λΆν¬μ λͺ¨μλ κ°λ§λΆν¬λ₯Ό λ°λ₯Έλ€. cookies λ°μ΄ν°λ₯Ό μ΄μ©νμ¬ μΉμ΄ κ³Όμκ° μμ°λ λ, λ°ν μλ μ΄μ½μΉ© κ°μμ λν λͺ¨λΈλ§μ μνν ν, λΆμνμμ€. νμ΄) β· μμ λ¬Έμ μ λν΄ ν¬κ² 3κ°μ§ λͺ¨λΈλ§ λ°©λ²μΌλ‘ μ κ·Όμ΄ κ°λ₯νλ€. (1) Fully independent model: λͺ¨λ λ°μ΄ν°κ° λ 립μ΄λΌ κ°μ νκ³ , νλμ ν¬μμ‘ λͺ¨λΈμ λ§λλ κ²μ΄λ€. μ΄λ κ° κ³΅μ₯λ³ μ°¨μ΄μ κ°μ 곡μ₯μμ μμ°λ μΉμ΄μ λΉμ·ν νΉμ±μ κ³ λ €νμ§ λͺ»νλ€λ νκ³κ° μλ€. .. λ² μ΄μ§μ ν¬μμ‘ νκ·(Bayesian poisson regression) λ¨μ ν¬μμ‘ νκ·(Simple poisson regression) λͺ¨λΈμ λ€μκ³Ό κ°λ€. β· κ°λ₯λλ ν¬μμ‘ λΆν¬λ‘ μ νκ³ , κ°λ₯λ λͺ¨μμ λ‘κ·Έλ₯Ό μ·¨ν κ²μ λνμ¬ λ 립λ³μμ μ νκ²°ν©μΌλ‘ μ μνλ€. μ΄λ, μμ μμμλ νλμ λ 립λ³μμ λν μ νκ²°ν©μΌλ‘ νννμμ§λ§, λ 립λ³μκ° μ¬λ¬κ°μ΄λ©΄ μ΄ λ 립λ³μλ€μ μ νκ²°ν©μ ν΅ν΄ κ°λ₯λ λͺ¨μμ λ‘κ·Έλ₯Ό μ·¨ν κ°μ λνμ¬ μ μνλ€. β· ν¬μμ‘ νκ· λͺ¨λΈμ μμΈ‘μ κ°λ³ λ°μ΄ν°μ κ°λ₯λμ νκ· , μ¦, κ°λ₯λμ λͺ¨μλ₯Ό ν΅ν΄ μ΄λ£¨μ΄μ§λ€. β· λ² μ΄μ§μ ν¬μμ‘ νκ·(Bayesian poisson regression)λ λ 립λ³μμ μ νκ²°ν©λ λͺ¨μ betaμ μ¬μ λΆν¬λ₯Ό μ μνλ€λ μ μμ ν¬μμ‘ νκ·μ μ°¨μ΄κ° μλ€. Rμ COUNT ν¨ν€μ§μ badhealth λ°μ΄ν°λ₯Ό μ΄μ©νμ¬ λ² μ΄μ§.. λ² μ΄μ§μ λ‘μ§μ€ν± νκ·(Bayesian logistic regression) λ¨μ λ‘μ§μ€ν± νκ·(Simple logistic regression) λͺ¨λΈμ λ€μκ³Ό κ°λ€. β· κ°λ₯λλ λ² λ₯΄λμ΄ λΆν¬λ‘ μ νκ³ , κ°λ₯λ λͺ¨μμ λ‘μ§(Logit)μ λ 립λ³μμ μ νκ²°ν©μΌλ‘ μ μνλ€. μ΄λ, μμ μμμλ νλμ λ 립λ³μμ λν μ νκ²°ν©μΌλ‘ νννμμ§λ§, λ 립λ³μκ° μ¬λ¬κ°μ΄λ©΄ μ΄ λ 립λ³μλ€μ μ νκ²°ν©μ ν΅ν΄ λ‘μ§μ μ μνλ€. β· λ‘μ§μ€ν± νκ·μ μμΈ‘μ κ°λ³ λ°μ΄ν°μ κ°λ₯λμ νκ· , μ¦ κ°λ₯λμ λͺ¨μλ₯Ό ν΅ν΄ μ΄λ£¨μ΄μ§λ€. βΆ λ² μ΄μ§μ λ‘μ§μ€ν± νκ·(Bayesian logistic regression)λ λ‘μ§μ λνλ΄λ λ 립λ³μμ μ νκ²°ν©λ λͺ¨μμ μ¬μ λΆν¬λ₯Ό μ μνλ€λ μ μμ λ‘μ§μ€ν± νκ·μ μ°¨μ΄κ° μλ€. Rμ boot ν¨ν€μ§μ urine λ°μ΄ν°λ₯Ό μ΄μ©νμ¬ λ² μ΄μ§μ λ‘μ§μ€ν± νκ· λΆμμ μνν κ².. DIC(Deviance Information Criterion) λ² μ΄μ§μ λͺ¨λΈμμλ λͺ¨λΈ μ νμ μν μ 보μ κΈ°μ€μΌλ‘μ¨ DIC(Deviance Information Criterion)μ μ μνκ³ μλ€. DICμ 곡μμ λ€μκ³Ό κ°λ€. β· theta hatμ κ° λͺ¨μμ μ¬ννκ· μ΄κ³ , μ¬νλΆν¬λ‘λΆν° μ»μ theta hatμ λ‘κ·Έ κ°λ₯λμ μ€μ§μ μΈ λͺ¨μμ κ°―μ(Effective number of parameters)λ₯Ό κ³ λ €νμ¬ DICλ₯Ό ꡬν μ μλ€. βΆ μ€μ§μ μΈ λͺ¨μμ κ°―μλ λͺ¨λΈμ μΆμ μΉ μ¬μ΄μ μκ΄(Correlation)μ κ³ λ €νκΈ° μν κ²μ΄λ€. μλ₯Ό λ€μ΄, λͺ¨λΈμ μΆμ μΉ μ¬μ΄μ 0.99μ μκ΄μ΄ μ‘΄μ¬νλ€λ©΄ μ΄λ₯Ό λ 립μ μΈ λͺ¨μλ‘ κ°μ£Όνλ€λ©΄ ν©λ¦¬μ μ΄μ§ μμ κ²μ΄λ€. Rμ car ν¨ν€μ§μ Leihardt λ°μ΄ν°λ₯Ό μ΄μ©ν λ λͺ¨λΈμ DICλ₯Ό ν΅ν΄ λΉκ΅νμ¬λ³΄μ. In: lib.. λ² μ΄μ§μ μ ν νκ·(Bayesian linear regression) Rμ car ν¨ν€μ§μ Leihardt λ°μ΄ν°λ₯Ό μ΄μ©νμ¬ λ² μ΄μ§μ μ ν νκ·(Bayesian linear regression) λΆμμ μνν κ²μ΄λ€. μν κ³Όμ μ λ€μκ³Ό κ°λ€. 1. λ°μ΄ν° νμΈ 2. λͺ¨λΈλ§ 3. λͺ¨λΈ νμΈ 4. μμ°¨ λΆμ 1. λ°μ΄ν° νμΈ In: library(car) data('Leinhardt') pairs(Leinhardt) Out: β· Leihardt λ°μ΄ν°λ μ°μν λ³μμΈ income, infant, regionκ³Ό λ²μ£Όν λ³μμΈ oilλ‘ κ΅¬μ±λμ΄ μλ€. β· incomeκ³Ό infantκ° λΉμ νμ μΈ κ΄κ³λ₯Ό λνλ΄κ³ μλ κ²μ νμΈν μ μλ€. μ νλͺ¨λΈμ κΈ°λ³Έμ μΌλ‘ λ³μκ° μ νκ΄κ³λ₯Ό κ°μ νκΈ° λλ¬Έμ λ λ³μμ λ‘κ·Έλ₯Ό μ·¨νμ¬ λ³νν κ²μ΄λ€. In: Leinhardt$log_income .. MCMC(Markov Chain Monte-Carlo)μ μλ ΄(Convergence) MCMC(Markov Chain Monte-Carlo)λ₯Ό ν΅ν΄ μμ±ν λ°μ΄ν°λ₯Ό νμ©νκΈ° μν΄μλ λ§λ₯΄μ½ν 체μΈ(Markov Chain)μ΄ μ μμν(Stationary)μ μλ ΄(Convergence)ν΄μΌ νλ€. μ΄λ₯Ό νμΈνκ³ λ€λ£¨λ λ°©λ²μ λν΄ λ€λ£° κ²μ΄λ€. 1. Trace plot 2. μκΈ°μκ΄μ±(Autocorrelation) 3. μ΄κΈ° λ¨κ³(Burn-in period) 1. Trace plot MCMCμ μλ ΄μ νμΈνλ κ°μ₯ μ§κ΄μ μΈ λ°©λ²μ λ°μ΄ν°μ μμ± κ³Όμ μ μ§μ κ·Έλ¦ΌμΌλ‘ λνλ΄λ κ²μ΄λ€. μννμμ λ°λ₯Έ μμ±λ λ°μ΄ν°μ λΆν¬λ₯Ό ν΅ν΄ μ΄λ₯Ό νμΈν μ μλ€. In: log_g κΉμ€ μνλ§(Gibbs sampling) κΉμ€ μνλ§(Gibbs sampling)μ λν΄ μμλ³Ό κ²μ΄λ€. λ€λ£° λ΄μ©μ λ€μκ³Ό κ°λ€. 1. κΉμ€ μνλ§ 2. κΉμ€ μνλ§μ μμ 1. κΉμ€ μνλ§ κΉμ€ μνλ§μ Metropolis Hastings(μ΄ν MH) μκ³ λ¦¬μ¦μ νΉλ³ν ννλ‘, μ μ λΆν¬(Proposal distribution)λ₯Ό μμ μ Full conditional distributionλ‘ λμ΄ μνλ§νλ λ°©λ²μ΄λ€. μ΄λ κ² ν¨μΌλ‘μ¨, κ° μνμμ λ°μνλ λ°μ΄ν°μ λν΄ Acceptance probabilityλ 1μ΄ λλ μ±μ§μ κ°μ§κ² λλ€. λ€μμ μ¦λͺ μ ν΅νμ¬ μ΄λ₯Ό νμΈν΄λ³΄μ. β· μ μ λΆν¬λ₯Ό Full conditional posteriorλ‘ λ μΌλ‘μ¨ λ―ΈμΈ κ· ν 쑰건(Detailed balance condition)μ΄ μ±λ¦½νκ² λλ€. .. JAGS(Just Another Gibbs Sampler) μ¬μ©λ² Rμ JAGS(Just Another Gibbs Sampler)μ μ¬μ©λ²μ λν΄ μμλ³Ό κ²μ΄λ€. JAGSλ₯Ό ν΅ν λ°μ΄ν° μμ± κ³Όμ μ 4λ¨κ³λ‘ λλ μ μλ€. 1. Specify the model 2. Set up the model 3. Run the MCMC(Markov Chain Monte Carlo) sampler 4. Post processing λ€μμ λͺ¨λΈμ μ΄λ₯Ό λ¨κ³λ³λ‘ μ μ©νμ¬ μ¬νλΆν¬λ‘λΆν° λ°μ΄ν°λ₯Ό μμ±νμ¬ λ³΄μ. 1. Specify the model In: library(rjags) mod_string = " model { for (i in 1:n) { y[i] ~ dnorm(mu, 1.0/sig2) } mu ~ dt(0.0, 1.0/1.0, 1) sig2 = 1.0 } " β· μμ μ½λμ.. λ©νΈλ‘ν΄λ¦¬μ€ ν€μ΄μ€ν μ€ μκ³ λ¦¬μ¦(Metropolis-Hastings algorithm) Metropolis-Hastings(μ΄ν MH) μκ³ λ¦¬μ¦μ λν΄ μμλ³Ό κ²μ΄λ€. MH μκ³ λ¦¬μ¦μ MCMC(Markov Chain Monte-Carlo)μ μΌλ°μ μΈ ννλ‘μ¨ νΉμ λΆν¬λ‘λΆν° μ μλΆν¬λ‘ κ°λ 체μΈμ λ°μμν¬ μ μλ€. μ΄λ₯Ό μ΄μ©νμ¬ νΉμ λΆν¬λ‘λΆν° λ°μ΄ν°λ₯Ό μμ±ν μ μλ€. λ€λ£° λ΄μ©μΌλ‘λ λ€μκ³Ό κ°λ€. 1. MH μκ³ λ¦¬μ¦ 2. Random walk MH μκ³ λ¦¬μ¦ κ΅¬ν 1. MH μκ³ λ¦¬μ¦ MH μκ³ λ¦¬μ¦μ λ€μκ³Ό κ°λ€. β· qλ μ μ λΆν¬(Proposal distribution)λ₯Ό μλ―Ένκ³ , gλ μ°λ¦¬μ λͺ©μ λΆν¬(Target distribution)μμ μ κ·ν μμ(Normalizing constant)λ₯Ό μ μΈν λΆλΆμ΄λ€. μ¦, λͺ©μ λΆν¬μ g(theta)λ λΉλ‘ κ΄κ³κ° μ±λ¦½νλ€. β· μ΄κΈ°κ°.. μ΄μ 1 2 λ€μ