function marg(M,Sig,SampleSize) %Where M is the mean %Sig is the standard deviation %SampleSize is the sample size %Random Number Generation %Generate norm dist samples with mean M, std Sig x=normrnd(M,Sig,1,SampleSize); x_sum=sum(x); %Sums the values of samples x2_sum=sum(x.^2); %Sums the square of the values of samples u1=[5.01:0.02:15]; %Declares u from 5 to 15 s1=[0.01:0.02:10]; %Declares s from 0 to 10 [u,s]=meshgrid(u1,s1); %Generates square matrix of u and s %Initialise normalisation constants norm_const=0; norm_const_u=0; norm_const_s=0; %To find posterior of mean and standard deviation for i=1:500 for j=1:500 pos(i,j)=(1/((s(i,j)^(SampleSize+1))*sqrt(2*pi)^SampleSize))*exp(-(x2_sum-(2*u(i,j)*x_sum)+SampleSize*u(i,j)^2)/(2*s(i,j)^2)); end; end; norm_const=sum(sum(pos))*0.0004; %To find normalisation constant pos=pos/norm_const; %To normalise area under graph = 1 pos_u=sum(pos,1); %Marginalises std to get posterior of mean norm_const_u=sum(pos_u*0.02); %To find normalisation constant pos_u=pos_u/norm_const_u; %To normalise area under graph = 1 pos_s=sum(pos,2); %Marginalises mean to get posterior of std norm_const_s=sum(pos_s*0.02); %To find normalisation constant pos_s=pos_s/norm_const_s; %To normalise area under graph = 1 subplot(2,2,1); mesh(u,s,pos); %Plots posterior of mean and std Title('Posterior of Mean and Std, P(u,s|x)'); xlabel('Mean'); ylabel('Standard Deviation'); zlabel('Probability Density'); grid on; subplot(2,2,2); contour(u,s,pos); %Contour plot of posterior of mean and std Title('Contour Plot of Posterior of Mean and Std P(u,s|x)'); xlabel('Mean'); ylabel('Standard Deviation'); grid on; subplot(2,2,3); plot(u1,pos_u); %Plots posterior of mean %hold on; postnormu(M,Sig,SampleSize); Title('Posterior of Mean, P(u|x)'); xlabel('Mean'); ylabel('Probability Density'); grid on; subplot(2,2,4); plot(s1,pos_s); %Plots posterior of std %hold on; postnorms(M,Sig,SampleSize); Title('Posterior of Std, P(s|x)'); xlabel('Standard Deviation'); ylabel('Probability Density'); grid on; %To Calculate Max, Mean and Variance max_pos=max(max(pos)); for i=1:500 for j=1:500 if pos(i,j)==max_pos %Locates maximum posterior u_max=u(i,j) %Finds lambda for max posterior s_max=s(i,j) end; end; end; u_mean=sum(sum(pos.*u)*0.0004) s_mean=sum(sum(pos.*s)*0.0004) u_std=sqrt(sum(sum((u.^2).*(pos*0.0004)))-u_mean.^2) s_std=sqrt(sum(sum((s.^2).*(pos*0.0004)))-s_mean.^2)