# findind marginal distribution

observedData <- read.csv(file="abcrejection1212overnight_observedData.csv", header=TRUE, sep=",")
observedDatats <- ts(data=observedData$x)

ma2likelihood <- function(theta, y, sigma2=1) {
  H <- rbind(c(1,theta[1],theta[2]))
  A <- rbind(c(0,0,0),
             c(1,0,0),
             c(0,1,0))
  Sigma <- matrix(0,3,3)
  Sigma[1,1] <- sigma2

  # initial states
  P <- sigma2*diag(3)
  x <- rep(0,3)
  suml <- 0
  for (t in 1:length(y)) {
    # likelihood contribution
    resid <- y[t] - H %*% x
    S <- H %*% P %*% t(H)
    suml <- suml + dnorm(resid, 0, sqrt(S), log=TRUE)
    # update step
    K <- P %*% t(H) %*% solve(S)
    x <- x + K %*% resid
    P <- (diag(3) - K %*% H) %*% P
    # prediction step
    P <- A %*% P %*% t(A) + Sigma
    x <- A %*% x
  }
  return(suml)
}

#loglikelihood = apply(X=theta_grid,MARGIN=1,FUN=ma2likelihood,y=observedDatats,sigma2=1)
#normconstant = 0.5*2*4*sum(exp(loglikelihood-max(loglikelihood))*prior_density)


#thetachain <- read.csv(file="abcmcmcresults1312.csv", header=TRUE, sep=",")
#uniquetheta = unique(thetachain)

# plot(density(thetachain[,1]),xlim=c(-2,2))
# abline(v=0.2,col="red",lty="dashed")
# abline(v=0.6,col="green",lty="dashed")
# lines(density(thetachain[,2]))

prior_density = 1/(0.5*2*4)
xseq = seq(from=-2,to=2,length.out = 200)
yseq = seq(from=-1,to=1,length.out = 200)
rect_grid = expand.grid(x = xseq, y= yseq)
theta1 = rect_grid[,1]; theta2 = rect_grid[,2]
theta_grid = rect_grid[(theta1+theta2)>(-1) & (theta1-theta2)<1,]

integrand <- function(t1,t2){ma2likelihood(theta=c(t1,t2),y=observedDatats,sigma2=1)}
marginal_theta1 = sapply(theta_grid[,1],function(x) integrate(f=Vectorize(integrand),
                                                              lower=-1, upper= 1, t1=x))

marginal_theta2 = sapply(theta_grid[,2],function(x) integrate(f=Vectorize(integrand),
                                                              lower=-2, upper= 2, t2=x))


prior_density = 1/(0.5*2*4)


marginal_theta1
mtheta1 = as.numeric(marginal_theta1[1,])
length(theta_grid[,1])


normconstant = 0.5*2*4*sum(exp(mtheta1-max(mtheta1))*prior_density)

plot(theta_grid[,1],exp(mtheta1-max(mtheta1)),pch=20,cex=0.5)
lines(exp(mtheta1-max(mtheta1)))
abline(v=0.6)





mtheta1 = as.numeric(marginal_theta1[1,])
length(theta_grid[,1])

xtheta = sort(theta_grid[,1])

plot(x=theta_grid[,1],y=exp(mtheta1-max(mtheta1)),pch=20)

marggtheta1 = as.numeric(marginal_theta1[1,])
normconstant = 0.5*2*4*sum(exp(marggtheta1-max(marggtheta1))*prior_density)
xes = as.numeric(marginal_theta1[2,])
length(marggtheta1)
x = seq(-2,2,length.out = 1216)

plot(x=x,y=exp(marggtheta1-max(marggtheta1)))



marginal_theta2 = rep(NA,nrow(theta_grid))
for (i in 1:nrow(theta_grid)){
  marginal_theta2[i] = integrate(Vectorize(integrand), lower = -2, upper =2, y= theta_grid[i,2])
}

marginal_theta2 <-lapply(theta_grid[,2],FUN=(integrate(Vectorize(integrand), lower = -2, upper =2)))