Commit dbf8f439 authored by oawile's avatar oawile

- add GaA project and give chrmuell rw access to it

parents
function [funVal,F,g1,g2] = EMProb1(x)
n=length(x);
if n~=2
error('wrong dimension')
end
x1=x(1);
x2=x(2);
g1 = 0.124*sqrt(1+x2^2)*(8/x1+1/(x1*x2))-1;
g2 = 0.124*sqrt(1+x2^2)*(8/x1-1/(x1*x2))-1;
penalty = 0;
if g1>0
penalty = n*g1;
end
if g2 > 0
penalty = penalty+n*g2;
end
F = x(1)*sqrt(1+x(2)^2);
funVal = F+penalty;
\ No newline at end of file
% Analyze the sampling data
HaarioData_pi1 = HaarioData_pi1_opt;
HaarioData_pi3 = HaarioData_pi3_opt;
HaarioData_pi4 = HaarioData_pi4_opt;
numIter = 100;
% load HaarioData.mat
burnInLen = 1e3;
% Analysis for Haario's pi1
meanEVec_pi1 = zeros(numIter,1);
for j=1:numIter
temp = mean(squeeze(HaarioData_pi1(burnInLen:end,:,j)));
meanEVec_pi1(j,1) = sqrt(sum(temp.^2,2));
end
meanE_pi1 = mean(meanEVec_pi1)
stdE_pi1 = std(meanEVec_pi1)
conf1 = 0.683;
k1 = sqrt(chi2inv(conf1,N)); % r is the number of dimensions (degrees of freedom)
C1 = k1*eye(N,N);
conf2 = 0.99;
k2 = sqrt(chi2inv(conf2,N)); % r is the number of dimensions (degrees of freedom)
C2 = k2*eye(N,N);
countInside=0;
countOutside=0;
MaxSample=2e4;
inside68Vec_pi1 = zeros(numIter,1);
for i=1:numIter
i
xTrans_pi1 = squeeze(HaarioData_pi1(:,:,i));
xTrans_pi1(:,1) = xTrans_pi1(:,1)./sqrt(100);
plot(xTrans_pi1(:,1),xTrans_pi1(:,2),'r.');
grid on
hold on
error_ellipse(C1(1:2,1:2),zeros(2,1),'style','b')
error_ellipse(C2(1:2,1:2),zeros(2,1),'style','k')
hold off
drawnow;
for j=burnInLen:MaxSample
temp = xTrans_pi1(j,:)*inv(C1)*xTrans_pi1(j,:)';
countInside=countInside+(temp<=k1);
temp = xTrans_pi1(j,:)*inv(C2)*xTrans_pi1(j,:)';
countOutside=countOutside+(temp>k2);
end
countInside
countOutside
inside68Vec_pi1(i,1)=countInside/MaxSample;
outside99Vec_pi1(i,1)=countOutside/MaxSample;
countInside = 0;
countOutside = 0;
end
mean(inside68Vec_pi1-0.683)
std(inside68Vec_pi1-0.683)
mean(outside99Vec_pi1-0.99)
std(outside99Vec_pi1-0.99)
% Analysis for Haario's pi3
meanEVec_pi3 = zeros(numIter,1);
for j=1:numIter
temp = mean(squeeze(HaarioData_pi3(burnInLen:end,:,j)));
meanEVec_pi3(j,1) = sqrt(sum(temp.^2,2));
end
meanE_pi3 = mean(meanEVec_pi3)
stdE_pi3 = std(meanEVec_pi3)
conf1 = 0.683;
k1 = sqrt(chi2inv(conf1,N)); % r is the number of dimensions (degrees of freedom)
C1 = k1*eye(N,N);
conf2 = 0.99;
k2 = sqrt(chi2inv(conf2,N)); % r is the number of dimensions (degrees of freedom)
C2 = k2*eye(N,N);
countInside=0;
countOutside=0;
MaxSample=4e4;
inside68Vec_pi3 = zeros(numIter,1);
for i=1:numIter
i
xTrans_pi3 = squeeze(HaarioData_pi3(:,:,i));
xTrans_pi3(:,2) = xTrans_pi3(:,2) + 0.03*xTrans_pi3(:,1).^2 - 3;
xTrans_pi3(:,1) = xTrans_pi3(:,1)./sqrt(100);
plot(xTrans_pi3(:,1),xTrans_pi3(:,2),'r.');
grid on
hold on
error_ellipse(C1(1:2,1:2),zeros(2,1),'style','b')
error_ellipse(C2(1:2,1:2),zeros(2,1),'style','k')
hold off
drawnow;
for j=burnInLen:MaxSample
temp = xTrans_pi3(j,:)*inv(C1)*xTrans_pi3(j,:)';
countInside=countInside+(temp<=k1);
temp = xTrans_pi3(j,:)*inv(C2)*xTrans_pi3(j,:)';
countOutside=countOutside+(temp>k2);
end
countInside
countOutside
inside68Vec_pi3(i,1)=countInside/MaxSample;
outside99Vec_pi3(i,1)=countOutside/MaxSample;
countInside = 0;
countOutside = 0;
end
mean(inside68Vec_pi3-0.683)
std(inside68Vec_pi3-0.683)
mean(outside99Vec_pi3-0.99)
std(outside99Vec_pi3-0.99)
% Analysis for Haario's pi3 with P=0.1
meanEVec_pi3_01 = zeros(numIter,1);
for j=1:numIter
temp = mean(squeeze(HaarioData_pi3(burnInLen:end,:,j)));
meanEVec_pi3_01(j,1) = sqrt(sum(temp.^2,2));
end
meanE_pi3_01 = mean(meanEVec_pi3_01)
stdE_pi3_01 = std(meanEVec_pi3_01)
conf1 = 0.683;
k1 = sqrt(chi2inv(conf1,N)); % r is the number of dimensions (degrees of freedom)
C1 = k1*eye(N,N);
conf2 = 0.99;
k2 = sqrt(chi2inv(conf2,N)); % r is the number of dimensions (degrees of freedom)
C2 = k2*eye(N,N);
countInside=0;
countOutside=0;
MaxSample=4e4;
inside68Vec_pi3_01 = zeros(numIter,1);
for i=1:numIter
i
xTrans_pi3_01 = squeeze(HaarioData_pi3(:,:,i));
xTrans_pi3_01(:,2) = xTrans_pi3_01(:,2) + 0.03*xTrans_pi3_01(:,1).^2 - 3;
xTrans_pi3_01(:,1) = xTrans_pi3_01(:,1)./sqrt(100);
plot(xTrans_pi3_01(:,1),xTrans_pi3_01(:,2),'r.');
grid on
hold on
error_ellipse(C1(1:2,1:2),zeros(2,1),'style','b')
error_ellipse(C2(1:2,1:2),zeros(2,1),'style','k')
hold off
drawnow;
for j=burnInLen:MaxSample
temp = xTrans_pi3_01(j,:)*inv(C1)*xTrans_pi3_01(j,:)';
countInside=countInside+(temp<=k1);
temp = xTrans_pi3_01(j,:)*inv(C2)*xTrans_pi3_01(j,:)';
countOutside=countOutside+(temp>k2);
end
countInside
countOutside
inside68Vec_pi3_01(i,1)=countInside/MaxSample;
outside99Vec_pi3_01(i,1)=countOutside/MaxSample;
countInside = 0;
countOutside = 0;
end
mean(inside68Vec_pi3_01-0.683)
std(inside68Vec_pi3_01-0.683)
mean(outside99Vec_pi3_01-0.01)
std(outside99Vec_pi3_01-0.01)
% Analysis for Haario's pi4
meanEVec_pi4 = zeros(numIter,1);
for j=1:numIter
temp = mean(squeeze(HaarioData_pi4(burnInLen:end,:,j)));
meanEVec_pi4(j,1) = sqrt(sum(temp.^2,2));
end
meanE_pi4 = mean(meanEVec_pi4)
stdE_pi4 = std(meanEVec_pi4)
conf1 = 0.683;
k1 = sqrt(chi2inv(conf1,N)); % r is the number of dimensions (degrees of freedom)
C1 = k1*eye(N,N);
conf2 = 0.99;
k2 = sqrt(chi2inv(conf2,N)); % r is the number of dimensions (degrees of freedom)
C2 = k2*eye(N,N);
countInside=0;
countOutside=0;
MaxSample=8e4;
inside68Vec_pi4 = zeros(numIter,1);
for i=1:numIter
i
xTrans_pi4 = squeeze(HaarioData_pi4(:,:,i));
xTrans_pi4(:,2) = xTrans_pi4(:,2) + 0.1*xTrans_pi4(:,1).^2 - 10;
xTrans_pi4(:,1) = xTrans_pi4(:,1)./sqrt(100);
plot(xTrans_pi4(:,1),xTrans_pi4(:,2),'r.');
grid on
hold on
error_ellipse(C1(1:2,1:2),zeros(2,1),'style','b')
error_ellipse(C2(1:2,1:2),zeros(2,1),'style','k')
hold off
drawnow;
for j=burnInLen:MaxSample
temp = xTrans_pi4(j,:)*inv(C1)*xTrans_pi4(j,:)';
countInside=countInside+(temp<=k1);
temp = xTrans_pi4(j,:)*inv(C2)*xTrans_pi4(j,:)';
countOutside=countOutside+(temp>k2);
end
countInside
countOutside
inside68Vec_pi4(i,1)=countInside/MaxSample;
outside99Vec_pi4(i,1)=countOutside/MaxSample;
countInside = 0;
countOutside = 0;
end
mean(inside68Vec_pi4-0.683)
std(inside68Vec_pi4-0.683)
mean(outside99Vec_pi4-0.01)
std(outside99Vec_pi4-0.01)
% Analyze the MH sampling data
numIter = 100;
N=8;
% Analysis for MH's pi3
maxLag = 1e3;
trajAutoCorrGaAVec_pi1 = zeros(maxLag+1,N,numIter);
trajAutoCorrGaAVec_pi3 = zeros(maxLag+1,N,numIter);
trajAutoCorrGaAVec_pi4 = zeros(maxLag+1,N,numIter);
burnIn = 1e3;
for i=1:numIter
temp1 = squeeze(HaarioData_pi1(burnIn:end,:,i));
temp3 = squeeze(HaarioData_pi3(burnIn:end,:,i));
temp4 = squeeze(HaarioData_pi4(burnIn:end,:,i));
for d=1:N
% pi1
[ACF1,Lags,Bounds] = autocorr(temp1(:,d),[maxLag],[],[]);
trajAutoCorrGaAVec_pi1(:,d,i) = ACF1;
% pi3
[ACF3,Lags,Bounds] = autocorr(temp3(:,d),[maxLag],[],[]);
trajAutoCorrGaAVec_pi3(:,d,i) = ACF3;
% pi4
[ACF4,Lags,Bounds] = autocorr(temp4(:,d),[maxLag],[],[]);
trajAutoCorrGaAVec_pi4(:,d,i) = ACF4;
end
end
figure(1)
plot(Lags,squeeze(mean(trajAutoCorrGaAVec_pi1,3)),'LineWidth',1.2)
grid on
hold on
ylim([-0.1 1.1])
figure(2)
plot(Lags,squeeze(mean(trajAutoCorrGaAVec_pi3,3)),'LineWidth',1.2)
grid on
hold on
ylim([-0.1 1.1])
figure(3)
plot(Lags,squeeze(mean(trajAutoCorrGaAVec_pi4,3)),'LineWidth',1.2)
grid on
hold on
ylim([-0.1 1.1])
% Analyze the MH sampling data
numIter = 100;
N=8;
burnIn = 1e4;
% Analysis for MH's pi3
maxLag = 1e3;
trajAutoCorrMHVec_pi1 = zeros(maxLag+1,N,numIter);
trajAutoCorrMHVec_pi3 = zeros(maxLag+1,N,numIter);
trajAutoCorrMHVec_pi4 = zeros(maxLag+1,N,numIter);
for i=1:numIter
burnIn = 1e4;
temp1 = squeeze(MHData_pi1(burnIn:end,:,i));
burnIn = 2e4;
temp3 = squeeze(MHData_pi3(burnIn:end,:,i));
burnIn = 4e4;
temp4 = squeeze(MHData_pi4(burnIn:end,:,i));
for d=1:N
% pi1
[ACF1,Lags,Bounds] = autocorr(temp1(:,d),[maxLag],[],[]);
trajAutoCorrMHVec_pi1(:,d,i) = ACF1;
% pi3
[ACF3,Lags,Bounds] = autocorr(temp3(:,d),[maxLag],[],[]);
trajAutoCorrMHVec_pi3(:,d,i) = ACF3;
% pi4
[ACF4,Lags,Bounds] = autocorr(temp4(:,d),[maxLag],[],[]);
trajAutoCorrMHVec_pi4(:,d,i) = ACF4;
end
end
figure(1)
plot(Lags,squeeze(mean(trajAutoCorrMHVec_pi1,3)),'--','LineWidth',1.2)
grid on
hold on
ylim([-0.1 1.1])
%xlabel('Lag k')
%ylabel('Auto-correlation')
figure(2)
plot(Lags,squeeze(mean(trajAutoCorrMHVec_pi3,3)),'--','LineWidth',1.2)
grid on
hold on
ylim([-0.1 1.1])
%xlabel('Lag k')
%ylabel('Auto-correlation')
figure(3)
plot(Lags,squeeze(mean(trajAutoCorrMHVec_pi4,3)))
grid on
hold on
ylim([-0.1 1.1])
% xlabel('Lag k')
% ylabel('Auto-correlation')
% Draw contour lines of distributins pi,pi3 and pi4
C=eye(2,2);
C(1,1) = 100;
conf1 = 0.683;
k1 = chi2inv(conf1,2); % r is the number of dimensions (degrees of freedom)
C1 = k1*C;
conf2 = 0.99;
k2 = chi2inv(conf2,2); % r is the number of dimensions (degrees of freedom)
C2 = k2*C;
[h,x1,y1,z1]=error_ellipse('C',C1,'mu',[0;0],'style','b');
[h,x2,y2,z2]=error_ellipse('C',C2,'mu',[0;0],'style','b');
close gcf
x1_init = x1;
x2_init = x2;
y1_init = y1;
y2_init = y2;
figure(1)
plot(x1,y1,'r','LineWidth',1.1)
hold on
plot(x2,y2,'r','LineWidth',1.1)
grid on
xlim([-51 51])
ylim([-86 16])
set(gca,'XTick',[-50:10:50])
y1 = y1_init - 0.03*x1_init.^2 + 3;
y2 = y2_init - 0.03*x2_init.^2 + 3;
figure(1)
plot(x1,y1,'b','LineWidth',1.1)
hold on
plot(x2,y2,'b','LineWidth',1.1)
grid on
xlim([-51 51])
ylim([-86 16])
set(gca,'XTick',[-50:10:50])
y1 = y1_init - 0.1*x1_init.^2 + 10;
y2 = y2_init - 0.1*x2_init.^2 + 10;
figure(1)
plot(x1,y1,'g','LineWidth',1.1)
hold on
plot(x2,y2,'g','LineWidth',1.1)
grid on
xlim([-51 51])
ylim([-86 16])
set(gca,'XTick',[-50:10:50])
% Draw an example of a random sample for each test case
close all
randInd = 50;%round(rand*100);
x1 = squeeze(HaarioData_pi1(:,1,randInd));
y1 = squeeze(HaarioData_pi1(:,2,randInd));
figure(1)
plot(x1,y1,'r.')
grid on
hold on
xlim([-51 51])
ylim([-86 16])
set(gca,'XTick',[-50:10:50])
x1 = squeeze(HaarioData_pi3(:,1,randInd));
y1 = squeeze(HaarioData_pi3(:,2,randInd));
figure(1)
plot(x1,y1,'b.')
grid on
xlim([-51 51])
ylim([-86 16])
set(gca,'XTick',[-50:10:50])
x1 = squeeze(HaarioData_pi4(:,1,randInd));
y1 = squeeze(HaarioData_pi4(:,2,randInd));
figure(1)
plot(x1,y1,'g.')
grid on
xlim([-51 51])
ylim([-86 16])
set(gca,'XTick',[-50:10:50])
\ No newline at end of file
function h=ellipse(ra,rb,ang,x0,y0,C,Nb)
% Ellipse adds ellipses to the current plot
%
% ELLIPSE(ra,rb,ang,x0,y0) adds an ellipse with semimajor axis of ra,
% a semimajor axis of radius rb, a semimajor axis of ang, centered at
% the point x0,y0.
%
% The length of ra, rb, and ang should be the same.
% If ra is a vector of length L and x0,y0 scalars, L ellipses
% are added at point x0,y0.
% If ra is a scalar and x0,y0 vectors of length M, M ellipse are with the same
% radii are added at the points x0,y0.
% If ra, x0, y0 are vectors of the same length L=M, M ellipses are added.
% If ra is a vector of length L and x0, y0 are vectors of length
% M~=L, L*M ellipses are added, at each point x0,y0, L ellipses of radius ra.
%
% ELLIPSE(ra,rb,ang,x0,y0,C)
% adds ellipses of color C. C may be a string ('r','b',...) or the RGB value.
% If no color is specified, it makes automatic use of the colors specified by
% the axes ColorOrder property. For several circles C may be a vector.
%
% ELLIPSE(ra,rb,ang,x0,y0,C,Nb), Nb specifies the number of points
% used to draw the ellipse. The default value is 300. Nb may be used
% for each ellipse individually.
%
% h=ELLIPSE(...) returns the handles to the ellipses.
%
% as a sample of how ellipse works, the following produces a red ellipse
% tipped up at a 45 deg axis from the x axis
% ellipse(1,2,pi/8,1,1,'r')
%
% note that if ra=rb, ELLIPSE plots a circle
%
% written by D.G. Long, Brigham Young University, based on the
% CIRCLES.m original
% written by Peter Blattner, Institute of Microtechnology, University of
% Neuchatel, Switzerland, blattner@imt.unine.ch
% Check the number of input arguments
if nargin<1,
ra=[];
end;
if nargin<2,
rb=[];
end;
if nargin<3,
ang=[];
end;
%if nargin==1,
% error('Not enough arguments');
%end;
if nargin<5,
x0=[];
y0=[];
end;
if nargin<6,
C=[];
end
if nargin<7,
Nb=[];
end
% set up the default values
if isempty(ra),ra=1;end;
if isempty(rb),rb=1;end;
if isempty(ang),ang=0;end;
if isempty(x0),x0=0;end;
if isempty(y0),y0=0;end;
if isempty(Nb),Nb=300;end;
if isempty(C),C=get(gca,'colororder');end;
% work on the variable sizes
x0=x0(:);
y0=y0(:);
ra=ra(:);
rb=rb(:);
ang=ang(:);
Nb=Nb(:);
if isstr(C),C=C(:);end;
if length(ra)~=length(rb),
error('length(ra)~=length(rb)');
end;
if length(x0)~=length(y0),
error('length(x0)~=length(y0)');
end;
% how many inscribed elllipses are plotted
if length(ra)~=length(x0)
maxk=length(ra)*length(x0);
else
maxk=length(ra);
end;
% drawing loop
for k=1:maxk
if length(x0)==1
xpos=x0;
ypos=y0;
radm=ra(k);
radn=rb(k);
if length(ang)==1
an=ang;
else
an=ang(k);
end;
elseif length(ra)==1
xpos=x0(k);
ypos=y0(k);
radm=ra;
radn=rb;
an=ang;
elseif length(x0)==length(ra)
xpos=x0(k);
ypos=y0(k);
radm=ra(k);
radn=rb(k);
an=ang(k)
else
rada=ra(fix((k-1)/size(x0,1))+1);
radb=rb(fix((k-1)/size(x0,1))+1);
an=ang(fix((k-1)/size(x0,1))+1);
xpos=x0(rem(k-1,size(x0,1))+1);
ypos=y0(rem(k-1,size(y0,1))+1);
end;
co=cos(an);
si=sin(an);
the=linspace(0,2*pi,Nb(rem(k-1,size(Nb,1))+1,:)+1);
% x=radm*cos(the)*co-si*radn*sin(the)+xpos;
% y=radm*cos(the)*si+co*radn*sin(the)+ypos;
h(k)=line(radm*cos(the)*co-si*radn*sin(the)+xpos,radm*cos(the)*si+co*radn*sin(the)+ypos);
set(h(k),'color',C(rem(k-1,size(C,1))+1,:));
end;
function [h,x,y,z]=error_ellipse(varargin)
% ERROR_ELLIPSE - plot an error ellipse, or ellipsoid, defining confidence region
% ERROR_ELLIPSE(C22) - Given a 2x2 covariance matrix, plot the
% associated error ellipse, at the origin. It returns a graphics handle
% of the ellipse that was drawn.
%
% ERROR_ELLIPSE(C33) - Given a 3x3 covariance matrix, plot the
% associated error ellipsoid, at the origin, as well as its projections
% onto the three axes. Returns a vector of 4 graphics handles, for the
% three ellipses (in the X-Y, Y-Z, and Z-X planes, respectively) and for
% the ellipsoid.
%
% ERROR_ELLIPSE(C,MU) - Plot the ellipse, or ellipsoid, centered at MU,
% a vector whose length should match that of C (which is 2x2 or 3x3).
%
% ERROR_ELLIPSE(...,'Property1',Value1,'Name2',Value2,...) sets the
% values of specified properties, including:
% 'C' - Alternate method of specifying the covariance matrix
% 'mu' - Alternate method of specifying the ellipse (-oid) center
% 'conf' - A value betwen 0 and 1 specifying the confidence interval.
% the default is 0.5 which is the 50% error ellipse.
% 'scale' - Allow the plot the be scaled to difference units.
% 'style' - A plotting style used to format ellipses.
% 'clip' - specifies a clipping radius. Portions of the ellipse, -oid,
% outside the radius will not be shown.
%
% NOTES: C must be positive definite for this function to work properly.
default_properties = struct(...
'C', [], ... % The covaraince matrix (required)
'mu', [], ... % Center of ellipse (optional)
'conf', 0.5, ... % Percent confidence/100
'scale', 1, ... % Scale factor, e.g. 1e-3 to plot m as km
'style', '', ... % Plot style
'clip', inf); % Clipping radius
if length(varargin) >= 1 & isnumeric(varargin{1})
default_properties.C = varargin{1};
varargin(1) = [];
end
if length(varargin) >= 1 & isnumeric(varargin{1})
default_properties.mu = varargin{1};
varargin(1) = [];
end
if length(varargin) >= 1 & isnumeric(varargin{1})
default_properties.conf = varargin{1};
varargin(1) = [];
end
if length(varargin) >= 1 & isnumeric(varargin{1})
default_properties.scale = varargin{1};
varargin(1) = [];
end
if length(varargin) >= 1 & ~ischar(varargin{1})
error('Invalid parameter/value pair arguments.')
end
prop = getopt(default_properties, varargin{:});
C = prop.C;
if isempty(prop.mu)
mu = zeros(length(C),1);
else
mu = prop.mu;
end
conf = prop.conf;
scale = prop.scale;
style = prop.style;
if conf <= 0 | conf >= 1
error('conf parameter must be in range 0 to 1, exclusive')
end
[r,c] = size(C);
if r ~= c | (r ~= 2 & r ~= 3)
error(['Don''t know what to do with ',num2str(r),'x',num2str(c),' matrix'])
end
x0=mu(1);
y0=mu(2);