Good catch!
Nick
[email protected]
-----Original Message-----
From: [email protected]
[mailto:[email protected]] On Behalf Of Maarten buis
Sent: 28 January 2008 11:18
To: [email protected]
Subject: RE: st: RE: how to constrain rowtotal?
--- Nick Cox <[email protected]> wrote:
> This may help:
>
> . gen which = cond(uniform() < 0.721, "A",
> cond(uniform() < 0.923, "B",
> cond(uniform() < 0.990, "C", "D")))
-- Maarten buis wrote:
> Even better would be:
> gen unif = uniform()
> gen which = cond(unif < 0.721, "A", ///
> cond(unif < 0.923, "B", ///
> cond(unif < 0.990, "C", "D")))
--- Nick Cox <[email protected]> wrote:
> I am not clear why you think that better.
The idea is to create a random variable with 4 categories whereby the
probability of falling in each category is know, in this case A= 0.721;
B= 0.202; C= 0.067; D= 0.010. This can be achieved by mapping a single
draw from a uniform distribution onto the four values. The logic is
that the probability that such a draw from a uniform distribution is
less than .721 is .721, between .721 and .923 is .202, etc.
In your code you used a new draw from a uniform distribution to assign
people to categories. So there are two reasons why a new draw from
uniform() could be less than .923: it is between .721 and .923 or it is
less then .721. In my code I exclude the latter possibility. The
consequences are shown in the simulation below. For simplicity I use
only three categories with probabilities 1=0.3; 2=0.3; 3=0.4.
-- Maarten
*------------------ end example ------------------------------
capture program drop sim
program define sim, rclass
drop _all
set obs 1000
gen unif = uniform()
gen byte buis = cond(unif < .3, 1, ///
cond(unif < .6, 2, 3))
count if buis == 1
return scalar buis1 = r(N)/1000
count if buis == 2
return scalar buis2 = r(N)/1000
count if buis == 3
return scalar buis3 = r(N)/1000
gen byte cox = cond(uniform() < .3, 1, ///
cond(uniform() < .6, 2, 3))
count if cox == 1
return scalar cox1 = r(N)/1000
count if cox == 2
return scalar cox2 = r(N)/1000
count if cox == 3
return scalar cox3 = r(N)/1000
end
simulate buis1=r(buis1) buis2=r(buis2) buis3=r(buis3) ///
cox1=r(cox1) cox2=r(cox2) cox3=r(cox3), ///
reps(1000): sim
sum
*------------------------- begin example ------------------------
(For more on how to use examples I sent to the Statalist, see
http://home.fsw.vu.nl/m.buis/stata/exampleFAQ.html )
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
