## Data loading
setwd("C:/Users/Hun/Desktop/")
data <- read.table("student_score.txt", header=T, sep=" ")


## Sample statistics
Sxy <- 0
Sx <- 0
Sy <- 0

for (i in 1:22) {
  Sxy <- Sxy+((data[i, 1]-mean(data[, 1]))*(data[i, 2]-mean(data[, 2])))
  Sx <- Sx+((data[i, 1]-mean(data[, 1]))^2)
  Sy <- Sy+((data[i, 2]-mean(data[, 2]))^2)
}


theta_hat <- Sxy / (sqrt(Sx) * sqrt(Sy))
# theta_hat <- cor(data$mech, data$vecs)

n <- nrow(data)


## Likelihood
h <- function(w, theta) {
  1/(cosh(w)-theta*theta_hat)^(n-1)
}

likelihood <- function(theta) {
  integ <- integrate(h, lower=0, upper=Inf, theta=theta)
  return(((n-2)*(1-theta^2)^{(n-1)/2}*(1-theta_hat^2)^{(n-4)/2}/pi) * integ$value)
}


## Prior
flat <- 1/2

jeff <- function(theta) {
  return(1/(1-theta^2))
}

tri <- function(theta) {
  return(1-abs(theta))
}


## Marginal
f_flat <- function(theta) {return(likelihood(theta) * 1/2)}
c_flat <- c()
for (k in seq(-0.999, 0.999, 0.001)) {c_flat <- c(c_flat, f_flat(k))}
m_flat <- sum(c_flat) / 1000

f_jeff <- function(theta) {return(likelihood(theta) * jeff(theta))}
c_jeff <- c()
for (k in seq(-0.999, 0.999, 0.001)) {c_jeff <- c(c_jeff, f_jeff(k))}
m_jeff <- sum(c_jeff) / 1000

f_tri <- function(theta) {return(likelihood(theta) * tri(theta))}
c_tri <- c()
for (k in seq(-0.999, 0.999, 0.001)) {c_tri <- c(c_tri, f_tri(k))}
m_tri <- sum(c_tri) / 1000


## Plotting the priors
x <- seq(-0.2, 1.0, 0.01)
y_flat <- c()
y_jeff <- c()
y_tri <- c()

for (j in c(1:length(x))) {
  y_flat <- c(y_flat, (likelihood(x[j])*flat)/m_flat)
}

for (j in c(1:length(x))) {
  y_jeff <- c(y_jeff, (likelihood(x[j])*jeff(x[j]))/m_jeff)
}

for (j in c(1:length(x))) {
  y_tri <- c(y_tri, (likelihood(x[j])*tri(x[j]))/m_tri)
}

plot(x, y_flat, type="l", lty=1, col="black", xlim=c(-0.2, 1.0), ylim=c(0.0, 2.6), xlab="theta_hat", ylab="posterior")
par(new=T)
plot(x, y_jeff, type="l", lty=2, col="red", xlim=c(-0.2, 1.0), ylim=c(0.0, 2.6), xlab="", ylab="")
par(new=T)
plot(x, y_tri, type="l", lty=3, col="blue", xlim=c(-0.2, 1.0), ylim=c(0.0, 2.6), xlab="", ylab="")