Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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/



© Copyright 1996–2024 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index