﻿

# 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.
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.
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