function initials=initialvalue(signalfile,zvaluefile,reglocafile,L)
%%L is number of the components in the normal mixture



signal=load(signalfile);%%%% signals, of the brain region of interest (ROI), gotten by BH procedure 
zvalue=load(zvaluefile);%%%% observed data on the whole image grid including the brain part and the non-brain part
regloca=load(reglocafile);%%%%% location of the ROI on the LX*LY*LZ image grid
 


n_reg=length(regloca);% the number of voxels in the ROI
zvalue_1=zeros(n_reg,1);% to store the z-value of the signals by BH procedure
num=0;
for i=1:n_reg
   pi=regloca(i);
   if signal(i)==1
     num=num+1;
     zvalue_1(num)=zvalue(pi);
     
   end
end
prop=1-num/n_reg;

%%%%%%%%%%%always use set the initial values for beta and h to be zeros
beta=0;
h=0;

initials=[beta,h];


temp=0;



try
  obj=gmdistribution.fit(zvalue_1(1:num,1),L);%% gmdistribution.fit is a Matlab built-in function, see http://www.mathworks.com/help/stats/gmdistribution.fit.html
  for j=1:L
      initials=[initials,obj.mu(j),obj.Sigma(j),obj.PComponents(j)];
  end
  
  
catch err %%%%%%% if gmdistribution.fit is divergent 
  mu=zeros(L,1);

  if L==1
     mu(1)=1;

  elseif mod(L,2)==0 %%%%%%if L is divisible by 2, then let mu_i be [-L,-L/2,-L/4,...,L/4,L/2,L] and thus symmetric around 0.
    temp=L/2;
    for mu_i=1:temp
      mu(mu_i)=mu_i*2;
      mu(L-mu_i+1)=-mu(mu_i);
    end
  else  %%%%%%%if L is not divisible by 2, then let mu_i be [-floor(L/2),-floor(L/2)+1,...,0,floor(L/2)-1,floor(L/2)] and thus symmetric around 0.
    temp=ceil(L/2);
      for mu_i=1:L
        mu(mu_i)=mu_i-temp;
      end

  end



  Sigma=ones(L,1);
  PComponents=(1/L)*ones(L,1);
  
  for j=1:L
      initials=[initials,mu(j),Sigma(j),PComponents(j)];
  end

end