Re: st: Iteratively re-weighted least squares (IRLS)

[Joseph Coveney <jcoveney@bigplanet.com>]

Just to add to the earlier reply, the -lnlsq(0)- option actually fits what
pharmacokineticists call the "exponential" variance model, which is (or at
least it used to be) often used by them to approximate a "proportional"
variance model such as yours.

Perhaps a better approach would be to explicitly and directly fit a
"proportional" error model, such with as the Nelder-Mead simplex algorithm.
The simplex algorithm is implemented in Stata as a user-written
command, -amoeba-.  I've illustrated its use in your problem in the do-file
below, and would recommend it over -nl , lnlsq()-, although the 
approximation with the latter doesn't look too shabby in the illustrative 
case and wouldn't need resampling techniques for confidence limits or 
hypothesis testing.

Joseph Coveney

clear
set more off
set seed `=date("2007-04-07", "ymd")'
set obs 500
// C = (a + b*X^g)*e
local a 2
local b 3
local g 0.75
generate double X = uniform()
generate double e = 0.15 * invnormal(uniform())
generate double C = (`a' + `b' * X^`g') * (1 + e)
*
capture program drop proerrmod
program define proerrmod
   version 9.2
   tempvar intermediate residual
   quietly {
       generate double `intermediate' = ///
         `1'[1,1] + `1'[1,2] * $xvar^(`1'[1,3])
       generate double `residual' = ($yvar / `intermediate') - 1
       replace `residual' = `residual' * `residual'
   }
   summarize `residual', meanonly
   scalar define `2' = -r(sum)
end
*
matrix define init = (1, 1, 1)
scalar define obj = .
matrix define B = (., ., .)
matrix colnames init = a:_cons b:_cons g:_cons
global yvar C
global xvar X
amoeba proerrmod init obj B
exit

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