% lecture 03
% Moo K. Chung mkchung@wisc.edu
% 
% Two-components Gaussian Mixture Model


% Simulating random numbers
% Y1 ~ N(0, 1), Y2 ~ N(10,5^2)

Y1=normrnd(0,1, [1000 1]);
Y2=normrnd(10,5,[2000 1]);

% We mix Y1 and Y2.
% The known mixing proportions 
% are 1/3 to 2/3. 

Y=[Y1;Y2]
hist(Y,101); set(get(gca,'Children'),'FaceColor',[.8 .8 1])

% The actual density function

f=inline('1/sqrt(2*pi*sigma)*exp(-0.5 * ((x-mu).^2/(sigma.^2)))')

f(0,1,0)
f(0,1,[-1,0,1])

x=-15:0.1:30
mixedf=  1/3*f(0,1,x)+ 2/3*f(10,5,x)
hold on;plot(x,500*mixedf,'r')

% Warning: the theoretical and simulated distributions do not 
% exactly match. Therefore our estimation is likely to be biased.


% Assuming all parameters except the mixing proportion "p" is unknown 
p=[]
p(1)=0.5;
for j=1:10
    p(j+1)= mean((p(j)*f(1,1,Y))./(p(j)*f(1,1,Y) + (1-p(j))*f(10,5,Y)));
end;

figure;plot(p)


% Assuming all parameters are unknown

[mask,mu0,sigma_square,p0]=EMSeg(Y,2)

    
% VBM data for NI 2004 paper
% segmentation is done using SPM'99
% download CCdensity.zip and unzip it

control1=[01 02 05 08 09]';
control2=[10 11 12 13 14 15 16]';

autism1=[01 02 03 04 06 07 08 09]';
autism2=[10 11 13 14 16 17 18]';


dir=''

read_control1=strcat(dir,'CCcontrol0',int2str(control1),'.txt')
read_control2=strcat(dir,'CCcontrol',int2str(control2),'.txt')
read_control=[read_control1; read_control2]

read_autism1=strcat(dir,'CCautism0',int2str(autism1),'.txt')
read_autism2=strcat(dir,'CCautism',int2str(autism2),'.txt')
read_autism=[read_autism1; read_autism2]

n_control=size(read_control,1)
n_autism=size(read_autism,1)

dim2=95;dim3=68;
control=zeros(n_control,dim2,dim3);
autism=zeros(n_autism,dim2,dim3);

for i=1:n_control
    i
    temp=load(read_control(i,:));
    temp=reshape(temp,dim2,dim3);
    %figure;imagesc(temp);colormap('bone');colorbar;
    control(i,:,:)=temp;
end;

% visualization
for i=1:n_autism
    i
    temp=load(read_autism(i,:));
    temp=reshape(temp,dim2,dim3);
    %figure;imagesc(temp);colormap('bone');colorbar;
    autism(i,:,:)=temp;
end;


autism_w=[];
for i=1:n_autism
    autism_w = [autism_w, squeeze(autism(i,:,:))];
end;
figure; imagesc(autism_w);colormap('bone');colorbar;


control_w=[];
for i=1:n_control
    control_w = [control_w, squeeze(control(i,:,:))];
end;
figure; imagesc(control_w);colormap('bone');colorbar;


% smoothing a single image

figure;imagesc(temp);colormap('bone');colorbar;

temps = gaussblur(temp,5)
figure;imagesc(temps);colormap('bone');colorbar;


%HW5: Following Chung et al. NI 2004 paper, perform Gaussian kernel 
%     smoothing and perform a two sample t-test on white matter density. 
%     Discusss your result. There is no need to do multiple comparisions. 

%HW6: Write an EM algorithm for estimating 2 unkonwn parameters mu1 and mu2
%     assuming the mixing proportion p1 and p2, and s.d. sigma1 and
%     sigma2 are known. Need to go through the textbook for detail.

%NOTE: You can submit your homework through email by submitting a single
%      PDF file. Please do not submit multiple files if you are submitting 
%      it by email. Auditing students are not required to submit homeworks.


