• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Dec 3rd, 2019 73 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. %%lab 2
2. %% 1
3. m=100;
4. n=100;
5. sztitan = [215,126];
6. b=sqrt(m*n);
7. sz = [m,n];
8. q1 = @(k)[0 -1 0;-1 4+k^2 -1;0 -1 0];
9. q2 = @(k)[0 -0.1 -10;-0.1 20.4+k^2 -0.1;-10 -0.1 0];
10. q3 = @(k)[0 -1 -2;-1 8+k^2 -1;-2 -1 0];
11. k = 0.01;
12. %imagesc(gmrfprec(sz,q2(k)));
13. Q = gmrfprec(sz,q3(k));     %q3 var bra med 2 eller 3 neighbour niv?er f?r data i labben
14. %% 1.a)
15. s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1);   %covariance for colon, when inverting Q
16. imagesc(reshape(s,b,[]));
17. colorbar
18. %% 1.b)
19. mu = 5;
20. R = chol(Q);
21. x = mu + R \ randn(size(R,1),1);
22. imagesc(reshape(x,b,[]))
23. colorbar
24. %% 1.c) higher order neigbourhoods
25. Q=Q^2;
26. s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1);
27. imagesc(reshape(s,b,[]));
28. %% 2 Interpolation
30. %% Image to restore
31. imagesc(reshape(xmiss,100,[]))
32. %%  Construct the observation matrix based on the location of the known pixels:
33. log_risk = log(Insurance(:,2)./Insurance(:,1));
34. Y=logrisk(known);
35.
36. A = sparse(1:length(Y), find(known), 1,length(Y), numel(xmiss));
37. spy(A)
38.
39. %% Given the observation matrix, a vector of covariates for all grid points,
40. %and precision matrices for x and ?; the joint matrices can be constructed
41. %as
42.
43. Aall = [A ones(length(Y),1)];
44. Qbeta = 1e-6 * speye(1); %use SPEYE to create a SPARSE-matrix
45. Qall = blkdiag(Q, Qbeta);
46.
47. %% c) Compute the posterior precision matrix
48. %assume a very small observation-uncertainty
49. Qeps = 1e5 * speye(length(Y));
50. %posterior precisionx
51. Qxy = Qall + Aall'*Qeps*Aall;
52. spy(Qxy);
53. %% reordering
54. p = amd(Qxy);
55. Qxy = Qxy(p,p);
56. Aall = Aall(:,p);
57. spy(Qxy);
58. %% Compute conditional expectation ... and do inverse reordering
59.
60. % E_xy = inv(Qxy)\Aall'*Qeps*Y;
61. R=chol(Qxy);
62. E_xy = R\(R'\(Aall'*Qeps*Y));
63.
64. E_xy(p) = E_xy;
65. E_zy = [speye(size(Q,1)) ones(size(Q,1),1)]*E_xy;
66. subplot(2,2,1);
67. imagesc(reshape(xmiss,100,[]));
68. colorbar
69. subplot(2,2,2);
70. imagesc(reshape(E_zy,100,[]));
71. colorbar
72. subplot(2,2,3);
73. imagesc(reshape(xtrue,100,[]));
74. colorbar
75. subplot(2,2,4);
76. x_rec = reshape(E_zy,100,[]);
77. imagesc(x_rec-xtrue);
78. colorbar
79. %%
81. imshow(titan);
82. %% reconstruction of image from titan
83. k = 0.001;
84. Q = gmrfprec(sztitan,q1(k));
85. Q=Q*Q;
86. miss=0.9;
88. titan=double(titan)/255;
89. xtrue = titan;
90. known=(rand(size(titan))>miss);
91. dknown=double(known);
92. xmiss=xtrue.*dknown;
93. Y=xmiss(known);
94. A = sparse(1:length(Y), find(known), 1,length(Y), numel(xmiss));
95. Aall = [A ones(length(Y),1)];
96. Qbeta = 1e-6 * speye(1);
97. Qall = blkdiag(Q, Qbeta);
98. %assume a very small observation-uncertainty
99. Qeps = 1e5 * speye(length(Y));
100. %posterior precision
101. Qxy = Qall + Aall'*Qeps*Aall;
102. p = amd(Qxy);
103. Qxy = Qxy(p,p);
104. Aall = Aall(:,p);
105. R=chol(Qxy);
106. E_xy = R\(R'\(Aall'*Qeps*Y));
107. E_xy(p) = E_xy;
108. E_zy = [speye(size(Q,1)) ones(size(Q,1),1)]*E_xy;
109. subplot(2,3,1);
110. % imagesc(reshape(xmiss,215,[]));
111. imshow(reshape(xmiss,215,[]));
112. colorbar
113. subplot(2,3,2);
114. % imagesc(reshape(E_zy,215,[]));
115. imshow(reshape(E_zy,215,[]));
116. colorbar
117. subplot(2,3,3);
118. % imagesc(reshape(xtrue,215,[]));
119. imshow(reshape(xtrue,215,[]));
120. colorbar
121. subplot(2,3,4);
122. x_rec = reshape(E_zy,215,[]);
123. plot(x_rec-xtrue,'xk');
124. subplot(2,3,5)
125. s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1);   %covariance for colon, when inverting Q
126. imagesc(reshape(s,215,[]));
127. colorbar
128.
129.
130.
131. %%
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.
Not a member of Pastebin yet?