Statalist


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

RE: st: help with my simple integration command in mata


From   Glenn Goldsmith <glenn.goldsmith@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   RE: st: help with my simple integration command in mata
Date   Tue, 16 Jun 2009 00:36:52 +0100

Hi Nathan,

1. The error is occurring because your test function -funky()- needs
to -return()- 10*x rather than just displaying it as it does now.
Replacing the -10*x- line with -return(10*x)- should do the trick.

2. It looks like there might be something else slightly odd going on
with your program though. Your loop calculates 15 values (0 through
14) but your return statement only uses 14 of these (2 through 15). I
assume it should use all 15, with the multipliers
(1,4,2,4,2,...,2,4,1) rather than the (1,2,4,2,...,2,4,1) you
currently have?

3. In terms of general debugging, mata can be a little difficult; but
inserting -mata set matalnum on- after your -mata clear- line will at
least give you line numbers for the errors. (You'll want to turn it
off after testing though, because when it's on it stops the compiler
from optimizing the code.)

4. There are a couple of other ways you could clean up your code to
achieve what you want slightly more efficiently. A modified version is
below:

**********************************************************************************
mata
mata clear
mata set matalnum off

scalar NDintegration(								// to perform integration via
simpson'sapproximation
pointer(real scalar function) scalar f,						//points to the
functionto be integrated
real scalar lower,									//lower bound of integration
real scalar higher)									//upper bound
{

  real scalar width
                  //I like declarations, your mileage may vary
  real colvector Y
  real rowvector X

  if (lower > higher) _error("limits of integration reversed")		//make
sure upper > lower NB: better to use _error() than return()

  width = (higher-lower) / 14								//calculates width of each "slice"

  Y = J(15,1,.)										//creates an empty 15 x 1 column vector
  for (i=1; i<=15; i++) {
    x = lower + ((i-1)*width)
    Y[i,1]  = (*f)(x)								        //places function evaluations
directly into Y, avoiding all the copying you were doing previously
  }

  X = (1,4,2,4,2,4,2,4,2,4,2,4,2,4,1)
            //amended 1 x 15 row vector of multipliers
  return(X*Y*width/3)
                 //matrix multiplication makes this quicker and more
readable

}

scalar funky(real scalar x) {
              //declares the output to be a scalar. otherwise you
could run into problems in the main program
  return(10*x)
}

f = &funky()

NDintegration(f, 0, 10)

end
***********************************************************

HTH,

Glenn.


Nathan Danneman <ndanneman@gmail.com> wrote:

Hi all,

I have written the following bit of code that uses Simpson's method to
approximate definite integrals.  It doesn't seem to work: it responds
that there is a conformability error, which I cannot seem to locate.

I would welcome any general tips on debugging programs in mata, and
especially welcome any specific tips on refining this code.

Thanks so much,
Nathan


mata
mata clear

function NDintegration(								// to perform integration via simpson's
approximation
pointer(real scalar function) scalar f,						//points to the function
to be integrated
real scalar lower,									//lower bound of integration
real scalar higher)									//upper bound
{

if (lower > higher) return("limits of integration reversed")		
//make sure upper > lower

width = (higher-lower) / 14								//calculates width of each "slice"

Y = (.)											//creates an empty matrix

	for (i=0; i<=14; i++) {
		x = lower + (i*width)
		yi = (*f)(x)								//is this part correct?  I'm not sure I've used
the pointer correctly
		d = yi
		Y = (Y, d)
	}
st_matrix("Y", Y)


return( (width/3) * (1*st_matrix("Y")[1,2] + 2*st_matrix("Y")[1,3] +
4*st_matrix("Y")[1,4] + 2*st_matrix("Y")[1,5] + 4*st_matrix("Y")[1,6]
+ 2*st_matrix("Y")[1,7] + 4*st_matrix("Y")[1,8] +
2*st_matrix("Y")[1,9] + 4*st_matrix("Y")[1,10] +
2*st_matrix("Y")[1,11] + 4*st_matrix("Y")[1,12] +
2*st_matrix("Y")[1,13] + 4*st_matrix("Y")[1,14] +
1*st_matrix("Y")[1,15]) )

//this last bit simply sums up the y values and multiplies them by the
appropriate fudge factors

}



function funky(real scalar x) {                       //I tried to use
this simple function to test the integration routine...no luck
	10*x
}

f = &funky()

NDintegration(f, 0, 10)


end
*
*   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