Comments (4)
I guess you mean code block 13.12? Without your code, I can't help you much.
In general, the correlation matrix will read its dimension off the dimension of the multi_normal distribution it is embedded in. If it's not embedded in one, then you'll need to manually specify it's dimension.
from rethinking.
Yes, exactly, block 13.12. I've been working through the examples in ch 12 & 13. Here is the code, including the data simulation:
simulate data
N_cafes <- 20
library(MASS)
set.seed(5) # used to replicate samples
vary_effects <- mvrnorm(N_cafes, Mu, Sigma)
a_cafe <- vary_effects[,1]
b_cafe <- vary_effects[,2]
N_visits <- 10
afternoon <- rep(0:1, N_visits*N_cafes/2)
cafe_id <- rep(1:N_cafes, each=N_visits)
mu <- a_cafe[cafe_id] + b_cafe[cafe_id] * afternoon
sigma <- 0.5 # sd within cafes
wait <- rnorm(N_visits * N_cafes, mu, sigma)
d <- data.frame(cafe=cafe_id, afternoon=afternoon, wait=wait)
fit model
m13.1 <- map2stan(
alist(
wait ~ dnorm(mu, Sigma),
mu <- a_cafe[cafe] + b_cafe[cafe] * afternoon,
c(a_cafe, b_cafe) ~ dmvnorm2(c(a,b), sigma_cafe, Rho),
a ~ dnorm(0, 10),
b ~ dnorm(0, 10),
sigma_cafe ~ dcauchy(0, 2),
sigma ~ dcauchy(0, 2),
Rho ~ dlkjcorr(2)
),
data = d, iter = 5000, warmup = 2000, chains = 2
)
from rethinking.
Your code is missing some variable assignments, like Mu
. This is the needed code, starting from the beginning of the chapter:
library(rethinking)
## R code 13.1
a <- 3.5 # average morning wait time
b <- (-1) # average difference afternoon wait time
sigma_a <- 1 # std dev in intercepts
sigma_b <- 0.5 # std dev in slopes
rho <- (-0.7) # correlation between intercepts and slopes
## R code 13.2
Mu <- c( a , b )
## R code 13.3
cov_ab <- sigma_a*sigma_b*rho
Sigma <- matrix( c(sigma_a^2,cov_ab,cov_ab,sigma_b^2) , ncol=2 )
## R code 13.5
sigmas <- c(sigma_a,sigma_b) # standard deviations
Rho <- matrix( c(1,rho,rho,1) , nrow=2 ) # correlation matrix
# now matrix multiply to get covariance matrix
Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas)
## R code 13.6
N_cafes <- 20
## R code 13.7
library(MASS)
set.seed(5) # used to replicate example
vary_effects <- mvrnorm( N_cafes , Mu , Sigma )
## R code 13.8
a_cafe <- vary_effects[,1]
b_cafe <- vary_effects[,2]
## R code 13.10
N_visits <- 10
afternoon <- rep(0:1,N_visits*N_cafes/2)
cafe_id <- rep( 1:N_cafes , each=N_visits )
mu <- a_cafe[cafe_id] + b_cafe[cafe_id]*afternoon
sigma <- 0.5 # std dev within cafes
wait <- rnorm( N_visits*N_cafes , mu , sigma )
d <- data.frame( cafe=cafe_id , afternoon=afternoon , wait=wait )
## R code 13.12
m13.1 <- map2stan(
alist(
wait ~ dnorm( mu , sigma ),
mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
c(a_cafe,b_cafe)[cafe] ~ dmvnorm2(c(a,b),sigma_cafe,Rho),
a ~ dnorm(0,10),
b ~ dnorm(0,10),
sigma_cafe ~ dcauchy(0,2),
sigma ~ dcauchy(0,2),
Rho ~ dlkjcorr(2)
) ,
data=d ,
iter=5000 , warmup=2000 , chains=2 )
If that doesn't work for you, then I suspect there is a gremlin in your global R namespace.
from rethinking.
I forgot to say: Doh! Thank you.
from rethinking.
Related Issues (20)
- Sampling issue for the Cat Adoptions Survival model.
- Random walk in ulam HOT 2
- problem in precis command
- Can multiple outcomes be regressed on the same linear model?
- ulam fits missing @stanfit slot? HOT 2
- WAIC not computing from ulam() model, even with log_lik = TRUE HOT 4
- Problem installing cmdstanr HOT 1
- ulam fits giving odd errors HOT 4
- errata 2nd Edition page 130
- log_lik fails with matrices in models fit by ulam HOT 6
- Where is sim_globe, and the rest of the code shown in the videos? HOT 2
- Panda_nuts Nut-Cracking Frequency? HOT 6
- The garden of forking data - 2 blue from 1? HOT 3
- quap allows incorrect(?) definition with categorical variable HOT 1
- Consistent syntax error on map2stan/ulam issues HOT 2
- Some typos here and there
- rethinking cannot be installed: Rcpp fails to compile HOT 1
- sim() not working for models with simplex
- Saving ulam models won't save cmdstanr csv HOT 3
- Error recovering results from cmdstanr HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from rethinking.