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)