Kit and Nicola,
many thanks for your suggestions. However, note that creating
the lowertriangle matrix is slow and uses lots of memory (try
to do this with, say, 100'000 cases).
Here are some timings with 2000 cases (code below):
1) Approach by Kit/Nicola: t=0.31
2) Same using -cross()-: t=0.35
3) My initial approach: t=0.01
4) My approach using -quadsum()-: t=0.14
Additional cases disproportionately increase computing time for
1, 2, and 4 (and linearly increase time for 3). Memory usage
increases disproportionately for 1 and 2, but not for 3 and 4.
So, the only practicable method for my purposes seems to be
approach 3. I guess I have to live with the precision problem
I mentionnend. It only occurs rarely.
ben
Code:
clear
mata:
real colvector runsum1(real colvector x)
{
return(lowertriangle(J(rows(x),rows(x),1))*x)
}
real colvector runsum2(real colvector x)
{
return(cross(uppertriangle(J(rows(x),rows(x),1)),x))
}
real colvector runsum3(real colvector x)
{
r = x
for (i=2; i<=rows(r); i++) r[i] = r[i-1] + r[i]
return(r)
}
real colvector runsum4(real colvector x)
{
r = J(rows(x),1,.)
for (i=1; i<=rows(r); i++) r[i] = quadsum(x[|1 \ i|])
return(r)
}
end
mata: x = uniform(2000,1)
mata: sum1 = runsum1(x)
mata: sum2 = runsum2(x)
mata: sum3 = runsum3(x)
mata: sum4 = runsum4(x)
> -----Original Message-----
> From: [email protected]
> [mailto:[email protected]] On Behalf Of
> Nicola Orsini
> Sent: Wednesday, August 17, 2005 5:54 PM
> To: [email protected]
> Subject: Re: st: precision problems with Mata? naaaahh...
>
>
> Hi Kit,
>
> to get the lower triangle -help mf_lowertriangle-
>
> mata:
> void function runsum(string scalar vname)
> {
> real matrix L
> real matrix runsum
> real scalar resindex
> real scalar vnew
> vnew = vname+"_sum"
> st_view(A,.,vname)
>
> L=lowertriangle(J(rows(A),rows(A),1))
>
> runsum = L*A
> resindex = st_addvar("double",vnew)
> st_store((1,rows(runsum)),resindex,runsum)
> }
> end
>
> Nicola
>
>
> Kit Baum wrote:
> > Here's a Mata solution (no doubt hideously inefficient, but
> I can't
> > find all the fns. I'm looking for) that appears to be just
> as precise
> > as a -generate, sum- at doing a running sum. It could
> readily operate
> > on a column vector instead of a Stata variable if that is desired.
> >
> > Note that a running sum can be calculated by premultiplying
> by a square
> > matrix with 1's in the lower triangle, 0's above...
> >
> > clear
> > set obs 1000
> > g double x = (uniform()-0.5)*10^(int(_n/10))
> > g sort = uniform()
> > sort sort
> > drop sort
> > g double sum = sum(x)
> >
> > mata:
> > void function runsum(string scalar vname)
> > {
> > real matrix L
> > real matrix runsum
> > real scalar resindex
> > real scalar vnew
> > vnew = vname+"_sum"
> > st_view(A,.,vname)
> > // if(cols(A)>1) A=A'
> > L = J(rows(A),rows(A),0)
> > for (i=1; i<=rows(A); i++) {
> > for(j=1; j<=i; j++) {
> > L[i,j] = 1
> > }
> > }
> > runsum = L*A
> > resindex = st_addvar("double",vnew)
> > st_store((1,rows(runsum)),resindex,runsum)
> > }
> > end
> >
> > mata runsum("x")
> > g double discrep = sum - x_sum
> > su
> > l, sep(0)
*
* 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/
