《Computational Statistics with Matlab》硬译2

T=500;
sigma=1;
thetamin=-30;thetamax=30;
theta=zeros(1,T);
seed=1;rand('state',seed);randn('state',seed);
theta(1)=unifrnd(thetamin,thetamax);
t=1;

while t<T
t=t+1;
theta_star=normrnd(theta(t-1),sigma);
alpha=min([1 cauchy(theta_star)/cauchy(theta(t-1))]);
u=rand;
if u<alpha
theta(t)=theta_star;
else
theta(t)=theta(t-1);
end
end

figure(1);clf;
subplot(3,1,1);
nbins=200;
thetabins=linspace(thetamin,thetamax,nbins);
counts=hist(theta,thetabins);
bar(thetabins,counts/sum(counts),'k');
xlim([thetamin thetamax]);
xlabel('	heta');ylabel('p(	heta)');
y=cauchy(thetabins);
hold on;
plot(thetabins,y/sum(y),'r--','LineWidth',3);
set(gca,'YTick',[]);
subplot(3,1,2:3);
stairs(theta,1:T,'k-');
ylabel('t');xlabel('	heta');
set(gca,'YDir','reverse');
xlim([thetamin thetamax]);

function y=bivexp(theta1,theta2)
lambda1=0.5;
lambda2=0.1;
lambda=0.01;
maxval=8;
y=exp(-(lambda1+lambda)*theta1-(lambda2+lambda)*theta2-lambda*maxval);

T=5000;
thetamin=[0 0];
thetamax=[8 8];
seed=1;
rand('state',seed);
randn('state',seed);
theta=zeros(2,T);
theta(1,1)=unifrnd(thetamin(1),thetamax(1));
theta(2,1)=unifrnd(thetamin(2),thetamax(2));

t=1;
while t<T
t=t+1;
theta_star=unifrnd(thetamin,thetamax);
pratio=bivexp(theta_star(1),theta_star(2))/bivexp(theta(1,t-1),theta(2,t-1));
alpha=min([1 pratio]);
u=rand;
if u<alpha
theta(:,t)=theta_star;
else
theta(:,t)=theta(:,t-1);
end
end

figure(1);clf;
subplot(1,2,1);
nbins=10;
thetabins1=linspace(thetamin(1),thetamax(1),nbins);
thetabins2=linspace(thetamin(2),thetamax(2),nbins);
hist3(theta','Edges',{thetabins1 thetabins2});
thetabins1=linspace(thetamin(1),thetamax(1),nbins);

xlabel('	heta_1');ylabel('	heta_2');zlabel('counts');
az=61;el=30;
view(az,el);
subplot(1,2,2);
nbins=20;
thetabins1=linspace(thetamin(1),thetamax(1),nbins);
thetabins2=linspace(thetamin(2),thetamax(2),nbins);
[theta1grid,theta2grid]=meshgrid(thetabins1,thetabins2);
ygrid=bivexp(theta1grid,theta2grid);
mesh(theta1grid,theta2grid,ygrid);
xlabel('	heta_1');ylabel('	heta_2');
zlabel('f(	heta_1,	heta_2)');
view(az,el);

function tau=kendalltau(order1,order2)
[dummy,ranking1]=sort(order1(:)',2,'ascend');
[dummy,ranking2]=sort(order2(:)',2,'ascend');
N=length(ranking1);
[ii,jj]=meshgrid(1:N,1:N);
ok=find(jj(:)>ii(:));
ii=ii(ok);
jj=jj(ok);
nok=length(ok);
sign1=ranking1(jj)>ranking1(ii);
sign2=ranking2(jj)>ranking2(ii);
tau=sum(sign1~=sign2);

lambda = 0.1; % scaling parameter
labels = { 'Washington' , 'Adams' , 'Jefferson' , 'Madison' , 'Monroe' };
omega = [ 1 2 3 4 5 ]; % correct ordering
L = length( omega ); % number of items in ordering
T = 500; % Set the maximum number of iterations
seed=1; rand( 'state' , seed ); randn('state',seed ); % set the random seed
theta = zeros( L , T ); % Init storage space for our samples
theta(:,1) = randperm( L ); % Random ordering to start with
t = 1;
while t < T % Iterate until we have T samples
t = t + 1;
lasttheta = theta(:,t-1); % Get the last theta
whswap = randperm( L ); whswap = whswap(1:2);
theta_star= lasttheta;
theta_star( whswap(1)) = lasttheta( whswap(2));
theta_star( whswap(2)) = lasttheta( whswap(1));
dist1 = kendalltau( theta_star , omega );
dist2 = kendalltau( lasttheta , omega );
pratio = exp(-dist1*lambda) / exp(-dist2*lambda);
alpha = min( [ 1 pratio ] );
u = rand; % Draw a uniform deviate from [ 0 1 ]
if u < alpha % Do we accept this proposal?
theta(:,t) = theta_star; % proposal becomes new theta
else
theta(:,t) = lasttheta; % copy old theta
end
if mod( t,10 ) == 0
fprintf( 't=%3d	' , t );
for j=1:L
fprintf( '%15s' , labels{ theta(j,t)} );
end
fprintf( '
' );
end
end