SHARE
TWEET

Untitled

a guest Dec 16th, 2019 86 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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top