function[dummy]=sample_plot(b,b_opt,k)

% This m-file plots samples parameters collected in the array b. b_opt is
% the optimal parameter vectors computed using an optimization method.

nparams  = length(b_opt);
nsamples = prod(size(b))/nparams;

figure(k)
for i = 1:nparams-1
  for l = (i-1)*nparams+1:i*(nparams-1)
    subplot(nparams-1,nparams-1,l)
    j = l-(i-1)*(nparams-1)+1;
    plot(b(j,:),b(i,:),'*')
    hold on, plot(b_opt(j),b_opt(i),'r*'), hold off
    xlabel(['b_',num2str(j)])
    ylabel(['b_',num2str(i)])
  end
end

% Plot the normal quantiles.
figure(k+1)
m = floor(sqrt(nparams)); n = ceil(sqrt(nparams));
if m*n < nparams, m=n, end
for i=1:nparams
  subplot(m,n,i)
    qqplot(b(i,:))
    title(['b_',num2str(i),' normal quantile plot'])
end

% Create normalized histograms of samples and compute confidence intervals.
figure(k+2)
for i=1:nparams
  [nb,tb]=nhist_fn(b(i,:));
  q=quantile(b(i,:),[0.025 0.975]);
  subplot(m,n,i)
    plot(tb,nb,'ok'), hold on
    jj = min( find( abs(tb-q(1))==min(abs(tb-q(1))) ) );
    tb_025 = tb(jj);
    plot(tb_025,0,'*r')
    mb = mean(b(i,:)); 
    plot(mb,0,'*m')
    jj = min( find( abs(tb-q(2))==min(abs(tb-q(2))) ) );
    tb_975 = tb(jj);
    plot(tb_975,0,'*g'), hold off
    title(['b_',num2str(i)])
    if i==1, legend('distribution','.025 quantile','mean','.975 quantile'), end
    fprintf('b_%d      quantile: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',i,tb_025,mb,tb_975)
end      

function[nhist_x,t]=nhist_fn(x)
%
% This function takes a vector x of realizations from a random variable and
% creates a normalized histogram normalized_x with independent variable t. 

[hist_x,t] = hist(x,ceil(length(x)/10));
nhist_x = hist_x/sum(hist_x)/(t(2)-t(1));