%% p = posterior(y,xs) %% computes posterior probability p(x|y) given distributions %% in problem 7, for each of the xs given the y. %% function p = posterior(y,xs) ys = repmat(y,1,length(xs)); posterior = exp(-0.5*(ys-xs).^2 - abs(xs));
function posterior = expectedposterior(y,xs) % "mean", or expected posterior ys = repmat(y,1,length(xs)); %% we will need to rescale scalefactor = sum(exp(-0.5*(ys-xs).^2 - abs(xs))); posterior = sum(1/scalefactor * xs.*exp(-0.5*(ys-xs).^2 - abs(xs)));
%% function [MAPs, MMSEs] = prob7(Ys) %% %% returns the best MAP and MMSE estimates for x for at each %% of the different values of Y. function [MAPs, MMSEs] = prob7(Ys) xs = [-5:0.005:20]; ys = []; MAPs = []; MMSEs = []; for Y=Ys, %% for each Y we are given %% compute the posterior distribution for this y ypost = prob7(Y,xs); ys = [ys; ypost/max(ypost)]; %% find the MAP estimate for this curve by taking %% the element ath cooreponds to the "peak" of the %% distribution [m ii] = max(ypost); MAPs = [MAPs xs(ii)]; %% compute the MMSE estimate for this curve by taking %% the expected value of the posterior over all X's MMSEs = [MMSEs expectedprob7(Y,xs)]; end; clf; hold on; plot(xs,ys) ypmaps = []; ypmmses = []; %% plot them, rescaling so that their curves are the same height %% (otherwise, the one at 10 is much much flatter) for ypi =1:length(MAPs), divisor = max( prob7(Ys(ypi),xs) ); ypmaps = [ypmaps prob7(Ys(ypi),MAPs(ypi))/divisor]; ypmmses = [ypmmses prob7(Ys(ypi),MMSEs(ypi))/divisor]; end plot(MAPs, ypmaps, '*'); plot(MMSEs, ypmmses, 'x'); hold off;
y=0.1 | y=1.0 | y=10.0 | |
XMAP | 0 | 0 | 9 |
XMMSE | 0.047 | 0.503 | 9.00 |
We can interpret these results as follows. When our observed value, y, is small, the MAP estimate attributes any deviations from the expecation of our signal, to the noise component. This intuitively 'makes sense', since given our noise model, the most likely cause for small fluctuations is the noise term. Meanwhile, the MMSE, is less eager to completely blame noise, since it is trying to achieve the "smallest error" explanation - and therefore splits its bets between the noise and the signal. For highly improbable situations when the deviation is extremely large, both estimations are equally perplexed, and assign most of the blame to the signal -- since the probability of signal (with an |x| instead of x^2) drops off more slowly than the probability of such large noise.
please see the attached graph. note that each of the "probabilities" are not the actual probabillities -they have been re-scaled to peak at 1.0 within the visible range. Otherwise, the highly y=10 ''peak'' becomes completely flat. We could have used a log-log plot to compensate. Maybe next time. :)