Skip to content

Instantly share code, notes, and snippets.

@rmcelreath

rmcelreath/prior_likelihood_conflict.r

Created September 11, 2023 11:29
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save rmcelreath/39dd410fc6bb758e54d79249b11eeb2f to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save rmcelreath/39dd410fc6bb758e54d79249b11eeb2f to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. rmcelreath revised this gist Sep 11, 2023. 1 changed file with 2 additions and 0 deletions.
    2 changes: 2 additions & 0 deletions prior_likelihood_conflict.r
    View file
    Original file line number Diff line number Diff line change
    @@ -1,5 +1,7 @@
    # prior - likelihood conflict

    library(rethinking)

    yobs <- 0

    mtt <- ulam(
  2. rmcelreath created this gist Sep 11, 2023.
    67 changes: 67 additions & 0 deletions prior_likelihood_conflict.r
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,67 @@
    # prior - likelihood conflict

    yobs <- 0

    mtt <- ulam(
    alist(
    y ~ dstudent(2,mu,1),
    mu ~ dstudent(2,10,1)
    ), data=list(y=yobs) , chains=4 , iter=2000 )
    mnn <- ulam(
    alist(
    y ~ dnorm(mu,1),
    mu ~ dnorm(10,1)
    ), data=list(y=yobs) , chains=4 , iter=2000)
    mtn <- ulam(
    alist(
    y ~ dstudent(2,mu,1),
    mu ~ dnorm(10,1)
    ), data=list(y=yobs) , chains=4 , iter=2000)
    mnt <- ulam(
    alist(
    y ~ dnorm(mu,1),
    mu ~ dstudent(2,10,1)
    ), data=list(y=yobs) , chains=4 , iter=2000)

    # plot
    par(mfrow=c(2,2),cex=1.05)
    ymax <- 0.53
    xlwd <- 1.5
    postcol <- 2
    xadj <- 0.8

    p <- extract.samples(mnn)
    dens(p$mu, xlim=c(-5,15), ylim=c(0,ymax), lwd=xlwd+1 , col=postcol, xlab="" ,adj=xadj )
    #mtext("normal prior, normal likelihood")
    curve( dnorm(yobs,x,1) , add=TRUE , lty=1 , lwd=xlwd ) # lik
    curve( dnorm(x,10,1) , add=TRUE , lty=2 , lwd=xlwd ) # prior

    text(0,0.42,"likelihood")
    text(10,0.42,"prior")

    p <- extract.samples(mtt)
    dens(p$mu , xlim=c(-5,15) , ylim=c(0,ymax) , lwd=xlwd+1 , col=postcol , xlab="" ,adj=xadj )
    #mtext("t prior, t likelihood")
    curve( dstudent(yobs,2,x,1) , add=TRUE , lty=1 , lwd=xlwd ) # lik
    curve( dstudent(x,2,10,1) , add=TRUE , lty=2 , lwd=xlwd ) # prior

    text(0,0.42,"likelihood")
    text(10,0.42,"prior")

    p <- extract.samples(mnt)
    dens(p$mu, xlim=c(-5,15), ylim=c(0,ymax), lwd=xlwd+1 , col=postcol, xlab="" ,adj=xadj )
    #mtext("t prior, normal likelihood")
    curve( dnorm(yobs,x,1) , add=TRUE , lty=1 , lwd=xlwd) # lik
    curve( dstudent(x,2,10,1) , add=TRUE , lty=2 , lwd=xlwd) # prior

    text(0,0.42,"likelihood")
    text(10,0.42,"prior")

    p <- extract.samples(mtn)
    dens(p$mu, xlim=c(-5,15), ylim=c(0,ymax), lwd=xlwd+1 , col=postcol, xlab="" ,adj=xadj )
    #mtext("normal prior, t likelihood")
    curve( dstudent(yobs,2,x,1) , add=TRUE , lty=1 , lwd=xlwd) # lik
    curve( dnorm(x,10,1) , add=TRUE , lty=2 , lwd=xlwd) # prior

    text(0,0.42,"likelihood")
    text(10,0.42,"prior")