theta = [0.8, -0.3, 0.4, 1.0;...
         0.0, 0.0, 0.0, 0.0];
R = [0.1, 0.0; 0.0, 1.0];


y = struct('class','RV','names',{{'y'}}); 
y.class = 'RV';
y.names= {'y'};

u = struct('class','RV','names',{{'u'}}); 
%yu = struct('class','RV','names',{{'y','u'}}); 
rgr= struct('class','RV',...
        'names',{{'y', 'y','y','u'}},...
		'times' , [-1, -2, -3, -1]);

%Data generating system
system = struct(...
	'class', 'ArxDS',...
	'y', struct('class','RV','names',{{'y', 'u'}}),...
	'u',  struct('class','RV','names',{{}}),...
	'rgr', rgr,...
	'theta', theta,...
	'offset',  [0.0, 0.0],...
	'r', R,...
	'opt', 'L_theta');

% store results
% logger = struct(...
% 	'class','dirfilelog',...
% 	'dirname','arx_test',...
% 	'maxlen',  1000);

A1.class = 'ARX';
A1.rv = y;
A1.rgr = rgr;
A1.log_level ='logbounds,logevidence';

estimators = { A1};

%experiment description
ndat = 90;
experiment.ndat = ndat;

M=estimator(system,estimators,experiment);%,logger);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot results

subplot(1,2,1);
hold off
plot([M.tth(1:ndat,1:2:end) ],'-.');
title(' Regression parameters \theta')
hold on
plot(M.mean_theta(1:ndat,:));
co = get(gca,'ColorOrder');
for i=1:4
    ind =1:10:ndat;
    h=errorbar(ind,M.mean_theta(ind,i),...
    M.mean_theta(ind,i)-M.lb_theta(ind,i),M.mean_theta(ind,i)-M.ub_theta(ind,i),'.');
    set(h,'color',co(i,:));
end

set(gca,'YLim',[-0.7,2]);

subplot(1,2,2);
hold off
plot([M.rR(1:ndat) ],'-.');
title('Variance parameters r')
hold on
plot(M.mean_r(1:ndat,:));
co = get(gca,'ColorOrder');
ind =1:10:ndat;
for i=1:size(M.mean_r, 2)
    h=errorbar(ind,M.mean_r(ind,i),...
    M.mean_r(ind,i)-M.lb_r(ind,i),M.mean_r(ind,i)-M.ub_r(ind,i),'.');
    set(h,'color',co(i,:));
end