mu0 <- 1.9; t20 <- 0.95^2;  s20 <- 0.01; nu0 <- 1
y <- c(1.64, 1.70, 1.72, 1.74, 1.82, 1.82, 1.82, 1.90, 2.08)

G <- 100; H <- 100;

mean.grid <- seq(1.505, 2.00, length = G)
prec.grid <- seq(1.75, 175, length = H)
post.grid <- matrix(nrow = G, ncol = H)

for (g in 1:G) {
  for (h in 1:H) {
    post.grid[g,h] <- 
      dnorm(mean.grid[g], mu0, sqrt(t20))*
      dgamma(prec.grid[h], nu0/2,  s20*nu0/2)*
      prod(dnorm(y, mean.grid[g], 1/sqrt(prec.grid[h])))
  }
}
post.grid <- post.grid/sum(post.grid)
plot(mean.grid, rowSums(post.grid), type = "l", col = "blue", xlab = expression(theta), ylab = expression(p(theta)), main = "Distribusi posterior marjinal dari parameter theta")

plot(prec.grid, colSums(post.grid),type = "l", col = "red", xlab = expression(sigma^2), ylab = expression(p(sigma^2)), main = "Distribusi posterior marjinal dari parameter sigma^2")

library(plotly)

fig <- plot_ly(
  y = prec.grid, 
  x = mean.grid, 
  z = t(post.grid), 
  type = "contour"
)
fig %>% layout(title = "Distribusi peluang gabungan theta dan sigma^2",xaxis = list(title = "theta"), yaxis = list(title = "sigma^2"), legend = NULL)