# RE: st: RE: how to constrain rowtotal?

 From "Nick Cox" <[email protected]> To <[email protected]> Subject RE: st: RE: how to constrain rowtotal? Date Mon, 28 Jan 2008 14:18:20 -0000

```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/
```