[time,y] = textread('pumps.txt','%f %f')
n=10;T=25000; lam=ones(n,1);beta=0.9*ones(1,T); scale=0.2;a.alph=0.1; nu=-0.4*ones(1,T);
a.beta=0.1; b.beta=1; alph(1)=exp(nu(1));
for t=1:T     for i=1:n
loglam(i,t)=log(lam(i,t));end
P=exp(nu(t)-a.alph*alph(t)+n*alph(t)*log(beta(t))...
    -n*gammaln(alph(t))+(alph(t)-1)*sum(loglam(1:n,t)));
nustar=nu(t)+ scale*randn;
alphstar=exp(nustar);
Pstar=exp(nustar-a.alph*alphstar+n*alphstar*log(beta(t))...
    -n*gammaln(alphstar)+(alphstar-1)*sum(loglam(1:n,t)));
if    (rand <= Pstar/P)     alph(t+1)=exp(nustar);
else                        alph(t+1)=alph(t);  end
% update parameters from full conditionals
for i=1:n
lam(i,t+1)=gamrnd(alph(t+1)+y(i),1/(beta(t)+time(i)));end
beta(t+1)=gamrnd(a.beta+n*alph(t+1),1/(b.beta+sum(lam(1:n,t+1))));
% accumulate draws for coda input
for i=1:n pars(t,i)=lam(i,t); end 
pars(t,n+1)=beta(t); pars(t,n+2)=alph(t);  end
  hist(beta,100); pause;   hist(alph,100); pause;
  [hbeta,smbeta,xbeta] = pltdens(beta);   plot(xbeta,smbeta); pause;
  [halph,smalph,xalph] = pltdens(alph);   plot(xalph,smalph); pause;
  coda(pars)