library(datasets)
data(attitude)
y <-attitude$rating
N <- length(y)
choose <- sample(seq(1:N),N,replace=FALSE)
holdout <- y[choose[1:10]]
training <- y[choose[11:N]]
library(rstan)
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
dat= list(attitude$complaints, attitude$privileges, attitude$learning, attitude$raises, attitude$critical, attitude$advance, attitude$rating)
com = (attitude$complaints - mean(attitude$complaints))/sd(attitude$complaints)
pri = (attitude$privileges - mean(attitude$privileges))/sd(attitude$privileges)
lea = (attitude$learning - mean(attitude$learning))/sd(attitude$learning)
rai = (attitude$raises - mean(attitude$raises))/sd(attitude$raises)
cri = (attitude$critical - mean(attitude$critical))/sd(attitude$critical)
adv = (attitude$advance - mean(attitude$advance))/sd(attitude$advance)
rat = (attitude$rating - mean(attitude$rating))/sd(attitude$rating)
# Find the least-squares solution
rat.lm <- lm(rat~com + pri + lea + rai + cri + adv,data=trees)
# show results
summary(rat.lm)
##
## Call:
## lm(formula = rat ~ com + pri + lea + rai + cri + adv, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.89889 -0.35781 0.02595 0.45533 0.95288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.717e-16 1.060e-01 0.000 1.000000
## com 6.707e-01 1.761e-01 3.809 0.000903 ***
## pri -7.343e-02 1.364e-01 -0.538 0.595594
## lea 3.089e-01 1.625e-01 1.901 0.069925 .
## rai 6.981e-02 1.892e-01 0.369 0.715480
## cri 3.120e-02 1.195e-01 0.261 0.796334
## adv -1.835e-01 1.506e-01 -1.218 0.235577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5806 on 23 degrees of freedom
## Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
## F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05
model.script = "
data {
int N;
real com[N];
real pri[N];
real lea[N];
real rai[N];
real cri[N];
real adv[N];
real rat[N];
}
parameters {
real beta0;
real beta1;
real beta2;
real beta3;
real beta4;
real beta5;
real ln_sigma;
}
transformed parameters {
real <lower=0> sigma;
sigma=exp(-ln_sigma);
}
model {
real mu[N];
beta0 ~ normal(0,1);
beta1 ~ normal(0,1);
beta2 ~ normal(0,1);
beta3 ~ normal(0,1);
beta4 ~ normal(0,1);
beta5 ~ normal(0,1);
ln_sigma ~ normal(0,1);
for (n in 1:N) {
mu[n] = beta0*com[n] + beta1*pri[n]+beta2*lea[n]+beta3*rai[n]+beta4*cri[n]+beta5*adv[n];
rat[n] ~ normal(mu[n],sigma);
}
}
"
brn <- 250
n.iter <- 1000
fit <- stan(model_code = model.script, data = dat,
iter = n.iter+brn, warmup=brn, chains = 4)
## In file included from C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/config.hpp:39:0,
## from C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file1f10429a5ec6.cpp:8:
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## In file included from C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/base.hpp:28:0,
## from C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array.hpp:21,
## from C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29,
## from C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/numeric/odeint.hpp:61,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/prim/arr.hpp:38,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/prim/mat.hpp:298,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/mat.hpp:12,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file1f10429a5ec6.cpp:8:
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<N>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index_range index_range;
## ^
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index index;
## ^
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<0ull>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index_range index_range;
## ^
## C:/Users/mahr/Documents/R/win-library/3.0/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index index;
## ^
## In file included from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/core.hpp:42:0,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file1f10429a5ec6.cpp:8:
## C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:
## C:/Users/mahr/Documents/R/win-library/3.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
##
## SAMPLING FOR MODEL '42401f2d51e3cddc0095c14c3ba65b11' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 1250 [ 0%] (Warmup)
## Iteration: 125 / 1250 [ 10%] (Warmup)
## Iteration: 250 / 1250 [ 20%] (Warmup)
## Iteration: 251 / 1250 [ 20%] (Sampling)
## Iteration: 375 / 1250 [ 30%] (Sampling)
## Iteration: 500 / 1250 [ 40%] (Sampling)
## Iteration: 625 / 1250 [ 50%] (Sampling)
## Iteration: 750 / 1250 [ 60%] (Sampling)
## Iteration: 875 / 1250 [ 70%] (Sampling)
## Iteration: 1000 / 1250 [ 80%] (Sampling)
## Iteration: 1125 / 1250 [ 90%] (Sampling)
## Iteration: 1250 / 1250 [100%] (Sampling)
##
## Elapsed Time: 0.033 seconds (Warm-up)
## 0.12 seconds (Sampling)
## 0.153 seconds (Total)
##
##
## SAMPLING FOR MODEL '42401f2d51e3cddc0095c14c3ba65b11' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 1250 [ 0%] (Warmup)
## Iteration: 125 / 1250 [ 10%] (Warmup)
## Iteration: 250 / 1250 [ 20%] (Warmup)
## Iteration: 251 / 1250 [ 20%] (Sampling)
## Iteration: 375 / 1250 [ 30%] (Sampling)
## Iteration: 500 / 1250 [ 40%] (Sampling)
## Iteration: 625 / 1250 [ 50%] (Sampling)
## Iteration: 750 / 1250 [ 60%] (Sampling)
## Iteration: 875 / 1250 [ 70%] (Sampling)
## Iteration: 1000 / 1250 [ 80%] (Sampling)
## Iteration: 1125 / 1250 [ 90%] (Sampling)
## Iteration: 1250 / 1250 [100%] (Sampling)
##
## Elapsed Time: 0.03 seconds (Warm-up)
## 0.109 seconds (Sampling)
## 0.139 seconds (Total)
##
##
## SAMPLING FOR MODEL '42401f2d51e3cddc0095c14c3ba65b11' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 1250 [ 0%] (Warmup)
## Iteration: 125 / 1250 [ 10%] (Warmup)
## Iteration: 250 / 1250 [ 20%] (Warmup)
## Iteration: 251 / 1250 [ 20%] (Sampling)
## Iteration: 375 / 1250 [ 30%] (Sampling)
## Iteration: 500 / 1250 [ 40%] (Sampling)
## Iteration: 625 / 1250 [ 50%] (Sampling)
## Iteration: 750 / 1250 [ 60%] (Sampling)
## Iteration: 875 / 1250 [ 70%] (Sampling)
## Iteration: 1000 / 1250 [ 80%] (Sampling)
## Iteration: 1125 / 1250 [ 90%] (Sampling)
## Iteration: 1250 / 1250 [100%] (Sampling)
##
## Elapsed Time: 0.031 seconds (Warm-up)
## 0.116 seconds (Sampling)
## 0.147 seconds (Total)
##
##
## SAMPLING FOR MODEL '42401f2d51e3cddc0095c14c3ba65b11' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 1250 [ 0%] (Warmup)
## Iteration: 125 / 1250 [ 10%] (Warmup)
## Iteration: 250 / 1250 [ 20%] (Warmup)
## Iteration: 251 / 1250 [ 20%] (Sampling)
## Iteration: 375 / 1250 [ 30%] (Sampling)
## Iteration: 500 / 1250 [ 40%] (Sampling)
## Iteration: 625 / 1250 [ 50%] (Sampling)
## Iteration: 750 / 1250 [ 60%] (Sampling)
## Iteration: 875 / 1250 [ 70%] (Sampling)
## Iteration: 1000 / 1250 [ 80%] (Sampling)
## Iteration: 1125 / 1250 [ 90%] (Sampling)
## Iteration: 1250 / 1250 [100%] (Sampling)
##
## Elapsed Time: 0.03 seconds (Warm-up)
## 0.122 seconds (Sampling)
## 0.152 seconds (Total)
print(fit)
## Inference for Stan model: 42401f2d51e3cddc0095c14c3ba65b11.
## 4 chains, each with iter=1250; warmup=250; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta0 0.66 0.00 0.18 0.31 0.54 0.66 0.78 1.00 2618 1
## beta1 -0.07 0.00 0.14 -0.34 -0.16 -0.07 0.02 0.21 3363 1
## beta2 0.31 0.00 0.16 -0.02 0.20 0.31 0.42 0.62 3152 1
## beta3 0.07 0.00 0.20 -0.32 -0.06 0.07 0.20 0.46 2632 1
## beta4 0.03 0.00 0.12 -0.21 -0.05 0.03 0.11 0.28 3510 1
## beta5 -0.18 0.00 0.15 -0.48 -0.28 -0.18 -0.08 0.12 3091 1
## ln_sigma 0.54 0.00 0.15 0.23 0.44 0.54 0.64 0.80 2919 1
## sigma 0.59 0.00 0.09 0.45 0.53 0.58 0.65 0.80 2910 1
## lp__ 0.77 0.05 2.13 -4.39 -0.35 1.16 2.32 3.77 1635 1
##
## Samples were drawn using NUTS(diag_e) at Tue Nov 14 10:09:41 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
post <- extract(fit)
traceplot(fit,int_sample=TRUE,pairs=c('alpha','beta'))
## Warning: Ignoring unknown parameters: int_sample, pairs
pairs(fit)
ttl = list(
expression(beta["com"]),
expression(beta["pri"]),
expression(beta["lea"]),
expression(beta["rai"]),
expression(beta["cri"]),
expression(beta["adv"]),
expression(sigma)
)
par(mfrow=c(2,2))
for (i in 1:6) {
plot(1:length(post[[i]]),post[[i]],xlab='Iterations',
ylab='Parameter values',type='l',main=ttl[[i]])
}
for (i in 1:6) {
prm = post[[i]]
den = density(prm)
mod = den$x[which(den$y==max(den$y))]
hist(prm,breaks=20,border='white',col='grey',
freq=F,main=ttl[[i]],xlab='Parameter values')
mod = den$x[which(den$y==max(den$y))]
segments(mod,0,mod,max(den$y),col='blue',lwd=2)
qnt = quantile(prm,c(.025,.975))
segments(qnt[1],0,qnt[2],0,col='blue',lwd=2)
val = seq(-4,4,length=1000)
if (i<4) lines(val,dnorm(val,0,1),lty=2)
else lines(val,dunif(val,0,10),lty=2)
}
