Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

RE: st: precision problems with Mata? naaaahh...


From   "Jann Ben" <ben.jann@soz.gess.ethz.ch>
To   <statalist@hsphsun2.harvard.edu>
Subject   RE: st: precision problems with Mata? naaaahh...
Date   Wed, 17 Aug 2005 19:55:44 +0200

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: owner-statalist@hsphsun2.harvard.edu 
> [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of 
> Nicola Orsini
> Sent: Wednesday, August 17, 2005 5:54 PM
> To: statalist@hsphsun2.harvard.edu
> 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/



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