Guest User

Untitled

a guest
Dec 16th, 2019
102
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # findind marginal distribution
  2.  
  3. observedData <- read.csv(file="abcrejection1212overnight_observedData.csv", header=TRUE, sep=",")
  4. observedDatats <- ts(data=observedData$x)
  5.  
  6. ma2likelihood <- function(theta, y, sigma2=1) {
  7. H <- rbind(c(1,theta[1],theta[2]))
  8. A <- rbind(c(0,0,0),
  9. c(1,0,0),
  10. c(0,1,0))
  11. Sigma <- matrix(0,3,3)
  12. Sigma[1,1] <- sigma2
  13.  
  14. # initial states
  15. P <- sigma2*diag(3)
  16. x <- rep(0,3)
  17. suml <- 0
  18. for (t in 1:length(y)) {
  19. # likelihood contribution
  20. resid <- y[t] - H %*% x
  21. S <- H %*% P %*% t(H)
  22. suml <- suml + dnorm(resid, 0, sqrt(S), log=TRUE)
  23. # update step
  24. K <- P %*% t(H) %*% solve(S)
  25. x <- x + K %*% resid
  26. P <- (diag(3) - K %*% H) %*% P
  27. # prediction step
  28. P <- A %*% P %*% t(A) + Sigma
  29. x <- A %*% x
  30. }
  31. return(suml)
  32. }
  33.  
  34. #loglikelihood = apply(X=theta_grid,MARGIN=1,FUN=ma2likelihood,y=observedDatats,sigma2=1)
  35. #normconstant = 0.5*2*4*sum(exp(loglikelihood-max(loglikelihood))*prior_density)
  36.  
  37.  
  38. #thetachain <- read.csv(file="abcmcmcresults1312.csv", header=TRUE, sep=",")
  39. #uniquetheta = unique(thetachain)
  40.  
  41. # plot(density(thetachain[,1]),xlim=c(-2,2))
  42. # abline(v=0.2,col="red",lty="dashed")
  43. # abline(v=0.6,col="green",lty="dashed")
  44. # lines(density(thetachain[,2]))
  45.  
  46. prior_density = 1/(0.5*2*4)
  47. xseq = seq(from=-2,to=2,length.out = 200)
  48. yseq = seq(from=-1,to=1,length.out = 200)
  49. rect_grid = expand.grid(x = xseq, y= yseq)
  50. theta1 = rect_grid[,1]; theta2 = rect_grid[,2]
  51. theta_grid = rect_grid[(theta1+theta2)>(-1) & (theta1-theta2)<1,]
  52.  
  53. integrand <- function(t1,t2){ma2likelihood(theta=c(t1,t2),y=observedDatats,sigma2=1)}
  54. marginal_theta1 = sapply(theta_grid[,1],function(x) integrate(f=Vectorize(integrand),
  55. lower=-1, upper= 1, t1=x))
  56.  
  57. marginal_theta2 = sapply(theta_grid[,2],function(x) integrate(f=Vectorize(integrand),
  58. lower=-2, upper= 2, t2=x))
  59.  
  60.  
  61. prior_density = 1/(0.5*2*4)
  62.  
  63.  
  64. marginal_theta1
  65. mtheta1 = as.numeric(marginal_theta1[1,])
  66. length(theta_grid[,1])
  67.  
  68.  
  69. normconstant = 0.5*2*4*sum(exp(mtheta1-max(mtheta1))*prior_density)
  70.  
  71. plot(theta_grid[,1],exp(mtheta1-max(mtheta1)),pch=20,cex=0.5)
  72. lines(exp(mtheta1-max(mtheta1)))
  73. abline(v=0.6)
  74.  
  75.  
  76.  
  77.  
  78.  
  79. mtheta1 = as.numeric(marginal_theta1[1,])
  80. length(theta_grid[,1])
  81.  
  82. xtheta = sort(theta_grid[,1])
  83.  
  84. plot(x=theta_grid[,1],y=exp(mtheta1-max(mtheta1)),pch=20)
  85.  
  86. marggtheta1 = as.numeric(marginal_theta1[1,])
  87. normconstant = 0.5*2*4*sum(exp(marggtheta1-max(marggtheta1))*prior_density)
  88. xes = as.numeric(marginal_theta1[2,])
  89. length(marggtheta1)
  90. x = seq(-2,2,length.out = 1216)
  91.  
  92. plot(x=x,y=exp(marggtheta1-max(marggtheta1)))
  93.  
  94.  
  95.  
  96. marginal_theta2 = rep(NA,nrow(theta_grid))
  97. for (i in 1:nrow(theta_grid)){
  98. marginal_theta2[i] = integrate(Vectorize(integrand), lower = -2, upper =2, y= theta_grid[i,2])
  99. }
  100.  
  101. marginal_theta2 <-lapply(theta_grid[,2],FUN=(integrate(Vectorize(integrand), lower = -2, upper =2)))
RAW Paste Data