Add conjugate priors plot

ch08.r

@@ -133,3 +133,40 @@ lines(t, prior, lwd = 4, col="green")
 lines(t, q, lwd = 4, col="purple")
 grid()
 
+## Conjugate priors
+
+# R code
+sigma0 = 0.25
+mu0 = 0.0
+
+mu = 1
+sigma = 0.25
+
+Nset = c(0, 1, 2, 5, 8, 12, 20)
+x0 = sigma * rnorm(100)
+posterior = list()
+
+for (i in 1:7) {
+    N = Nset[i]
+    x = x0[1:N]
+    t = seq(-1, 1.5, 2.5/1000)
+
+    p0 = dnorm(t, 0, 1)
+    theta = mu*(N*sigma0^2)/(N*sigma0^2+sigma^2) + mu0*(sigma^2)/(N*sigma0^2+sigma^2)
+    sigmaN = sqrt(1/(1/sigma0^2+N/sigma^2));
+    posterior[[i]] = dnorm(t, theta, sigmaN)
+}
+
+plot(t, posterior[[7]], type="n", xlab="", ylab="")
+lines(t, posterior[[1]], lwd=3, col="red")
+lines(t, posterior[[2]], lwd=3, col="orange")
+lines(t, posterior[[3]], lwd=3, col="yellow")
+lines(t, posterior[[4]], lwd=3, col="green")
+lines(t, posterior[[5]], lwd=3, col="turquoise")
+lines(t, posterior[[6]], lwd=3, col="blue")
+lines(t, posterior[[7]], lwd=3, col="purple")
+
+legend(-0.9, 7, legend=c("N = 0", "N = 1", "N = 2", "N = 5", "N = 8", "N = 12", "N = 20"), col=c("red", "orange", "yellow", "green", "turquoise", "blue", "purple"), lty=1:1, lwd=3)
+grid()
+
+#############