Thuis
Contacten

    Hoofdpagina


Docent: drs. Rob Flohr

Dovnload 2 Mb.

Docent: drs. Rob Flohr



Pagina23/25
Datum05.12.2018
Grootte2 Mb.

Dovnload 2 Mb.
1   ...   17   18   19   20   21   22   23   24   25

node mean sd MC error 2.5% median 97.5% start sample

lambda[1] 219.6 14.78 0.2398 192.3 219.3 249.8 1 3000

lambda[2] 199.8 14.16 0.269 173.1 199.4 228.7 1 3000

lambda[3] 300.1 17.64 0.3531 267.5 299.4 336.0 1 3000

lambda[4] 179.9 13.09 0.242 155.3 179.4 208.1 1 3000

lambda[5] 240.1 15.6 0.3109 210.3 240.2 271.5 1 3000

predict[1] 220.4 21.28 0.332 180.0 220.0 262.0 1 3000

predict[2] 199.6 19.75 0.3913 162.0 199.0 239.0 1 3000

predict[3] 299.8 24.57 0.4608 253.0 300.0 349.0 1 3000

predict[4] 179.7 18.58 0.3586 144.0 179.0 218.0 1 3000

predict[5] 240.4 21.58 0.4263 200.0 240.0 284.0 1 3000

--------------------------------------------------------------------------------------------------------------------------------------


Voorbeeld: meerdere 'chains'

- bespreken odc-bestand Lesaffre hs. 8 i.v.m. 8 sets initial values (startwaarden voor de priors) en het volgen via 'Trace'.


-let op: voordat je op compile klikt eerst aantal chains invullen!
-let op: voordat je op Update klikt eerst op 'Trace' drukken!
# Osteoporosis regression example
model

{

for (i in 1:N){



tbbmc[i] ~dnorm(mu[i],tau)

mu[i] <- beta0+beta1*bmi[i]

}

sigma2 <- 1/tau



sigma <- sqrt(sigma2)
beta0 ~ dnorm(0,1.0E-6)

beta1 ~ dnorm(0,1.0E-6)

tau ~ dgamma(1.0E-3,1.0E-3)

}
# Data


list(tbbmc=c(1.798, 2.588, 2.325, 2.236, 1.925, 2.304, 2.183, 2.010, 1.579, 1.911, 1.681, 1.581, 2.138, 2.535, 1.480, 1.873, 1.680, 1.697, 1.861, 2.122, 1.870, 1.767, 2.501, 1.953, 1.340, 1.309, 1.711, 1.702, 2.176, 1.587, 1.756, 2.367, 1.714, 1.824, 2.074, 2.177, 2.004, 2.209, 1.432, 2.156, 2.410, 1.925, 2.050, 1.540, 1.478, 2.062, 1.531, 2.001, 1.562, 2.527, 1.606, 1.676, 1.849, 2.309, 2.504, 2.000, 2.052, 1.863, 1.771, 1.757, 1.721, 1.513, 1.770, 2.645, 1.709, 2.146, 1.904, 1.505, 1.407, 1.774, 1.709, 2.140, 1.894, 1.857, 1.978, 1.987, 1.942, 1.735, 1.851, 1.570, 1.455, 2.729, 2.032, 2.436, 1.614, 1.303, 2.125, 1.762, 1.513, 1.967, 1.470, 2.837, 1.336, 1.663, 1.778, 2.215, 1.675, 1.740, 2.021, 2.018, 2.051, 1.634, 1.559, 1.661, 1.683, 1.296, 1.571, 2.274, 1.911, 1.544, 1.684, 1.474, 2.016, 1.915, 1.581, 2.267, 2.210, 2.143, 1.424, 1.794, 1.499, 2.022, 2.205, 2.039, 1.748, 2.247, 2.243, 1.742, 1.809, 1.798, 1.975, 1.813, 1.435, 1.548, 2.107, 1.699, 1.781, 1.681, 1.901, 2.215, 2.003, 1.574, 2.699, 1.566, 2.133, 2.091, 1.490, 1.590, 1.401, 2.672, 2.415, 2.256, 1.775, 1.566, 2.159, 2.566, 1.874, 2.152, 1.820, 2.027, 1.498, 2.381, 2.002, 1.624, 1.470, 1.917, 2.595, 1.390, 1.779, 2.265, 1.885, 1.080, 1.837, 2.481, 1.898, 1.977, 2.022, 1.901, 0.579, 1.948, 1.863, 2.397, 1.717, 2.029, 1.924, 2.319, 1.638, 1.235, 1.769, 1.830, 1.745, 1.829, 1.934, 2.331, 1.756, 2.259, 1.439, 1.872, 1.796, 1.958, 1.829, 2.378, 1.826, 1.612, 2.485, 1.879, 2.032, 1.846, 1.467, 1.739, 1.597, 1.442, 1.999, 1.838, 1.773, 2.158, 1.458, 2.002, 1.738, 1.653, 1.462, 1.638, 2.201, 1.995, 1.704, 1.780, 2.463, 1.732, 1.270, 1.728, 2.183, 1.703, 1.505, 1.850),

bmi=c(23.61, 30.48, 27.18, 34.68, 26.72, 25.78, 29.24, 30.76, 21.64, 23.53, 22.58, 29.38,37.18, 32.60, 19.91, 28.00, 21.64, 25.08, 23.92, 27.34, 23.83, 32.44, 32.49, 29.05,20.17, 18.73, 22.31, 29.34, 29.52, 24.44, 28.40, 27.89, 25.11, 36.44, 26.40, 23.53,34.25, 29.59, 20.08, 28.52, 28.88, 23.24, 24.22, 23.56, 19.56, 30.78, 23.24, 22.77,22.37, 31.22, 24.39, 24.22, 24.84, 27.02, 35.96, 25.88, 23.07, 31.53, 22.77, 27.34,33.33, 25.57, 23.24, 39.84, 25.71, 34.96, 24.64, 21.79, 22.03, 30.46, 33.78, 25.31,25.00, 23.87, 28.16, 35.06, 27.13, 22.18, 23.80, 25.00, 17.97, 24.77, 26.93, 36.33,24.65, 28.54, 29.52, 28.48, 26.38, 23.81, 19.56, 30.47, 24.32, 24.03, 33.33, 32.89,21.23, 20.20, 23.53, 25.85, 26.22, 24.12, 26.22, 29.33, 23.53, 32.27, 25.21, 35.49,22.48, 26.38, 25.00, 28.62, 29.97, 30.78, 27.03, 27.78, 29.40, 39.26, 24.65, 24.39,34.10, 23.50, 26.67, 25.33, 19.38, 25.97, 24.30, 24.35, 24.84, 25.89, 29.00, 20.13,

20.83, 21.76, 26.44, 30.44, 24.46, 24.39, 31.14, 28.67, 24.84, 27.27, 28.00, 19.96,24.46, 28.55, 24.86, 21.93, 22.37, 31.22, 26.35, 24.91, 25.96, 23.62, 27.47, 25.43,27.34, 27.94, 27.63, 27.10, 20.55, 28.26, 25.81, 25.64, 23.07, 22.48, 30.08, 22.37,23.61, 25.10, 27.69, 17.83, 25.64, 29.75, 25.64, 27.24, 22.86, 25.89, 23.23, 22.64,31.64, 37.32, 22.55, 22.10, 27.34, 32.05, 22.66, 20.88, 24.31, 27.25, 22.67, 35.52,28.84, 29.30, 25.22, 30.84, 17.71, 27.34, 32.39, 22.89, 24.69, 28.96, 18.90, 22.66,33.12, 26.62, 30.06, 29.05, 20.93, 19.47, 20.39, 19.17, 29.17, 29.62, 21.45, 27.24,20.70, 29.27, 27.34, 21.78, 21.10, 27.47, 26.78, 23.12, 24.35, 28.62, 27.82, 21.76,31.35, 24.65, 37.46, 21.79, 18.99, 28.30), N=234)
# Initial values
list(beta0=0.4,beta1=0.025,tau=1/0.05)

list(beta0=0.4,beta1=0.055tau=1/0.05)

list(beta0=0.4,beta1=0.025,tau=1/0.12)

list(beta0=0.4,beta1=0.055,tau=1/0.12)


list(beta0=1.2,beta1=0.025,tau=1/0.05)

list(beta0=1.2,beta1=0.055,tau=1/0.05)

list(beta0=1.2,beta1=0.025,tau=1/0.12)

list(beta0=1.2,beta1=0.055,tau=1/0.12)

model is syntactically correct

data loaded

model compiled

chain initialized but other chain(s) contain uninitialized variables

chain initialized but other chain(s) contain uninitialized variables

chain initialized but other chain(s) contain uninitialized variables

chain initialized but other chain(s) contain uninitialized variables

chain initialized but other chain(s) contain uninitialized variables

chain initialized but other chain(s) contain uninitialized variables

chain initialized but other chain(s) contain uninitialized variables

model is initialized

Time series









Kernel density




Node statistics



node mean sd MC error 2.5% median 97.5% start sample

beta0 0.8225 0.1166 5.884E-4 0.5929 0.8225 1.05 1 40000

beta1 0.04004 0.004368 2.237E-5 0.03153 0.04004 0.04865 1 40000

tau 12.05 1.117 0.005411 9.95 12.01 14.36 1 40000


Huiswerkopgaven voor les 8 op 17 juni 2014

(1)
Bestudeer nogmaals hoofdstuk 5 uit De Bayesiaanse benadering.


(2)
Maak in Winbugs een model voor het volgende vraagstuk:

De Britse natuurkundige Henry Cavendish wilde in 1798 de dichtheid van de aarde bepalen en voerde daartoe 23 experimenten uit die de volgende gegevens m.b.t. de dichtheid van verschillende gesteenten opleverden:
5.36 5.29 5.58 5.65 5.57 5.53 5.62 5.29 5.44 5.34 5.79 5.10 5.27 5.39 5.42 5.47 5.63 5.34 5.46 5.30 5.78 5.68 5.85

a) maak een histogram van de data (als je dit in R doet, gebruik dan:


y=c(5.36, 5.29, ..............., 5.68,5.85)
hist(y, breaks=c(5, 5.16, 5.32, 5.48, 5.80, 5.96))

Foutje: er had moeten staan:
hist(y, breaks=c(5, 5.16, 5.32, 5.48, 5,64, 5.80, 5.96))

b) baseer je model op een normale verdeling (d.w.z., ga ervan uit dat de data gegenereerd worden door een normale verdeling: likelihood)

c) voor de eenvoud gaan we uit van een constante populatievariantie: 0.04. Omdat WinBUGS werkt met de zogeheten 'precision' (=reciproke van de variantie) moeten we eerst de precision van de likelihood bepalen: 1/0.04=25

d) men wist in die tijd al (dus voordat de experimenten uitgevoerd werden) dat de dichtheid varieerde tussen 2.5 en 7.5. Op basis daarvan mag je uitgaan van een prior die normaal verdeeld is met gemiddelde 5 en 'precision' 2 (variantie: 0.5)

e) bepaal op basis van Bayesiaanse statistiek en met behulp van WinBUGS een 95% credible interval voor de dichtheid van de aarde.

f) formuleer in woorden wat dit interval betekent

-----------------------------------------------------------------------------------------------------------------------------------



a)


e)

model {


for (i in 1:n) { y[i]~dnorm(theta,25) }

theta~dnorm(5,2)

}

list(y=c(5.36,5.29,5.58,5.65,5.57,5.53,5.62,5.29,5.44,5.34,5.79,5.10,5.27,5.39,5.42,5.47,5.63,5.34,5.46,5.30,5.78,5.68,5.85), n=23)



#experimenten om de dichtheid van de aarde te bepalen (Cavendish 1798)

model is syntactically correct

data loaded

model compiled



initial values generated, model initialized

Kernel density




Node statistics
1   ...   17   18   19   20   21   22   23   24   25

  • node mean sd MC error 2.5% median 97.5% start sample

  • Dovnload 2 Mb.