function [] = plot_marginal_is_2(IN,num_bins,varargin)

%PLOT_MARGINAL_IS		Plots marginal ppd's from Gibbs Sampling algorithm
%                       with Importance Sampling correction as defined in
%                       Brooks & Frazer, 2005
%
%USAGE          [Num_IS_out,Num_out] = plot_marginal_is(IN,num_bins,varargin)
%
%INPUT
%       IN: the output structure from Gibbs Sampling routines
%       num_bins: number of bins
%       varargin: optional model bounds for axis scaling
%
%OUTPUT
%       Num_IS_out: corrected number in each bin
%       Num_out: total number in each bin
%
%
%EJEMPLO        plot_marginal_is(FB_cmed_tc10,30,xlabels,'c','r')
%
%BAB 4/6/05

%--------------------------
%CONSTANTS & INITIALIZATION
%--------------------------
num_it=length(IN.HB.M)-IN.IN.GS.ts_num;    %use only  those models samples at Ts
fac=1-(1/IN.IN.GS.ts);                  %the correction factor
M=IN.HB.M(:,end-num_it:end);    %sample only from those models found at Tc
E=IN.HB.E(end-num_it:end);
Num=[];
col_corr='r';   %corrected color string
col_uncorr='k'; %uncorrected color string
num_params=length(IN.IN.OBJ.xlabels);   %number of parameters
sp_m=ceil(num_params/3);                %the number of rows for subplotting

%---------------------------------
%DEFINE GRID FOR MARGINAL PLOTTING
%---------------------------------
MB=IN.IN.OBJ.MB;
if isfield(IN.HB,'Xgrid_mid')==1
    Xgrid_mid=IN.HB.Xgrid_mid;
    Xgrid=IN.HB.Xgrid;
else  
    for i=1:length(MB)
        Xgrid(i,:)=linspace(MB(i,1),MB(i,2),num_bins);  %define bins
        for j=1:length(Xgrid(i,:))-1
            Xgrid_mid(i,j)=(Xgrid(i,j+1)+Xgrid(i,j))/2; %bin midpoints
        end
    end
end

%-------------
%FORM THE PPDs
%-------------
denom = sum(exp(-(fac*E))); %denominator for the correction, this makes the marginals PDFs not just DFs
%LOCATION
loc_i=[1 2 5];
for i=1:length(loc_i)
  plot_smooth(loc_i(i),i,1,Xgrid,Xgrid_mid,M,E,fac,denom,num_it)
end
legend('dx (18.9760:19.4334)','dy ( -155.4694:-155.0002)','d(0:20km)','Location','NorthEast')

%GEOMETRY
geom_i=[3 4 6 7];
for i=1:length(geom_i)
    plot_smooth(geom_i(i),i,2,Xgrid,Xgrid_mid,M,E,fac,denom,num_it)
end
legend('l(0:40km)','w(0:40km)','strike(60:120,240:300)','dip(0:90)','Location','NorthEast')

%SLIP
slip_i=[8 9];
for i=1:length(slip_i)
    plot_smooth(slip_i(i),i,3,Xgrid,Xgrid_mid,M,E,fac,denom,num_it)
end
legend('strike-slip(-2:2m)','dip-slip (-2:2m)','Location','NorthEast')

%***********************
%SUBFUNCTION
%*********************
function []=plot_smooth(i,col_num,plot_num,Xgrid,Xgrid_mid,M,E,fac,denom,num_it)

col_is=['r';'g';'b';'m';'k';'c'];
col_str=['r-';'g-';'b-';'m-';'k-';'c-'];
title_str={'LOCATION','GEOMETRY','SLIP'};

for j=1:length(Xgrid(i,:))-1	%for each bin
    k=find(M(i,:)<Xgrid(i,j+1) & M(i,:)>Xgrid(i,j));	%indices of samples in each bin
    if k>0
        numerator=sum(exp(-fac*E(k)));
        Num_IS(j)=numerator/denom; %THE CORRECTION; number in each corrected bin
        corr_num=numerator/denom;
        Num(j)=length(k);	%number in each bin
    else
        Num_IS(j)=0;
        Num(j)=0;
    end
end

%--------
%PLOTTING
%--------
subplot(1,3,plot_num)
hold on
grid on
set(gca,'Box','on');
Xtmp=(Xgrid_mid(i,:)-min(Xgrid_mid(i,:)))./(max(Xgrid_mid(i,:))-min(Xgrid_mid(i,:)));
%H_is=plot(Xtmp,smooth(Num_IS,20),col_is(col_num));
H_is=stairs(Xtmp,Num_IS,col_is(col_num));
%H=plot(Xtmp,smooth(Num/num_it,20),col_str(col_num,:));
set(H_is,'LineWidth',1.5)
%set(H,'LineWidth',2)
axis([0 1 0 1])     %set the axes

title(title_str(plot_num))
    



