Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- %%lab 2
- %% 1
- m=100;
- n=100;
- sztitan = [215,126];
- b=sqrt(m*n);
- sz = [m,n];
- q1 = @(k)[0 -1 0;-1 4+k^2 -1;0 -1 0];
- q2 = @(k)[0 -0.1 -10;-0.1 20.4+k^2 -0.1;-10 -0.1 0];
- q3 = @(k)[0 -1 -2;-1 8+k^2 -1;-2 -1 0];
- k = 0.01;
- %imagesc(gmrfprec(sz,q2(k)));
- Q = gmrfprec(sz,q3(k)); %q3 var bra med 2 eller 3 neighbour niv?er f?r data i labben
- %% 1.a)
- s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1); %covariance for colon, when inverting Q
- imagesc(reshape(s,b,[]));
- colorbar
- %% 1.b)
- mu = 5;
- R = chol(Q);
- x = mu + R \ randn(size(R,1),1);
- imagesc(reshape(x,b,[]))
- colorbar
- %% 1.c) higher order neigbourhoods
- Q=Q^2;
- s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1);
- imagesc(reshape(s,b,[]));
- %% 2 Interpolation
- load lab3.mat
- %% Image to restore
- imagesc(reshape(xmiss,100,[]))
- %% Construct the observation matrix based on the location of the known pixels:
- log_risk = log(Insurance(:,2)./Insurance(:,1));
- Y=logrisk(known);
- A = sparse(1:length(Y), find(known), 1,length(Y), numel(xmiss));
- spy(A)
- %% Given the observation matrix, a vector of covariates for all grid points,
- %and precision matrices for x and ?; the joint matrices can be constructed
- %as
- Aall = [A ones(length(Y),1)];
- Qbeta = 1e-6 * speye(1); %use SPEYE to create a SPARSE-matrix
- Qall = blkdiag(Q, Qbeta);
- %% c) Compute the posterior precision matrix
- %assume a very small observation-uncertainty
- Qeps = 1e5 * speye(length(Y));
- %posterior precisionx
- Qxy = Qall + Aall'*Qeps*Aall;
- spy(Qxy);
- %% reordering
- p = amd(Qxy);
- Qxy = Qxy(p,p);
- Aall = Aall(:,p);
- spy(Qxy);
- %% Compute conditional expectation ... and do inverse reordering
- % E_xy = inv(Qxy)\Aall'*Qeps*Y;
- R=chol(Qxy);
- E_xy = R\(R'\(Aall'*Qeps*Y));
- E_xy(p) = E_xy;
- E_zy = [speye(size(Q,1)) ones(size(Q,1),1)]*E_xy;
- subplot(2,2,1);
- imagesc(reshape(xmiss,100,[]));
- colorbar
- subplot(2,2,2);
- imagesc(reshape(E_zy,100,[]));
- colorbar
- subplot(2,2,3);
- imagesc(reshape(xtrue,100,[]));
- colorbar
- subplot(2,2,4);
- x_rec = reshape(E_zy,100,[]);
- imagesc(x_rec-xtrue);
- colorbar
- %%
- titan = imread('titan.jpg');
- imshow(titan);
- %% reconstruction of image from titan
- k = 0.001;
- Q = gmrfprec(sztitan,q1(k));
- Q=Q*Q;
- miss=0.9;
- titan = imread('titan.jpg');
- titan=double(titan)/255;
- xtrue = titan;
- known=(rand(size(titan))>miss);
- dknown=double(known);
- xmiss=xtrue.*dknown;
- Y=xmiss(known);
- A = sparse(1:length(Y), find(known), 1,length(Y), numel(xmiss));
- Aall = [A ones(length(Y),1)];
- Qbeta = 1e-6 * speye(1);
- Qall = blkdiag(Q, Qbeta);
- %assume a very small observation-uncertainty
- Qeps = 1e5 * speye(length(Y));
- %posterior precision
- Qxy = Qall + Aall'*Qeps*Aall;
- p = amd(Qxy);
- Qxy = Qxy(p,p);
- Aall = Aall(:,p);
- R=chol(Qxy);
- E_xy = R\(R'\(Aall'*Qeps*Y));
- E_xy(p) = E_xy;
- E_zy = [speye(size(Q,1)) ones(size(Q,1),1)]*E_xy;
- subplot(2,3,1);
- % imagesc(reshape(xmiss,215,[]));
- imshow(reshape(xmiss,215,[]));
- colorbar
- subplot(2,3,2);
- % imagesc(reshape(E_zy,215,[]));
- imshow(reshape(E_zy,215,[]));
- colorbar
- subplot(2,3,3);
- % imagesc(reshape(xtrue,215,[]));
- imshow(reshape(xtrue,215,[]));
- colorbar
- subplot(2,3,4);
- x_rec = reshape(E_zy,215,[]);
- plot(x_rec-xtrue,'xk');
- subplot(2,3,5)
- s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1); %covariance for colon, when inverting Q
- imagesc(reshape(s,215,[]));
- colorbar
- %%
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement