Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: Random walk simplex


From   Austin Nichols <austinnichols@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Random walk simplex
Date   Fri, 23 Oct 2009 11:02:18 -0400

Akshay Shanker <Akshay.Shanker@standards.org.au> :
Try out:

clear
set obs 1000
g long i=_n
g byte t=1
loc k=3
forv i=1/`k' {
g double x`i'=-ln(uniform()) if t==1
}
loc maxt=6
qui forv plust=2/`maxt' {
set obs `=_N+1000'
replace i=mod(_n-1,1000)+1 if mi(i)
bys i (t): replace t=t[_n-1]+1 if mi(t)
tsset i t
forv i=1/`k' {
loc N=_N
loc n=0
while `n'<`N' {
g double x=l.x`i'*exp(rnormal(0,1)) if mi(x`i')
replace x`i'=x if mi(x`i')&exp(-x)/exp(-l.x`i')*(x/l.x`i')>uniform()
drop x
count if x`i'<.
loc n=r(N)
}
}
}
egen double s=rowtotal(x*)
forv i=1/`k' {
g double t`i'=x`i'/s
}
sc t1 t2, by(t)
ren t3 origt3
g double t3=1-t1-t2
egen double sum=rowtotal(t?)
su


On Thu, Oct 22, 2009 at 10:55 PM, Akshay Shanker
<Akshay.Shanker@standards.org.au> wrote:
> Hi all,
>
> I was hoping someone would be able to  help me out.
>
> I am trying to create a random walking simplex in STATA. The procedure
> to create a random K dimension simplex seems straightforward:
>
> - generate K uniform RVs Xi...Xk in (0,1)
> - each point in the random simplex is then i/(sum(-logXi)).
>
> But how would I go about getting this simplex to do a random walk in
> Stata? More specifically, how can I choose a markov like draw of
> -logXi,t based on an already given -logXi,t-1. There seems to be a
> metropolis-hastings algo which selects this draw from a gamma
> distribution with a mean of Xi,t-1 - but not sure how to do this in
> STATA. Ideally I'd like to be able to generate a KxN matrix where the
> coloumns are the simplex, which walk over the rows.
>
> Regards,
> Akshay

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



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