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

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/

- Prev by Date:
**Re: st: Re: Generating index numbers with panel data in long format** - Next by Date:
**Re: st: xtreg using independent variables from previous wave** - Previous by thread:
**st: help with my simple integration command in mata** - Next by thread:
**Re: st: Re: upper function within the loop** - Index(es):

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