Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: creating a new user function in Stata


From   "Carlos Madeira" <carlos_m_madeira@hotmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: creating a new user function in Stata
Date   Sat, 10 Feb 2007 16:24:47 +0000

Hi Maarten,

thanks for your answer. I think it was very clear.

Actually, I didn't think before of the problem you pointed out with the small numbers, but I think I can find a way around it. I think the problem is that in general we define a beta function in the space [0-1], so if I think about it that will necessarily imply that the numbers involved for the standard-deviations and the correlation are very small. However, I think there is a simple way around that. If I translate my natural "beta-space" to be equivalent to [0-100], then perhaps things will work well.

Best,

Carlos






From: Maarten buis <maartenbuis@yahoo.co.uk>
Reply-To: statalist@hsphsun2.harvard.edu
To: statalist@hsphsun2.harvard.edu
Subject: Re: st: creating a new user function in Stata
Date: Sat, 10 Feb 2007 11:47:28 +0000 (GMT)

--- Carlos Madeira <carlos_m_madeira@hotmail.com> wrote:
> I had a problem writing a user-written function in Stata and was
> wondering if you could give me a help about how to write it. My goal
> is to create a program that generates a variable with cdf values of a
> bivariate-beta function in a similar way that the function i
> beta(a1,b1,x1) would do.

Short answer is that functions are hard coded into Stata, so you can't
write your own function. You can write an -egen- program. But given
your example I suspect you want this to be part of a ml program. In
that case I wouldn't bother and just break the equation up in manageble
parts and create temporary variables. So looking at the equation you
gave us:
> cdf_biv_beta =
> ibeta(a1,b1,x1)(1+w(x1-a1/(a1+b1))(x2-a2/(a2+b2))ibeta(a2,b2,x2) +
> (....)

I would do:
tempvar mean1 mean2 temp
gen double `mean1' = `a1'/(`a1' + `b1')
gen double `mean2' = `a2'/(`a2' + `b2')
gen double `temp' = 1 + `w'*(`x1'-`mean1')*(`x2'-`mean2')
replace `lnf' = ln(ibeta(`a1',`b1',`x1') * /*
*/ `temp' * /*
*/ ibeta(`a2',`b2',`x2') + /*
*/ ... )

I have two additional comments:
1. When writing down the log likelihood it is worth the effort to see
if taking the log simplifies the equation (turns multiplications in
additions, and the like). An example is shown in a short technical text
on http://home.fsw.vu.nl/m.buis/software/betafit.html .
2. the part 1 + `w'*(`x1'-`mean1')*(`x2'-`mean2') worries me from a
numerical point of view. `x1'-`mean1' is the difference between an
observed and expected proportion, so this is usually a small number.
You multiply it with a similarly small number: `x2'-`mean2'. Adding
such a small number to 1 can lead to big rounding errors. The last Mata
matters column in the Stata Journal (volume 6, number 4) by William
Gould was very clear piece on this issue.

Hope this helps,
Maarten


-----------------------------------------
Maarten L. Buis
Department of Social Research Methodology
Vrije Universiteit Amsterdam
Boelelaan 1081
1081 HV Amsterdam
The Netherlands

visiting address:
Buitenveldertselaan 3 (Metropolitan), room Z434

+31 20 5986715

http://home.fsw.vu.nl/m.buis/
-----------------------------------------



___________________________________________________________
All New Yahoo! Mail Tired of unwanted email come-ons? Let our SpamGuard protect you. http://uk.docs.yahoo.com/nowyoucan.html
*
* 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/
_________________________________________________________________
Don't just search. Find. Check out the new MSN Search! http://search.msn.com/

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