function colors = bsizerC2(Y,T,T_vec,nu,W,a_l0,b_l0,Psi,meth,alpha,Sflag,l_vec,nsB,nsM,fname,rs)
%
% function MAP = bsizerC2(Y,T,T_vec,nu,W,a_l0,b_l0,Psi,meth,alpha,Sflag,l_vec,nsB,nsM,fname,rs)
%
% Bayesian Sizer in the Continuous case with dependent observations and errors in predictors.
%                 
% 
% INPUTS:
% Y     - data points (n*1 -vector)
% T     - measure points, e.g. time points (n*1 -vector)
% T_vec - points at where differences are calculated. Can be the same as T. (n*1 -vector)
% nu    - scalar parameter nu
% W     - Matrix prior parameter for Sigma
% a_l0  - prior shape parameter for lambda_0
% b_l0  - prior scale parameter for lambda_0
% b_l0  - prior scale parameter for lambda_0
% Psi   - vector prior parameter of standard deviations of t_i's (n*1 -vector).
%         Set very small values if time points are treated fixed.
% meth  - 0 independent maps, 
%         1 conservative simultaneous, 
%         2 simultaneous using highest pointwise probabilities
%         [OPTIONAL: default 0]
% alpha - Confidence level (0.5 < alpha < 1)
%         [OPTIONAL: default 0.8]
% Sflag - 0 - binary coloring 
%         1 - color maps with shades (works only with independent and MP maps)
%         [OPTIONAL: default 0]          
% l_vec - column vector containing the smoothing parameters in increasing order
%         [OPTIONAL: default logspace(log10(.5),6,80)]
% nsB   - Length of burnin in MCMC simulation
%         [OPTIONAL: default 1000]
% nsM   - Length of monitor phase in MCMC simulation
%         [OPTIONAL: default 3000]
%
% fname - Name of the file, where the generations of mu_lambda-sampler are stored
%
% rs    - 0 do not generate new mu-lambda sample (read it from the file 'fname') 
%         1 sample new (mu,lambda,Sigma,tau)-sample
%         [OPTIONAL: default 1]
% OUTPUT: Bayesian Sizer map (the color matrix)
%
% The routine uses Statistics toolbox (to generate random numbers from Gamma, multivariate
% normal and Inverse -Wishart distributions)
%


[rY,cY]=size(Y);
[rT,cT]=size(T);
T=T';

error(nargchk(8,16,nargin))

if cY>1 | cT>1
  error('Y and T must be column vectors')
end;

if rY~=rT
  error('Y and T must be of same length')
end;

if nu<=rY 
   disp('Note: Parameter nu must be larger than n (to get a proper prior)')
end;

%if nu-round(nu)~=0
%  error('nu must be an integer')
%end;

[rW,cW]=size(W);
if rW~=cW | rW~=rY
  error('Matrix W has to be a square matrix of size n*n')
end;

if det(W)< 0 
  error('Matrix W has to be positively (semi)definite')
end;

if a_l0<=0 
  error('Shape parameter a_l0 must be positive')
end;

if b_l0<=0 
  error('Scale parameter b_l0 must be positive')
end;

[rPsi,cPsi]=size(Psi);
if cPsi>1 | rPsi~=rY
  error('Psi must be a column vector of size n*1')
end;

if min(Psi)<0
  error('Every component of Psi must be non negative')
end;  

if nargin <= 8 ;    %  
  meth_u = 0 ;
else ;              %  meth was specified, use that
  if length(meth)>0
     meth_u = meth ;
  else
    meth_u = 0 ;
  end;  
end ;
meth=meth_u;
clear meth_u;

if ~(ismember(meth,[0 1 2 3]))
 error('meth must be either 0, 1, 2 or 3') 
end;

if nargin <= 9;    
  alpha_u = 0.80 ;
else ;              %  alpha was specified, use that
  if length(alpha)>0
    alpha_u = alpha ;
  else
    alpha_u = 0.8;
  end; 
end ;
alpha=alpha_u;
clear alpha_u;

if ~(alpha>0.5 & alpha<1)
  error('alpha must be between 0.5 and 1')
end;

if nargin <= 10 ;    
  Sflag_u = 0;    
else ;              %  Sflag was specified, use that
  Sflag_u = Sflag ;
end ;
Sflag=Sflag_u;
clear Sflag_u;

if ~(ismember(Sflag,[0 1]))
 error('Sflag must be either 0 or 1') 
end;

if nargin <= 11 ;    
  l_vec_u = [0 logspace(log10(.5),6,80)];    
else               %  l_vec was specified, use that
  if length(l_vec)>0  % The input was not empty matrix
     l_vec_u = l_vec ;
  else
     l_vec_u = [0 logspace(log10(.5),6,80)];
  end;  
end ;
l_vec=l_vec_u;
clear l_vec_u;

[sv,si]=sort(l_vec);
if min(l_vec)<0 | sum(abs(si-[1:length(l_vec)]))>1
  error('elements of l_vec should be non negative and in increasing order')
end;

if nargin <= 12 ;      
  nsB_u = 1000 ;  
else ;              %  nsB was specified, use that
  nsB_u = nsB ;
end ;
nsB=nsB_u;
clear nsB_u;

if nsB<0 | (nsB-ceil(nsB)>0)
  error('length of burnin phase must be non negative integer')
end;  

if nargin <= 13 ;      
  nsM_u = 1000 ;  
else ;              %  nsM was specified, use that
  nsM_u = nsM ;
end ;
nsM=nsM_u;
clear nsM_u;

if nsM<=100 | (nsM-ceil(nsM)>0)
  error('length of monitor phase must be positive integer and larger than 100')
end;  

if nargin <= 14 ;      
  fname_u = ['data_run_cont_1'] ;  
else ;              %  nsM was specified, use that
  fname_u = fname ;
end ;
fname=fname_u;
clear fname_u;

if nargin <= 15 ;      
  rs_u = 1 ;  
else ;              %  rs was specified, use that
  rs_u = rs ;
end ;
rs=rs_u;
clear rs_u;

if ~(ismember(rs,[0 1]))
 error('rs must be either 0 or 1') 
end;

% Lets us first define some variables that are needed later
Psi=diag(Psi.^2); % Note: components of Psi were standard deviations, so square them.
t_prop_sigma=Psi; % Proposal distribution of \tau's. Edit this if in doubt!
pinv_Psi=pinv(Psi);
beta_t=0;
n_lvalues=length(l_vec);
n=length(Y);
%D=diag(-1*ones(n,1))+diag(ones(n-1,1),1);
%D=D(1:n-1,:);
%D(n,n-1:n)=[-1 1];
%D2=2*diag(ones(n,1))+diag(ones(n-1,1),1);
%D2(n,n-1:n)=[-1 1];
%T_lag=diag([diff(T) T(end)-T(end-1)]);
ce=5; % Exponent used in color mapping (in case of maps with color shades)
nT=length(T_vec);


% Background coloring of the maps
colors(:,:,1)=.29*ones(n_lvalues,nT); % Initialize red channel   %.5
colors(:,:,2)=.42*ones(n_lvalues,nT); % Initialize green channel %.2 
colors(:,:,3)=.59*ones(n_lvalues,nT); % Initialize blue channel  %.5


% Set initial values for MCMC sampler
lambda_0=a_l0/b_l0;
wish_cov=W;
Sigma=inv(wishrnd(inv(wish_cov),nu));
t_old=T';
K1=diag(diff(t_old(2:end)).^(-1));
K2=diag(diff(t_old(1:(end-1))).^(-1));
A=diag(-1*ones(n-3,1),1)+diag(ones(n-4,1),2);
A=[A' ; [zeros(1,n-4) 1 -1]; [zeros(1,n-3) 1] ]'; 
B=diag(-1*ones(n-2,1))+ diag(ones(n-3,1),1);
B=[B' ; [zeros(1,n-3)  1]; [zeros(1,n-2) ] ]';
C=(K1*A-K2*B);
h_vec=diff(t_old);
R=diag(1/3.*(h_vec(1:(n-2))+h_vec(2:(n-1))))...
  +diag(1/6.*h_vec(2:(n-2)),1)+diag(1/6.*h_vec(2:(n-2)),-1);
K=C'*pinv(R)*C;
mu=pinv(eye(n)+lambda_0*Sigma*K)*Y;


if rs==1

disp('     Generating samples, please wait ....')
 for s=1:(nsM+nsB)

  if (mod(s,(nsM+nsB)/10)==0) 
    fprintf(1,'%2.0f percent ready \n',(s/(nsM+nsB)*100)); 
  end;
  wish_cov=W+(Y-mu)*(Y-mu)';
  %wish_cov=0.5*(wish_cov+wish_cov');  % Guarantee that the wish_cov is symmetric 
  Sigma=inv(wishrnd((inv(wish_cov)+inv(wish_cov)')/2,nu+1));
  %if (s>nsB) Sigma_sample(s-nsB,:,:)=Sigma;  end;
  lambda_0=gamrnd(a_l0+(n-2)/2,1/(b_l0+1/2*mu'*K*mu));
  %if (s>nsB) lambda_0_sample(s-nsB)=lambda_0;  end;
  mu=mvnrnd((pinv(eye(n)+lambda_0*Sigma*K)*Y)', pinv(pinv(Sigma)+lambda_0*K),1)';
  %if (s>nsB) mu_sample(s-nsB,:)=mu';  end;

  
  % Update the \tau-vector with Metropolis-Hastings. 
  
  t_prop=mvnrnd(t_old,t_prop_sigma,1)';
  % The commented end of the next line in for case, when you trust your interval on \tau's
  if (all(diff(t_prop)>0)) %&((min(t_prop)>=T(1))&(max(t_prop)<=T(end))))
    K1_prop=diag(diff(t_prop(2:end)).^(-1));
    K2_prop=diag(diff(t_prop(1:(end-1))).^(-1));
    C_prop=(K1_prop*A-K2_prop*B);
    h_vec_prop=diff(t_prop); 
    R_prop=diag(1/3.*(h_vec_prop(1:(n-2))+h_vec_prop(2:(n-1))))...
	+diag(1/6.*h_vec_prop(2:(n-2)),1)+diag(1/6.*h_vec_prop(2:(n-2)),-1);
    K_prop=C_prop'*pinv(R_prop)*C_prop;
    prop=-.5*(t_prop-T')'*pinv_Psi*(t_prop-T')-.5*lambda_0*mu'*K_prop*mu;
    old=-.5*(t_old -T')'*pinv_Psi*(t_old -T') -.5*lambda_0*mu'*K*mu;
   
    ratio=min(1,exp(prop-old));
    if rand<=ratio
      %disp('accepted')
      t_old=t_prop;
      if (s>nsB)
	%T_sample(s-nsB,:)=t_old';
      end;	
      K=K_prop;
      C=C_prop;
      beta_t=beta_t+1;
      R=R_prop;
    else
      if (s>nsB)
	%T_sample(s-nsB,:)=t_old';
      end;
    end;
  else
      if (s>nsB)
         %T_sample(s-nsB,:)=t_old'; 
      end;	
  end;
    
  if (s>nsB)
   temp=zeros(n,n);
   temp(2:(end-1),(2:end-1))=pinv(R);
   for smth_par=1:n_lvalues
     smooth=pinv(eye(n)+l_vec(smth_par)*K)*mu;
     gamma=temp*[zeros(1,n);C;zeros(1,n)]*smooth;
     cntr=1;
     for t=1:length(T_vec)
       if ~(T_vec(t)>=t_old(cntr)&T_vec(t)<=t_old(cntr+1))
	 cntr=cntr+1;
       end;	 
       if (T_vec(t)>=t_old(cntr)&T_vec(t)<=t_old(cntr+1)) 
	 h_c=t_old(cntr+1)-t_old(cntr);
	 sample_0(t,s-nsB,smth_par)=(smooth(cntr+1)-smooth(cntr))/h_c...
	                          -1/6*(t_old(cntr+1)+t_old(cntr)-2*T_vec(t))...
	                          *((1+(T_vec(t)-t_old(cntr))/h_c)*gamma(cntr+1)...
				  +((1+(t_old(cntr+1)-T_vec(t))/h_c)*gamma(cntr)))...
				  -1/6*(T_vec(t)-t_old(cntr))*(t_old(cntr+1)-T_vec(t))...
				  *(gamma(cntr+1)-gamma(cntr))/h_c;
       end;     
     end;
     %keyboard;
   end;
  end;
 
 
 
 end;  
 fprintf(1,'acceptance ratio for=%f\n',beta_t/(nsM+nsB));
disp('     Sample generation finished. Now calculating credible regions with different values of lambda. Please wait ....')


save(fname,'sample_0');
% Uncomment the lines below and above to investigate the samples of \Sigma, \mu and \tau.
%save Sigma_sample Sigma_sample  
%save mu_sample mu_sample
%save lambda_0_sample lambda_0_sample
%if beta_t>0
 %save T_sample T_sample
%end;

else
  load(fname);
end;

%wb = waitbar(0,'Progress');

 for smth_par=1:n_lvalues    % loop over smoothing parameters 

    disp( smth_par /n_lvalues )
    %waitbar( smth_par /n_lvalues);
     
    
    p_vec=zeros(nT,1);
    sample=sample_0(:,:,smth_par);
    mean_dif=mean(sample');
    
    p_vec(find(mean_dif<0))=sum(sample(find(mean_dif<0),:)'<0)./nsM;
    p_vec(find(mean_dif>0))=sum(sample(find(mean_dif>0),:)'>0)./nsM;
    p_vec=p_vec.*(p_vec>=alpha);
    
  if meth==0                  % Independent coloring scheme

       %%%%%%%%%%%% Independent coloring scheme %%%%%%%%%%%%%%%  
      
     if Sflag==0  % The binary coloring using only two colors
       
       for r=1:nT            % loop over measure points
          if (mean_dif(r)>0 & p_vec(r)>=alpha)
	       colors(smth_par,r,1)=0;
	       colors(smth_par,r,2)=0;
	       colors(smth_par,r,3)=1;
 	     end;
	    if (mean_dif(r)<0 & p_vec(r)>=alpha)
               colors(smth_par,r,1)=1;
	       colors(smth_par,r,2)=0;
	       colors(smth_par,r,3)=0;
	    end;  
       end;	
       
      else  % Coloring for maps using blue and red palette
       
        for r=1:nT            % loop over measure points
          if (mean_dif(r)>0 & p_vec(r)>=alpha)
	        colors(smth_par,r,1)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,3)=1;
 	     end;
	     if (mean_dif(r)<0 & p_vec(r)>=alpha)
	       colors(smth_par,r,1)=1; 
	       colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	       colors(smth_par,r,3)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	     end;  
       end;	

     end;   
       
   end;   
  
   if meth==1           % Simultaneous credible intervals I (conservative version)
   
   %%%%%%%%%%%%%% Simultaneous inference in \tau (conservative version) %%%%%%%%%%%%%%%
      
	 mean_dif=mean_dif';
	 std_dif=sqrt(diag(cov(sample')));
	 stded_sample=abs(((sample-repmat(mean_dif,1,nsM)))./repmat(std_dif,1,nsM));
         d=[floor(min(max(stded_sample'))*10)/10:.025:floor(min(max(stded_sample'))*10)/10+200*.025];
         for d_ind=1:length(d);
	   if (sum(min(stded_sample<d(d_ind)))/nsM)>alpha
	      break
	   end;
	 end; 
	  
	 
	 for r=1:nT         % loop over measure points
	   
	   if (mean_dif(r)-d(d_ind)*std_dif(r)>0)
	      if Sflag==0
	        colors(smth_par,r,1)=0;
	        colors(smth_par,r,2)=0;
	        colors(smth_par,r,3)=1;    
	      else       
                colors(smth_par,r,1)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,3)=1;
	      end; 
	   end;
	   if (mean_dif(r)+d(d_ind)*std_dif(r)<0)
	      if Sflag==0
   	  	colors(smth_par,r,1)=1;
	        colors(smth_par,r,2)=0;
	        colors(smth_par,r,3)=0;
	      else	 
  		colors(smth_par,r,1)=1; 
	        colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,3)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	      end;	 
	   end;  
     	   
         end;	 
	 
    end;	 % meth==1

    
   if meth==2              % Simultaneous credible intervals in \tau

    %%%%%%%%%%%%%%%%%%%%%%% Highest pointwise probability method  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
    [b,k]=sort(-p_vec);
    for r=1:nT     
             
   	     cum_sum=0;
	     ind_vec=zeros(length(k(1:r)),1);          % Initialize the event
             ind_vec=mean_dif(k(1:r))>0;               % Define the event; positive slopes are coded with ones, 
	                                               % negative with zeros
	     for (s_ind=1:nsM)
	       cum_sum=cum_sum+isequal(sample(k(1:r),s_ind)>0,ind_vec'); % Count the cumulative sum of favourable cases 
	                                                                  % (event takes place)
         end;
	   
	  if ((cum_sum/nsM)<=alpha)
         break;                     
	  
	  
	  else
	    
	     if (p_vec(k(r)) >= alpha & mean_dif(k(r)) > 0 )           
	        if Sflag==0
		  colors(smth_par,k(r),1)=0;
	          colors(smth_par,k(r),2)=0;
	          colors(smth_par,k(r),3)=1;
	        else
 	          colors(smth_par,k(r),1)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(smth_par,k(r),2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(smth_par,k(r),3)=1;
	        end;	 
          end;

         if (p_vec(k(r)) >=alpha & mean_dif(k(r)) < 0 )         
           if Sflag==0
		 colors(smth_par,k(r),1)=1;
	         colors(smth_par,k(r),2)=0;
	         colors(smth_par,k(r),3)=0;
	       else
	         colors(smth_par,k(r),1)=1; 
	         colors(smth_par,k(r),2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	         colors(smth_par,k(r),3)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	       end; 
         end;    
	 end; % ((cum_sum/ns)<=alpha) 
	
   end;                  % for r=1:nT
   
  end;        % of method 2  
      
 end; % smoothing parameter

 %%%%%%%%%%%%%%%%%%%%%%%%% Twice simultaneous features %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
  if meth==3
  
    p_vec=zeros(nT*n_lvalues,1);

    for smth_par=1:n_lvalues    % loop over smoothing parameters 
      sample=sample_0(:,:,smth_par);
      mean_dif=mean(sample');

      p_vec((smth_par-1)*nT+find(mean_dif<0))=sum(sample(find(mean_dif<0),:)'<0)./nsM;
      p_vec((smth_par-1)*nT+find(mean_dif>0))=sum(sample(find(mean_dif>0),:)'>0)./nsM;
      p_vec=p_vec.*(p_vec>=alpha);
      
      sample_vec(((smth_par-1)*nT+1):((smth_par)*nT),:)=sample;
      mean_dif_vec(((smth_par-1)*nT+1):((smth_par)*nT))=mean(sample');
      
    end;
   
    [b,k]=sort(-p_vec);
      
    for r=1:(nT*n_lvalues)
             
   	  cum_sum=0;
	  ind_vec=zeros(length(k(1:r)),1);            % Initialize the event
          ind_vec=mean_dif_vec(k(1:r))>0;             % Define the event; positive slopes are coded with ones, 
	                                              % negative with zeros
						       
	  for (s_ind=1:nsM)
              cum_sum=cum_sum+isequal(sample_vec(k(1:r),s_ind)>0,ind_vec'); 
                                            	      % Count the cumulative sum of favourable cases 
	                                              % (event takes place)
          end;

        if ((cum_sum/nsM)<=alpha)
            break;                     
        else
	    if (p_vec(k(r)) >= alpha & mean_dif_vec(k(r)) > 0 )           
	        if Sflag==0
		  colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,1)=0;
	          colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,2)=0;
	          colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,3)=1;
	        else
		  colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,1)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,3)=1;
	        end;	 
          end;

	  if (p_vec(k(r)) >=alpha & mean_dif_vec(k(r)) < 0 )         
            if Sflag==0
		 colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,1)=1;
	         colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,2)=0;
	         colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,3)=0;
	    else
	         colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,1)=1; 
	         colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	         colors(floor((k(r)-1)/nT)+1,mod(k(r)-1,nT)+1,3)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	    end; 
          end;    

       end; % ((cum_sum/ns)<=alpha) 
	
   end;                 
      
 
end;  % meth==3
 
 % Drawing the output

if Sflag==1

  col_inf_R(:,:,1)=ones(1,101);
  col_inf_R(:,:,2)=[0:.01:1].^ce*.6;
  col_inf_R(:,:,3)=[0:.01:1].^ce*.6;
  col_inf_B(:,:,1)=[0:.01:1].^ce*.6;
  col_inf_B(:,:,2)=[0:.01:1].^ce*.6;
  col_inf_B(:,:,3)=ones(1,101);
  
  subplot('Position',[.1 .25 .8 .65])
  image(T_vec,[log10(min(l_vec)) log10(max(l_vec))],colors)
  set(gca,'YDir','normal')
  ylabel('log_{10}(\lambda)','FontSize',13);
  %xlabel('Data index')
  if (meth==0)
     title(['Pointwise credible intervals. \alpha= ',num2str(alpha)])
  end;
  if (meth==1)
     title(['Simultaneous credible intervals (conservative method). \alpha= ',num2str(alpha)])
  end;
  if (meth==2)
     title(['Simultaneous credible intervals. \alpha= ',num2str(alpha)])
  end;  
  if (meth==3)
     title(['Twice simultaneous credible intervals (\tau,\lambda). \alpha= ',num2str(alpha)])
  end;  

  subplot('Position',[.52 .07 .37 .03])
  image(col_inf_R);
  set(gca,'YTick',[])
  ticks=[1 20 40 60 80 101];
  set(gca,'XTick',ticks);
  set(gca,'XTickLabel',[alpha+[0:.2:1].*(1-alpha)])
  
  subplot('Position',[.1 .07 .37 .03])
  image(col_inf_B);
  set(gca,'YTick',[])
  set(gca,'XTick',ticks);
  set(gca,'XTickLabel',[alpha+[0:.2:1].*(1-alpha)])
  
else
  
  image(T_vec,[log10(min(l_vec)) log10(max(l_vec))],colors)
  set(gca,'YDir','normal')
  ylabel('log_{10}(\lambda)','FontSize',13);
  %xlabel('Data index')
   if (meth==0)
     title(['Pointwise credible intervals. \alpha= ',num2str(alpha)])
  end;
  if (meth==1)
     title(['Simultaneous credible intervals (conservative method). \alpha= ',num2str(alpha)])
  end;
  if (meth==2)
     title(['Simultaneous credible intervals. \alpha= ',num2str(alpha)])
  end;  
  if (meth==3)
     title(['Twice simultaneous credible intervals (\tau,\lambda). \alpha= ',num2str(alpha)])
  end;  

end;
 
 
%close(wb);
disp('     Finished!')