程序见后文，运行结果是
       1.0000000        1.0000000        1.7320508
       1.0473677       0.87799707      -0.22331818
前两个值是系数beta1和beta2，最后一个是sigma。前两个值还可以，最后一个一定是哪里有问题，请高手帮忙解答，先表示感谢！

new;
library pgraph;
graphset;
pqgwin many;
rndseed 19890930;   //set the seed
/***********data generation***********/
b1 = 1;     b2 = 1;     beta = b1 | b2;     
sig_sq = 3;     sig = sqrt(sig_sq);
t_par = b1| b2 | sig;
n = 300;
x1 = ones(n,1);
x2 = rndu(n,1);
x = x1 ~ x2;
e = sig*rndn(n,1);
y = x*beta+e;
/**********initial values***********/
beta = y/x;
sig = stdc(y - x*beta);
p_f1 = pos(y, x, beta, sig);  //get initial posterior density
/***********MCMC loop***********/
burn = 100;     //burn the first 100 samples
iter = 3000;
nrep = burn + iter;
mcmc = zeros(nrep,3);
/*********function definition*********/
proc pos(y,x,beta,sig);
    local n,u,f;
    n = rows(y);
    u = (y-x*beta)' * (y-x*beta);
    f=((2*pi)^(-n/2))*(sig^(-n))*exp(-u/(2*(sig^2)));
retp(ln(f));
endp;   
fn drawpar(mean, sigma) = mean+chol(sigma)'*rndn(rows(mean),1);


ixx = inv(x'*x);
sig2 = sig^2;
beta_sig = sig2*ixx;     //var(beta)=sigma^2*(x'*x)^(-1)

for i(1,nrep,1);
    betatem = drawpar(beta,beta_sig);
    p_f = pos(y,x,betatem,sig);
    qb01 = 0;
    temp = 1 | exp(p_f - p_f1 + qb01);
    if (minc(temp) > rndu(1,1));
        beta = betatem;
        p_f1 = p_f;
    endif;
    mcmc[i,1:2] = beta';
endfor;

for i(1,nrep,1);
    sig_sig = 0.3;//setting by experience
    sig_tem = drawpar(sig,sig_sig);
    p_f = pos(y,x,beta,sig_tem);
    qb01 = 0;
    temp = 1 | exp(p_f - p_f1 + qb01);
    if (minc(temp) > rndu(1,1));
        sig = sig_tem;
        p_f1 = p_f;
    endif;
    mcmc[i,3] = sig;
endfor;

index = seqa(0,1,iter);
title("MCMC draws for beta1");
xy(index,mcmc[burn+1:nrep,1]);
title("MCMC draws for beta2");
xy(index,mcmc[burn+1:nrep,2]);
title("MCMC draws for sig");
xy(index,mcmc[burn+1:nrep,3]);
m_mcmc = meanc(mcmc);
print t_par';
print m_mcmc';
end;