Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Statistical Rethinking 2023 - Lecture 08

Statistical Rethinking 2023 - Lecture 08

Richard McElreath

January 25, 2023
Tweet

More Decks by Richard McElreath

Other Decks in Education

Transcript

  1. Statistical Rethinking
    8. Markov chain Monte Carlo
    2023

    View Slide

  2. View Slide

  3. 15 20 25 30
    6 8 10 12
    Marriage rate
    Divorce rate
    AL
    AK
    AR
    CO
    CT
    DE
    DC
    GA
    HI
    ID
    KY
    ME
    MN
    NJ
    ND
    OK
    RI
    TN
    UT
    VA
    WY
    Age at marriage
    M D
    A
    ?
    Marriage Divorce
    23 24 25 26 27 28 29
    6 8 10 12
    Median age of marriage
    Divorce rate
    AL
    AR
    CT
    DC
    ID
    ME
    MA
    MN
    NJ
    ND
    OK
    RI
    UT WY
    23 24 25 26 27 28 29
    15 20 25 30
    Median age of marriage
    Marriage rate
    AK
    AR
    DE
    DC
    HI
    ID
    ME
    MA
    MN
    NJ
    NY
    ND
    OK
    PA
    RI
    UT
    WY

    View Slide

  4. M D
    A
    M* D*
    A*
    Marriage
    observed
    Age at marriage
    observed
    Divorce
    observed

    View Slide

  5. M D
    A
    M* D*
    A*
    Marriage
    observed
    Age at marriage
    observed
    Divorce
    observed
    P
    Population

    View Slide

  6. K
    Knowledge
    Real, latent modeling problems

    View Slide

  7. K S
    Score
    Knowledge
    Real, latent modeling problems

    View Slide

  8. K S
    Score
    Q
    Knowledge
    Test/Instrument
    Real, latent modeling problems

    View Slide

  9. K S
    Score
    Q
    Knowledge
    Test/Instrument
    X
    Person
    feature
    Real, latent modeling problems

    View Slide

  10. Q Y
    Evaluation
    R
    Quality
    Referee
    X
    Feature
    Real, latent modeling problems

    View Slide

  11. Social networks
    YAB
    RAB
    Relationship A→B
    Behavior A→B

    View Slide

  12. Social networks
    YAB
    RAB
    RBA
    Relationship B→A
    Relationship A→B
    Behavior A→B

    View Slide

  13. Social networks
    YAB
    XA
    XB
    RAB
    RBA
    Relationship B→A
    Relationship A→B
    Features of B
    Features of A

    View Slide

  14. Social networks
    YAB
    XA
    XB
    RAB
    RBA
    Relationship B→A
    Relationship A→B
    Features of B
    Features of A

    View Slide

  15. Problems & solutions
    Real research problems:
    Many unknowns
    Nested relationships
    Bounded outcomes
    Di cult calculations
    K S
    Q
    X
    YAB
    XA
    XB
    RAB
    RBA

    View Slide

  16. aNAlYzE tHe DatA

    View Slide

  17. Computing the posterior
    1. Analytical approach (o en impossible)
    2. Grid approximation (very intensive)
    3. Quadratic approximation (limited)
    4. Markov chain Monte Carlo (intensive)

    View Slide

  18. King Markov

    View Slide

  19. e Metropolis Archipelago

    View Slide

  20. Contract: King Markov must visit each island in
    proportion to its population size
    Here’s how he does it...

    View Slide

  21. (1) Flip a coin to choose island on le or right.
    Call it the “proposal” island.
    1
    2
    1
    2

    View Slide

  22. 1 2 3 4 5 6 7
    p5
    proposal
    (1) Flip a coin to choose island on le or right.
    Call it the “proposal” island.
    (2) Find population of proposal island.

    View Slide

  23. p5
    (1) Flip a coin to choose island on le or right.
    Call it the “proposal” island.
    p4
    1 2 3 4 5 6 7
    proposal
    (3) Find population of current island.
    (2) Find population of proposal island.

    View Slide

  24. (2) Find population of proposal island.
    (1) Flip a coin to choose island on le or right.
    Call it the “proposal” island.
    p5
    p4
    1 2 3 4 5 6 7
    proposal
    (3) Find population of current island.
    (4) Move to proposal, with probability = p5
    p4

    View Slide

  25. (2) Find population of proposal island.
    (1) Flip a coin to choose island on le or right.
    Call it the “proposal” island.
    (3) Find population of current island.
    1 2 3 4 5 6 7
    (4) Move to proposal, with probability = p5
    p4
    (5) Repeat from (1)

    View Slide

  26. (2) Find population of proposal island.
    (1) Flip a coin to choose island on le or right.
    Call it the “proposal” island.
    (3) Find population of current island.
    (4) Move to proposal, with probability = p5
    p4
    1 2 3 4 5 6 7
    (5) Repeat from (1)
    is procedure ensures visiting each island in
    proportion to its population, in the long run.

    View Slide

  27. View Slide

  28. View Slide

  29. Markov chain Monte Carlo
    Usual use: Draw samples from a posterior
    distribution
    “Islands”: parameter values
    “Population size”: posterior probability
    Visit each parameter value in proportion to its
    posterior probability
    Any number of dimensions (parameters)

    View Slide

  30. “Markov chain Monte Carlo”
    Chain: Sequence of draws from distribution
    Markov chain: History doesn’t matter, just
    where you are now
    Monte Carlo: Random simulation

    View Slide

  31. “Markov chain Monte Carlo”
    Metropolis algorithm: Simple version of Markov
    chain Monte Carlo (MCMC)
    Easy to write, very general, o en ine cient

    View Slide

  32. Metropolis, Rosenbluth, Rosenbluth, Teller and Teller (1953)
    OURNAL OF CHEMICAL PHYSICS VOLUME 21, NUMBER 6 JUNE
    Equation of State Calculations by Fast Computing Machines
    NICHOLAS METROPOLIS, ARIANNA W. ROSENBLUTH, MARSHALL N. ROSENBLUTH, AND AUGUSTA H. TELLER,
    Los Alamos Scientific Laboratory, Los Alamos, New Mexico
    AND
    EDWARD TELLER, * Department of Physics, University of Chicago, Chicago, Illinois
    (Received March 6, 1953)
    A general method, suitable for fast computing machines, for investigatiflg such properties as equations of
    state for substances consisting of interacting individual molecules is described. The method consists of a
    modified Monte Carlo integration over configuration space. Results for the two-dimensional rigid-sphere
    system have been obtained on the Los Alamos MANIAC and are presented here. These results are compared
    to the free volume equation of state and to a four-term virial coefficient expansion.

    View Slide

  33. Metropolis, Rosenbluth, Rosenbluth, Teller and Teller (1953)
    OURNAL OF CHEMICAL PHYSICS VOLUME 21, NUMBER 6 JUNE
    Equation of State Calculations by Fast Computing Machines
    NICHOLAS METROPOLIS, ARIANNA W. ROSENBLUTH, MARSHALL N. ROSENBLUTH, AND AUGUSTA H. TELLER,
    Los Alamos Scientific Laboratory, Los Alamos, New Mexico
    AND
    EDWARD TELLER, * Department of Physics, University of Chicago, Chicago, Illinois
    (Received March 6, 1953)
    A general method, suitable for fast computing machines, for investigatiflg such properties as equations of
    state for substances consisting of interacting individual molecules is described. The method consists of a
    modified Monte Carlo integration over configuration space. Results for the two-dimensional rigid-sphere
    system have been obtained on the Los Alamos MANIAC and are presented here. These results are compared
    to the free volume equation of state and to a four-term virial coefficient expansion.

    View Slide

  34. Arianna Rosenbluth
    (1927–2020)

    View Slide

  35. Arianna Rosenbluth
    (1927–2020)

    View Slide

  36. View Slide

  37. Mathematical Analyzer, Numerical
    Integrator, and Computer
    MANIAC:
    1000 pounds
    5 kilobytes of memory
    70k multiplications/sec
    Your laptop:
    4–7 pounds
    8+ million kilobytes memory
    Billions of multiplications/sec

    View Slide

  38. MCMC is diverse
    Metropolis has yielded to newer, more
    e cient algorithms
    Many innovations in the last decades
    Best methods use gradients
    We’ll use Hamiltonian Monte Carlo

    View Slide

  39. Basic Rosenbluth (aka Metropolis) algorithm

    View Slide

  40. Basic Rosenbluth (aka Metropolis) algorithm
    large step size small step size

    View Slide

  41. -3 -2 -1 0 1 2 3
    0.0 0.1 0.2 0.3 0.4
    value
    probability

    View Slide

  42. -3 -2 -1 0 1 2 3
    0.0 0.1 0.2 0.3 0.4
    value
    probability
    -3 -2 -1 0 1 2 3
    1 2 3 4 5
    value
    -log probability

    View Slide

  43. x
    y
    z
    x
    y
    z

    View Slide

  44. View Slide

  45. Low probability
    High probability

    View Slide

  46. Hamiltonian Monte Carlo

    View Slide

  47. Hamiltonian Monte Carlo

    View Slide

  48. View Slide

  49. 0WFSUIJOLJOH )BNJMUPOJBO .POUF $BSMP JO UIF SBX ćF ).$ BMHPSJUIN OFFET ĕWF UIJOHT UP HP

    B GVODUJPO U UIBU SFUVSOT UIF OFHBUJWF MPHQSPCBCJMJUZ PG UIF EBUB BU UIF DVSSFOU QPTJUJPO QBSBNFUFS
    WBMVFT

    B GVODUJPO grad_U UIBU SFUVSOT UIF HSBEJFOU PG UIF OFHBUJWF MPHQSPCBCJMJUZ BU UIF DVSSFOU
    QPTJUJPO
    B TUFQ TJ[F epsilon
    B DPVOU PG MFBQGSPH TUFQT L BOE
    B TUBSUJOH QPTJUJPO current_q
    ,FFQ JO NJOE UIBU UIF QPTJUJPO JT B WFDUPS PG QBSBNFUFS WBMVFT BOE UIBU UIF HSBEJFOU BMTP OFFET UP
    SFUVSO B WFDUPS PG UIF TBNF MFOHUI 4P UIBU UIFTF U BOE grad_U GVODUJPOT NBLF NPSF TFOTF MFUT
    QSFTFOU UIFN ĕSTU CVJMU DVTUPN GPS UIF % (BVTTJBO FYBNQMF ćF U GVODUJPO KVTU FYQSFTTFT UIF MPH
    QPTUFSJPS BT TUBUFE CFGPSF JO UIF NBJO UFYU
    J
    MPH Q(ZJ|µZ, ) +
    J
    MPH Q(YJ|µY, ) + MPH Q(µZ|, .) + MPH Q(µY, , .)
    4P JUT KVTU GPVS DBMMT UP dnorm SFBMMZ
    3 DPEF
    # U needs to return neg-log-probability
    U <- function( q , a=0 , b=1 , k=0 , d=1 ) {
    muy <- q[1]
    mux <- q[2]
    U <- sum( dnorm(y,muy,1,log=TRUE) ) + sum( dnorm(x,mux,1,log=TRUE) ) +
    dnorm(muy,a,b,log=TRUE) + dnorm(mux,k,d,log=TRUE)
    return( -U )
    }
    /PX UIF HSBEJFOU GVODUJPO SFRVJSFT UXP QBSUJBM EFSJWBUJWFT -VDLJMZ (BVTTJBO EFSJWBUJWFT BSF WFSZ
    DMFBO ćF EFSJWBUJWF PG UIF MPHBSJUIN PG BOZ VOJWBSJBUF (BVTTJBO XJUI NFBO B BOE TUBOEBSE EFWJBUJPO
    C XJUI SFTQFDU UP B JT
    ∂ MPH /(Z|B, C)
    ∂B =
    Z − B
    C
    "OE TJODF UIF EFSJWBUJWF PG B TVN JT B TVN PG EFSJWBUJWFT UIJT JT BMM XF OFFE UP XSJUF UIF HSBEJFOUT
    ∂6
    ∂µY
    =
    ∂ MPH /(Y|µY, )
    ∂µY
    +
    ∂ MPH /(µY|, .)
    ∂µY
    =
    J
    YJ − µY

    +
    − µY
    .
    Pages 276–278

    View Slide

  50. Calculus is a superpower
    Hamiltonian Monte Carlo needs gradients
    How does it get them? Write them yourself or…
    Auto-di : Automatic di erentiation
    Symbolic derivatives of your model code
    Used in many machine learning approaches;
    “Backpropagation” is special case

    View Slide

  51. Stanisław Ulam (1909–1984)
    mc-stan.org

    View Slide

  52. Stanisław Ulam and his daughter Claire with MANIAC

    View Slide

  53. PAUSE

    View Slide

  54. Example: e Judgment at Princeton

    View Slide

  55. library(rethinking)
    data(Wines2012)
    Example: e Judgment at Princeton

    View Slide

  56. library(rethinking)
    data(Wines2012)
    Example: e Judgment at Princeton
    20 wines (10 French, 10 NJ)
    9 French & American judges
    180 scores

    View Slide

  57. Q S
    Score
    J
    Wine
    quality
    Judge

    View Slide

  58. Q S
    Score
    J
    Wine
    quality
    Judge
    X
    Wine origin

    View Slide

  59. Q S
    Score
    J
    Wine
    quality
    Judge
    X
    Wine origin
    Z
    Judge origin

    View Slide

  60. Q S
    Score
    J
    Wine
    quality
    Judge
    X
    Wine origin
    Z
    Judge origin

    View Slide

  61. Q S
    Score
    J
    Wine
    quality
    Judge
    X
    Wine origin
    Z
    Judge origin

    View Slide

  62. Q S
    Score
    J
    Wine
    quality
    Judge
    X
    Wine origin
    Z
    Judge origin
    Estimand: Association between wine quality and
    wine origin. Stratify by judge for e ciency.

    View Slide

  63. library(rethinking)
    data(Wines2012)
    d <- Wines2012
    dat <- list(
    S=standardize(d$score),
    J=as.numeric(d$judge),
    W=as.numeric(d$wine),
    X=ifelse(d$wine.amer==1,1,2),
    Z=ifelse(d$judge.amer==1,1,2)
    )
    mQ <- ulam(
    alist(
    S ~ dnorm(mu,sigma),
    mu <- Q[W],
    Q[W] ~ dnorm(0,1),
    sigma ~ dexp(1)
    ) , data=dat , chains=4 , cores=4 )
    precis(mQ,2)
    Q S
    J
    X Z
    Si
    ⇠ Normal(µi, )
    µi = QW [i]
    Qj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAACZnicbVHLS+NAHJ5kfdZXtSwevAwWRUFKooJ7WSiK4EksWis0IUym03Z0HmHml8US8k/ubc9e/DN2Gnvw9YOBL99jMvNNmgluIQj+ef6PufmFxaXl2srq2vpGfXPr3urcUNalWmjzkBLLBFesCxwEe8gMIzIVrJc+XUz13h9mLNfqDiYZiyUZKT7klICjknp5m3C8H1kucQTsGYprbSQR5UEk84QfOWEkyWEU1apvvP8bd5Ki1+dx6bhO8vhtNjgKq0gV/uC4fM60Ygr41BYeJvVm0AqqwV9BOANNNJubpP43GmiaS7cFFcTafhhkEBfEAKeClbUotywj9ImMWN9BRSSzcVHVVOI9xwzwUBu3FOCKfZ8oiLR2IlPnlATG9rM2Jb/T+jkMf8UFV1kOTNG3Hw1zgUHjaed4wA2jICYOEGq4OyumY2IIBfcyNVdC+PnKX8H9cSs8aR13Tpvt81kdS2gH7aIDFKIz1EZX6AZ1EUUv3rK35TW8V3/d/+lvv1l9b5ZpoA/j4//dXrZe

    View Slide

  64. library(rethinking)
    data(Wines2012)
    d <- Wines2012
    dat <- list(
    S=standardize(d$score),
    J=as.numeric(d$judge),
    W=as.numeric(d$wine),
    X=ifelse(d$wine.amer==1,1,2),
    Z=ifelse(d$judge.amer==1,1,2)
    )
    mQ <- ulam(
    alist(
    S ~ dnorm(mu,sigma),
    mu <- Q[W],
    Q[W] ~ dnorm(0,1),
    sigma ~ dexp(1)
    ) , data=dat , chains=4 , cores=4 )
    precis(mQ,2)
    Q S
    J
    X Z
    Si
    ⇠ Normal(µi, )
    µi = QW [i]
    Qj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAACZnicbVHLS+NAHJ5kfdZXtSwevAwWRUFKooJ7WSiK4EksWis0IUym03Z0HmHml8US8k/ubc9e/DN2Gnvw9YOBL99jMvNNmgluIQj+ef6PufmFxaXl2srq2vpGfXPr3urcUNalWmjzkBLLBFesCxwEe8gMIzIVrJc+XUz13h9mLNfqDiYZiyUZKT7klICjknp5m3C8H1kucQTsGYprbSQR5UEk84QfOWEkyWEU1apvvP8bd5Ki1+dx6bhO8vhtNjgKq0gV/uC4fM60Ygr41BYeJvVm0AqqwV9BOANNNJubpP43GmiaS7cFFcTafhhkEBfEAKeClbUotywj9ImMWN9BRSSzcVHVVOI9xwzwUBu3FOCKfZ8oiLR2IlPnlATG9rM2Jb/T+jkMf8UFV1kOTNG3Hw1zgUHjaed4wA2jICYOEGq4OyumY2IIBfcyNVdC+PnKX8H9cSs8aR13Tpvt81kdS2gH7aIDFKIz1EZX6AZ1EUUv3rK35TW8V3/d/+lvv1l9b5ZpoA/j4//dXrZe
    mean sd 5.5% 94.5% n_eff Rhat4
    Q[1] 0.14 0.30 -0.34 0.64 2962 1
    Q[2] 0.11 0.32 -0.40 0.61 3033 1
    Q[3] 0.27 0.31 -0.21 0.75 2608 1
    Q[4] 0.55 0.33 0.04 1.09 2598 1
    Q[5] -0.13 0.32 -0.63 0.37 2726 1
    Q[6] -0.37 0.32 -0.88 0.12 2710 1
    Q[7] 0.29 0.32 -0.22 0.81 2679 1
    Q[8] 0.27 0.32 -0.24 0.79 3248 1
    Q[9] 0.08 0.30 -0.40 0.56 3410 1
    Q[10] 0.12 0.32 -0.40 0.63 2697 1
    Q[11] -0.02 0.31 -0.52 0.50 2740 1
    Q[12] -0.03 0.32 -0.53 0.48 2809 1
    Q[13] -0.11 0.32 -0.62 0.39 3225 1
    Q[14] 0.01 0.33 -0.51 0.52 2690 1
    Q[15] -0.22 0.31 -0.71 0.28 2754 1
    Q[16] -0.21 0.32 -0.72 0.30 2120 1
    Q[17] -0.14 0.31 -0.64 0.37 3831 1
    Q[18] -0.86 0.32 -1.37 -0.35 3093 1
    Q[19] -0.16 0.30 -0.65 0.31 2783 1
    Q[20] 0.38 0.32 -0.12 0.88 2873 1
    sigma 1.00 0.06 0.91 1.09 2759 1

    View Slide

  65. Drawing the Markov Owl
    Complex machinery, but lots of diagnostics
    (1) Trace plots
    (2) Trace rank plots
    (3) R-hat convergence measure
    (4) Number of e ective samples
    (5) Divergent transitions

    View Slide

  66. 0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2962
    Q[1]
    0 200 400 600 800 1000
    -1.5 0.0 1.0
    n_eff = 3033
    Q[2]
    0 200 400 600 800 1000
    -2 -1 0 1
    n_eff = 2608
    Q[3]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2598
    Q[4]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2726
    Q[5]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2710
    Q[6]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2679
    Q[7]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3248
    Q[8]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3410
    Q[9]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2697
    Q[10]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2740
    Q[11]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2809
    Q[12]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 3225
    Q[13]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2690
    Q[14]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 2754
    Q[15]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2120
    Q[16]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3831
    Q[17]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 3093
    Q[18]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2783
    Q[19]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2873
    Q[20]
    Trace plots

    View Slide

  67. 0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2962
    Q[1]
    0 200 400 600 800 1000
    -1.5 0.0 1.0
    n_eff = 3033
    Q[2]
    0 200 4
    -2 -1 0 1
    Q[3]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2726
    Q[5]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2710
    Q[6]
    0 200 4
    -1.0 0.0 1.0
    Q[7]
    1.0
    n_eff = 3410
    Q[9]
    1.0
    n_eff = 2697
    Q[10]
    1.0
    Q[11]

    View Slide

  68. 800 1000 0 200 400 600 800 1000
    -1.5
    0 20
    -2 -
    800 1000
    eff = 2726
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2710
    Q[6]
    0 20
    -1.0 0.0 1.0
    Q[7]
    eff = 3410 n_eff = 2697
    Q[10]
    1.0
    Q[11]
    Warmup Samples

    View Slide

  69. Need more than 1 chain to
    check convergence
    Convergence: Each chain
    explores the right distribution
    and every chain explores the
    same distribution
    library(rethinking)
    data(Wines2012)
    d <- Wines2012
    dat <- list(
    S=standardize(d$score),
    J=as.numeric(d$judge),
    W=as.numeric(d$wine),
    X=ifelse(d$wine.amer==1,1,2),
    Z=ifelse(d$judge.amer==1,1,2)
    )
    mQ <- ulam(
    alist(
    S ~ dnorm(mu,sigma),
    mu <- Q[W],
    Q[W] ~ dnorm(0,1),
    sigma ~ dexp(1)
    ) , data=dat , chains=4 , cores=4 )
    precis(mQ,2)

    View Slide

  70. 0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2962
    Q[1]
    0 200 400 600 800 1000
    -1.5 0.0 1.0
    n_eff = 3033
    Q[2]
    0 200 400 600 800 1000
    -2 -1 0 1
    n_eff = 2608
    Q[3]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2598
    Q[4]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2726
    Q[5]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2710
    Q[6]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2679
    Q[7]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3248
    Q[8]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3410
    Q[9]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2697
    Q[10]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2740
    Q[11]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2809
    Q[12]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 3225
    Q[13]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2690
    Q[14]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 2754
    Q[15]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2120
    Q[16]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3831
    Q[17]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 3093
    Q[18]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2783
    Q[19]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2873
    Q[20]

    View Slide

  71. 0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2962
    Q[1]
    0 200 400 600 800 1000
    -1.5 0.0 1.0
    n_eff = 3033
    Q[2]
    0 200 400 600 800 1000
    -2 -1 0 1
    n_eff = 2608
    Q[3]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2598
    Q[4]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2726
    Q[5]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2710
    Q[6]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2679
    Q[7]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3248
    Q[8]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3410
    Q[9]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2697
    Q[10]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2740
    Q[11]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2809
    Q[12]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 3225
    Q[13]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2690
    Q[14]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 2754
    Q[15]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2120
    Q[16]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3831
    Q[17]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 3093
    Q[18]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2783
    Q[19]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2873
    Q[20]

    View Slide

  72. 0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2962
    Q[1]
    0 200 400 600 800 1000
    -1.5 0.0 1.0
    n_eff = 3033
    Q[2]
    0 200 400 600 800 1000
    -2 -1 0 1
    n_eff = 2608
    Q[3]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2598
    Q[4]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2726
    Q[5]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2710
    Q[6]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2679
    Q[7]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3248
    Q[8]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3410
    Q[9]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2697
    Q[10]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 2740
    Q[11]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2809
    Q[12]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 3225
    Q[13]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 2690
    Q[14]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5 1.5
    n_eff = 2754
    Q[15]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2120
    Q[16]
    0 200 400 600 800 1000
    -1.0 0.0 1.0
    n_eff = 3831
    Q[17]
    0 200 400 600 800 1000
    -2.0 -0.5 1.0
    n_eff = 3093
    Q[18]
    0 200 400 600 800 1000
    -1.5 -0.5 0.5
    n_eff = 2783
    Q[19]
    0 200 400 600 800 1000
    -1.0 0.0 1.0 2.0
    n_eff = 2873
    Q[20]

    View Slide

  73. 0 200 400 600 800 1000
    -10 0 5 10
    n_eff = 3946
    Q[1]
    0 200 400 600 800 1000
    -10 0 5 10
    n_eff = 3842
    Q[2]
    0 200 400 600 800 1000
    -200 0 200
    n_eff = 114
    Q[3]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 178
    Q[4]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 141
    Q[5]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 145
    Q[6]
    0 200 400 600 800 1000
    -200 0 200
    n_eff = 114
    Q[7]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 200
    Q[8]
    0 200 400 600 800 1000
    -200 0 200
    n_eff = 133
    Q[9]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 135
    Q[10]
    0 200 400 600 800 1000
    -200 200
    n_eff = 146
    Q[11]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 230
    Q[12]
    0 200 400 600 800 1000
    -400 -100 200
    n_eff = 112
    Q[13]
    0 200 400 600 800 1000
    -400 0 200
    n_eff = 179
    Q[14]
    0 200 400 600 800 1000
    -400 0 200
    n_eff = 205
    Q[15]
    0 200 400 600 800 1000
    -400 -100 200
    n_eff = 186
    Q[16]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 108
    Q[17]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 119
    Q[18]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 196
    Q[19]
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 92
    Q[20]

    View Slide

  74. 600 800 1000 0 200 400 600 800 1000
    -300 0
    0 200 400 6
    -200 0 2
    600 800 1000
    n_eff = 133
    0 200 400 600 800 1000
    -300 0 200
    n_eff = 135
    Q[10]
    0 200 400 6
    -200 200
    Q[11]
    n_eff = 112
    0 200
    n_eff = 179
    Q[14]
    0 200
    Q[15]

    View Slide

  75. Trace rank (Trank) plots
    n_eff = 2962
    Q[1] n_eff = 3033
    Q[2] Q
    n_eff = 2726
    Q[5] n_eff = 2710
    Q[6] Q
    Rank orders of samples. No
    chain should tend to be below/
    above others.

    View Slide

  76. n_eff = 2962
    Q[1] n_eff = 3033
    Q[2] n_eff = 2608
    Q[3] n_eff = 2598
    Q[4]
    n_eff = 2726
    Q[5] n_eff = 2710
    Q[6] n_eff = 2679
    Q[7] n_eff = 3248
    Q[8]
    n_eff = 3410
    Q[9] n_eff = 2697
    Q[10] n_eff = 2740
    Q[11] n_eff = 2809
    Q[12]
    n_eff = 3225
    Q[13] n_eff = 2690
    Q[14] n_eff = 2754
    Q[15] n_eff = 2120
    Q[16]
    n_eff = 3831
    Q[17] n_eff = 3093
    Q[18] n_eff = 2783
    Q[19] n_eff = 2873
    Q[20]

    View Slide

  77. n_eff = 3946
    Q[1] n_eff = 3842
    Q[2] n_eff = 114
    Q[3] n_eff = 178
    Q[4]
    n_eff = 141
    Q[5] n_eff = 145
    Q[6] n_eff = 114
    Q[7] n_eff = 200
    Q[8]
    n_eff = 133
    Q[9] n_eff = 135
    Q[10] n_eff = 146
    Q[11] n_eff = 230
    Q[12]
    n_eff = 112
    Q[13] n_eff = 179
    Q[14] n_eff = 205
    Q[15] n_eff = 186
    Q[16]
    n_eff = 108
    Q[17] n_eff = 119
    Q[18] n_eff = 196
    Q[19] n_eff = 92
    Q[20]

    View Slide

  78. n_eff = 145
    Q[6] n_eff = 114
    Q[7]
    n_eff = 135
    Q[10] n_eff = 146
    Q[11]

    View Slide

  79. R-hat
    When chains converge:
    (1) Start and end of each chain explores
    same region
    (2) Independent chains explore same region
    R-hat is a ratio of variances:
    As total variance shrinks to average
    variance within chains, R-hat approaches 1
    NO GUARANTEES; NOT A TEST
    0 200 400 600 800 1000
    0.00 0.10 0.20
    sample
    variance
    between
    within
    mean sd 5.5% 94.5% n_eff Rhat4
    Q[1] 0.14 0.30 -0.34 0.64 2962 1
    Q[2] 0.11 0.32 -0.40 0.61 3033 1
    Q[3] 0.27 0.31 -0.21 0.75 2608 1

    View Slide

  80. n_e
    Estimate of number of e ective
    samples
    “How long would the chain be, if
    each sample was independent of
    the one before it?”
    When samples are autocorrelated,
    you have fewer e ective samples
    0 5 10 15 20 25 30
    0.0 0.4 0.8
    Lag
    ACF
    Series post[, 1, 1]
    mean sd 5.5% 94.5% n_eff Rhat4
    Q[1] 0.14 0.30 -0.34 0.64 2962 1
    Q[2] 0.11 0.32 -0.40 0.61 3033 1
    Q[3] 0.27 0.31 -0.21 0.75 2608 1

    View Slide

  81. mQO <- ulam(
    alist(
    S ~ dnorm(mu,sigma),
    mu <- Q[W] + O[X],
    Q[W] ~ dnorm(0,1),
    O[X] ~ dnorm(0,1),
    sigma ~ dexp(1)
    ) , data=dat , chains=4 , cores=4 )
    plot(precis(mQO,2))
    Q S
    J
    X Z
    Si
    ⇠ Normal(µi, )
    µi = QW [i]
    + OX[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAACkHicfVHbSiNBEO0ZL+vGW1Yf96UxKIoSZlTQhxVdRZB9WA0aE8gMQ0+nE3vty9BdI4Yh37P/45t/Y0+SB29Y0HD6nFNV3VVpJriFIHj2/Knpmdlvc98r8wuLS8vVHyu3VueGsibVQpt2SiwTXLEmcBCsnRlGZCpYK70/K/XWAzOWa3UDg4zFkvQV73FKwFFJ9f91wvFGZLnEEbBHKP5qI4kYbkYyT/iOE/qSbEVRZXTHG0e4kRStDo+HeBtfJkW7hE5uJP8+LRPshGX25dfyuM0bx/ljphVTwEtbuJVUa0E9GAX+CMIJqKFJXCXVp6iraS5dCSqItZ0wyCAuiAFOBRtWotyyjNB70mcdBxWRzMbFaKBDvO6YLu5p444CPGJfZxREWjuQqXNKAnf2vVaSn2mdHHqHccFVlgNTdNyolwsMGpfbwV1uGAUxcIBQw91bMb0jhlBwO6y4IYTvv/wR3O7Ww736bmO/dnI6Gccc+onW0CYK0QE6QRfoCjUR9Ra9Pe+Xd+Sv+If+sf97bPW9Sc4qehP+nxfJIcPT

    View Slide

  82. Q S
    J
    X Z
    Si
    ⇠ Normal(µi, )
    µi = QW [i]
    + OX[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAACkHicfVHbSiNBEO0ZL+vGW1Yf96UxKIoSZlTQhxVdRZB9WA0aE8gMQ0+nE3vty9BdI4Yh37P/45t/Y0+SB29Y0HD6nFNV3VVpJriFIHj2/Knpmdlvc98r8wuLS8vVHyu3VueGsibVQpt2SiwTXLEmcBCsnRlGZCpYK70/K/XWAzOWa3UDg4zFkvQV73FKwFFJ9f91wvFGZLnEEbBHKP5qI4kYbkYyT/iOE/qSbEVRZXTHG0e4kRStDo+HeBtfJkW7hE5uJP8+LRPshGX25dfyuM0bx/ljphVTwEtbuJVUa0E9GAX+CMIJqKFJXCXVp6iraS5dCSqItZ0wyCAuiAFOBRtWotyyjNB70mcdBxWRzMbFaKBDvO6YLu5p444CPGJfZxREWjuQqXNKAnf2vVaSn2mdHHqHccFVlgNTdNyolwsMGpfbwV1uGAUxcIBQw91bMb0jhlBwO6y4IYTvv/wR3O7Ww736bmO/dnI6Gccc+onW0CYK0QE6QRfoCjUR9Ra9Pe+Xd+Sv+If+sf97bPW9Sc4qehP+nxfJIcPT
    sigma
    O[2]
    O[1]
    Q[20]
    Q[19]
    Q[18]
    Q[17]
    Q[16]
    Q[15]
    Q[14]
    Q[13]
    Q[12]
    Q[11]
    Q[10]
    Q[9]
    Q[8]
    Q[7]
    Q[6]
    Q[5]
    Q[4]
    Q[3]
    Q[2]
    Q[1]
    -1.5 -1.0 -0.5 0.0 0.5 1.0
    Value
    wines
    origins

    View Slide

  83. mQOJ <- ulam(
    alist(
    S ~ dnorm(mu,sigma),
    mu <- (Q[W] + O[X] - H[J])*D[J],
    Q[W] ~ dnorm(0,1),
    O[X] ~ dnorm(0,1),
    H[J] ~ dnorm(0,1),
    D[J] ~ dexp(1),
    sigma ~ dexp(1)
    ) , data=dat , chains=4 )
    plot(precis(mQOJ,2))
    Q S
    J
    X Z
    Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==

    View Slide

  84. Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==
    expected score

    View Slide

  85. Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==
    expected score
    wine quality

    View Slide

  86. Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==
    expected score
    wine quality
    origin

    View Slide

  87. Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==
    expected score
    wine quality
    origin judge
    harshness
    judge
    discrimination

    View Slide

  88. Q S
    J
    X Z
    Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==
    sigma
    D[9]
    D[8]
    D[7]
    D[6]
    D[5]
    D[4]
    D[3]
    D[2]
    D[1]
    H[9]
    H[8]
    H[7]
    H[6]
    H[5]
    H[4]
    H[3]
    H[2]
    H[1]
    O[2]
    O[1]
    Q[20]
    Q[19]
    Q[18]
    Q[17]
    Q[16]
    Q[15]
    Q[14]
    Q[13]
    Q[12]
    Q[11]
    Q[10]
    Q[9]
    Q[8]
    Q[7]
    Q[6]
    Q[5]
    Q[4]
    Q[3]
    Q[2]
    Q[1]
    -3 -2 -1 0 1 2
    harshness
    discrimination

    View Slide

  89. Q S
    J
    X Z
    Si
    ⇠ Normal(µi
    , )
    µi = (QW [i]
    + OX[i]
    HJ[i]
    )DJ[i]
    Qj
    ⇠ Normal(0, 1)
    Oj
    ⇠ Normal(0, 1)
    Hj
    ⇠ Normal(0, 1)
    Dj
    ⇠ Normal(0, 1)
    ⇠ Exponential(1)
    AAAC5XicfVJNbxMxEPVuC5QUaIAjF4uIKhEl2m0rwaVSBa0UcaCNSppI2ZXldZzU1B8r24sarXLnwgFU9dr/1Fv/Sw94N3ugH8pIlp7fezMej52knBkbBNeev7T86PGTlae11WfPX6zVX746NirThPaI4koPEmwoZ5L2LLOcDlJNsUg47SenXwq9/5Nqw5T8bqcpjQWeSDZmBFtHofrNEWJwPTJMwMjSM5t/U1pgPmtGIkNswwkTgVtRVCv3cH0HNmEX5f0hi2fwPTxA+aCEH2AH5V9L2NqrkMvqoh8PVg82wqLowWK5s1jeWyzPe7/l2D9LlaTSssIWtlC9EbSDMuB9EFagAao4RPWraKRIJlwJwrExwzBIbZxjbRnhdFaLMkNTTE7xhA4dlFhQE+flK83gO8eM4Fhpt6SFJft/Ro6FMVOROKfA9sTc1QryIW2Y2fGnOGcyzSyVZH7QOOPQKlg8ORwxTYnlUwcw0cz1CskJ1phY9zFqbgjh3SvfB8eb7XCrvdndbux+rsaxAt6At6AJQvAR7IIOOAQ9QLzE++X98f76E/+3f+5fzK2+V+W8BrfCv/wHQKDiqw==
    -2 -1 0 1
    0.0 0.2 0.4 0.6 0.8
    origin contrast
    Density

    View Slide

  90. e Folk eorem of
    Statistical Computing
    “When you have computational
    problems, o en there’s a problem
    with your model.”
    Andrew Gelman
    Spider-Man of Bayesian data analysis

    View Slide

  91. View Slide

  92. Divergent transitions
    Divergent transition: A kind of
    rejected proposal
    Simulation diverges from true path
    Many DTs: poor exploration &
    possible bias
    Will discuss again in later lecture

    View Slide

  93. El Pueblo Unido
    Desktop MCMC has been a revolution in
    scienti c computing
    Custom scienti c modeling
    High-dimension
    Propagate measurement error
    Sir David Spiegelhalter

    View Slide

  94. Course Schedule
    Week 1 Bayesian inference Chapters 1, 2, 3
    Week 2 Linear models & Causal Inference Chapter 4
    Week 3 Causes, Confounds & Colliders Chapters 5 & 6
    Week 4 Over tting / MCMC Chapters 7, 8, 9
    Week 5 Generalized Linear Models Chapters 10, 11
    Week 6 Integers & Other Monsters Chapters 11 & 12
    Week 7 Multilevel models I Chapter 13
    Week 8 Multilevel models II Chapter 14
    Week 9 Measurement & Missingness Chapter 15
    Week 10 Generalized Linear Madness Chapter 16
    https://github.com/rmcelreath/stat_rethinking_2023

    View Slide

  95. View Slide