Code covered by the BSD License  

Highlights from
Bayesian robust hidden Markov model

image thumbnail

Bayesian robust hidden Markov model

by

 

25 Sep 2013 (Updated )

MatLab object for segmenting sequences of real-valued data with noise, outliers and missing values.

TestBRHMM()
function TestBRHMM()
%TestBRHMM    Test the Bayesian robust hidden Markov model
%   
%   This routine tests the capabilities of the BRHMM on synthetic data. It
%   generates a set of sequences of data from a BRHMM with unknown
%   parameters and, from these data alone, it estimates the parameters of
%   the model responsible for generating them. Results are plotted as a
%   color-coded time series, where each color denotes an assignment to a
%   particular class.
%   
%   Copyright (c) 2013 Gabriel Agamennoni.

% Set number of states, symbols and features.
NumState=2;
NumSym=3;
NumFeat=5;

% Set number of sequences, points per sequence and missing values.
NumSeq=2;
NumPoint=100;
NumMiss=20;

% Set parameter generation options.
TransParam=1/5;
EmissParam=1/5;
LocParam=2;
DispParam=5;

% Set sampling options.
NumDeg=5;
NumObs=1000;

% Print header and display status.
fprintf('\n')
fprintf('Sampling data ... ')

% Generate parameters for sampling.
[Trans,Emiss,Loc,Disp]=GenParam(NumState,NumSym,NumFeat,...
    TransParam,EmissParam,LocParam,DispParam);

% Create model for sampling.
Obj=BRHMM(NumState,NumSym,NumFeat);

% Set hyper-parameters.
Obj.TransWeight=Trans;
Obj.TransStren(:)=NumObs;
Obj.EmissWeight=Emiss;
Obj.EmissStren(:)=NumObs;
Obj.CompLoc=Loc;
Obj.CompScale(:)=NumObs;
Obj.CompDisp=Disp;
Obj.CompPrec(:)=max(NumObs,NumFeat);

% Sample data and remove values at random.
[~,~,~,Feat]=Obj.Sim(NumPoint(ones(NumSeq,1)),'NumDeg',NumDeg);
for i=1:NumSeq
    Feat{i}(randi(NumFeat*NumPoint,NumMiss,1))=nan();
end

% Update status.
fprintf('Done\n')
fprintf('Estimating model ... ')

% Create model for estimation.
Obj=BRHMM(NumState,NumSym,NumFeat);

% Constrain transition parameters.
Obj.TransWeight=Trans;
Obj.TransStren(:)=NumObs;

% Estimate model and state probabilities.
[Obj,Bound,Prob]=Obj.Estim(Feat,'NumDeg',NumDeg);

% Update status.
fprintf('Done\n')
fprintf('Plotting results ... ')

% Close any existing figure.
close('all')

% Plot results.
PlotBound(Bound)
PlotSeg([Feat{:}],[Prob{:}])

% Update status and print footer.
fprintf('Done\n')
fprintf('\n')

end



function [Trans,Emiss,Loc,Disp]=GenParam(NumState,NumSym,NumFeat,...
    TransParam,EmissParam,LocParam,DispParam)

% Sample transition parameters.
Trans=rand(NumState,NumState);
Trans=bsxfun(@rdivide,Trans,sum(Trans,1));
Trans=(1-TransParam)*eye(NumState)+TransParam^2*Trans+...
    (1-TransParam)*TransParam/NumState;

% Sample emission parameters.
Emiss=rand(NumSym,NumState);
Emiss=bsxfun(@rdivide,Emiss,sum(Emiss,1));
Emiss=(1-EmissParam)*eye(NumSym,NumState)+EmissParam^2*Emiss+...
    (1-EmissParam)*EmissParam/NumSym;

% Sample location parameters.
Loc=LocParam*randn(NumFeat,NumSym);

% Allocate space.
Disp=zeros(NumFeat,NumFeat,NumSym);

% Sample dispersion parameters.
for i=1:NumSym
    Fact=randn(NumFeat,NumFeat+DispParam)/sqrt(NumFeat+DispParam);
    Disp(:,:,i)=Fact*Fact';
end

end



function PlotBound(Bound)

% Set options.
FontName='times';
FontSize=20;
LineColor=[0,0,1];
LineWidth=2;
MarkerSize=20;

% Create figure.
Fig=figure(...
    'NumberTitle','off',...
    'Name','Variational Lower Bound on the Model Evidence');

% Create axes.
Axes=axes(...
    'Parent',Fig,...
    'NextPlot','add',...
    'Box','on',...
    'Layer','top',...
    'FontName',FontName,...
    'FontSize',FontSize);

% Annotate axes.
set(get(Axes,'XLabel'),...
    'String','Iteration',...
    'FontName',FontName,...
    'FontSize',FontSize)
set(get(Axes,'YLabel'),...
    'String','Lower bound',...
    'FontName',FontName,...
    'FontSize',FontSize)
set(get(Axes,'Title'),...
    'String','Variational lower bound on the model evidence',...
    'FontName',FontName,...
    'FontSize',FontSize)

% Plot lower bound.
line(...
    'Parent',Axes,...
    'XData',1:numel(Bound),...
    'YData',Bound,...
    'Color',LineColor,...
    'LineWidth',LineWidth,...
    'Marker','.',...
    'MarkerSize',MarkerSize)

% Adjust axes.
set(Axes,...
    'XLim',[0,numel(Bound)+1])

end



function PlotSeg(Feat,Prob)

% Set options.
FontName='times';
FontSize=20;
ColorDilut=1/2;
LineWidth=2;
MarkerSize=20;
Margin=3/20;

% Store number of states, features and points.
[NumFeat,NumPoint]=size(Feat);
[NumState,~]=size(Prob);

% Build color map.
ColorMap=hsv(NumState);

% Create figure.
Fig=figure(...
    'NumberTitle','off',...
    'Name','Bayesian Robust Sequence Segmentation');

% Create invisible axes.
Axes=axes(...
    'Parent',Fig,...
    'Position',[Margin/2,Margin/2,1-Margin,1-Margin],...
    'Visible','off');

% Annotate axes.
set(get(Axes,'Title'),...
    'String','Bayesian robust hidden Markov model',...
    'FontName',FontName,...
    'FontSize',FontSize,...
    'Visible','on')
set(get(Axes,'YLabel'),...
    'String','Feature values',...
    'FontName',FontName,...
    'FontSize',FontSize,...
    'Visible','on')

% Allocate space for axis handles.
Axes=zeros(NumFeat,1);

% Plot features.
for i=1:NumFeat
    
    % Set position.
    Pos(1)=Margin/2;
    Pos(2)=Margin/2+(1-Margin)*(1-i/NumFeat);
    Pos(3)=1-Margin;
    Pos(4)=1/NumFeat-Margin/NumFeat;
    
    % Create axes.
    Axes(i)=axes(...
        'Parent',Fig,...
        'Position',Pos,...
        'NextPlot','add',...
        'Box','on',...
        'Layer','top',...
        'FontName',FontName,...
        'FontSize',FontSize);
    
    % Plot lines. """
    line(...
        'Parent',Axes(i),...
        'XData',1:NumPoint,...
        'YData',Feat(i,:),...
        'Color',ColorDilut*get(Axes(i),'Color'),...
        'LineStyle',':',...
        'LineWidth',LineWidth)
    
    % Plot features.
    for j=1:NumPoint
        line(...
            'Parent',Axes(i),...
            'XData',j,...
            'YData',Feat(i,j),...
            'Color',(1-ColorDilut)*Prob(:,j)'*ColorMap+...
                ColorDilut*get(Axes(i),'Color'),...
            'LineStyle','none',...
            'Marker','.',...
            'MarkerSize',MarkerSize)
    end
    
    % Adjust axes.
    if i<NumFeat
        set(Axes(i),...
            'XTickLabel',{})
    end
    set(Axes(i),...
        'XLim',[0,NumPoint+1])
    
end

% Allocate space for line handles.
Line=zeros(NumState,1);

% Plot pure lines.
for i=1:NumState
    Line(i)=line(...
        'Parent',Axes(1),...
        'XData',[],...
        'YData',[],...
        'Color',ColorMap(i,:),...
        'Marker','.',...
        'MarkerSize',MarkerSize);
end

% Allocate space for labels.
Label=cell(NumState,1);

% Create labels.
for i=1:NumState
    Label{i}=sprintf('State %d: %2.1f%%',i,100*(sum(Prob(i,:))/NumPoint));
end

% Annotate plot.
set(legend(Line,Label{:}),...
    'FontName',FontName,...
    'FontSize',FontSize,...
    'Location','northeast')

% Adjust axis limits.
set(Axes,...
    'XLim',[0,NumPoint+1])

% Link axes.
linkaxes(Axes,'x')

end

Contact us