I ran Program 1 in WinBUGS 1.3 and it worked as expected:
- xsum=3 every time in one million iterations.
- X and pr results looked okay
Finn Krogstad
U. Washington, Forest Engineering
http://students.washington.edu/fkrogsta
---------- Forwarded message ----------
Date: Wed, 14 Jun 2000 16:09:12 -0400 (EDT)
From: Belisle Patrick <[log in to unmask]>
To: [log in to unmask]
Subject: Bug with dbern??
Dear all,
this is related to Jukka Ranta's last message to the Bugs list (yesterday).
Suppose you want to sample n subjects out of N (SI sample):
Fan, Muller and Rezucha (1962) have described the following algorithm
to do it:
1/ define a[k] the number of subjects included amongst k-1 first subjects
2/ define pr[k] <- ( n - a[k] ) / ( N + 1 - k )
3/ Generate U[k] ~ uniform(0,1)
4/ Include subject k if U[k] <= pr[k]
[ Define X[k] <- step(pr[k] - U[k]); ]
It can be shown that prob{X[j] = 1} = n/N for j=1,2,..,N
and sum(X[]) = n
I tried to implement it in Bugs, using X[i] ~ dbern(pr[i]),
but even if at some point pr[]'s become 0.0 --- suppose that
you want to pick 3 subjects out of 8 and that the first 3 are selected,
then pr[4] = pr[5] = .. = pr[8] = 0.0 --- , Bugs will still
return X[i] = 1 for some of the i's concerned.
Unless I am missing a point here, this seems like a bug, doesn't it?
(I am using classic Bugs 0.6, see program 1 below)
However, if I apply more directly the algorithm below (using dunif()
rather than dbern(), see pgm 2 below), it works!!! I am puzzled...
Thanks in advance for any help / comments.
Patrick Bélisle
Montréal, Qc, CAN
Dept of Clin Epidemiology, Mtl Gen Hospital
--- Program 1 ----------------------
model si;
const n=3, N=8;
var X[N], pr[N], a[N], xsum;
{
a[1] <- 0;
pr[1] <- n/N;
X[1] ~ dbern(pr[1]);
for (i in 2:N)
{
a[i] <- a[i-1] + X[i-1];
pr[i] <- (n-a[i])/(N+1-i);
X[i] ~ dbern(pr[i]);
}
xsum <- sum(X[]);
# If you monitor xsum, X and pr, update the program a few iterations,
# and sooner or later you will get a value of xsum different from 3:
# then list the values of X and pr, and you will some inconsistencies...
}
--- Program 2 ----------------------
model si;
const n=3, N=8;
var X[N], pr[N], a[N], xsum, U[N];
{
a[1] <- 0;
pr[1] <- n/N;
U[1] ~ dunif(0,1);
X[1] <- step(pr[1]-U[1]);
for (i in 2:N)
{
a[i] <- a[i-1] + X[i-1];
pr[i] <- (n-a[i])/(N+1-i);
U[i] ~ dunif(0,1);
X[i] <- step(pr[i]-U[i]);
}
xsum <- sum(X[]);
# This program gives the expected results.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|