You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
130 lines
2.9 KiB
Plaintext
130 lines
2.9 KiB
Plaintext
*! version 1.0.0 7 June 1995 sg41: STB-26
|
|
program define rfpr_ll1 /* b log_likelihood grad */
|
|
version 4.0
|
|
local b "`1'"
|
|
local f "`2'"
|
|
local g "`3'"
|
|
macro shift 3
|
|
local options "FIRSTIT LASTIT FAST(string)"
|
|
parse "`*'"
|
|
|
|
local y = trim("$S_mldepn")
|
|
local doit "$S_sample"
|
|
local i "$S_ivar"
|
|
|
|
tempname rho s2su beta u
|
|
tempvar xb F p
|
|
|
|
local rhocol = colnumb(`b',"rho:_cons")
|
|
scalar `rho' = abs(`b'[1,`rhocol'])
|
|
|
|
if `rho' >= 1 {
|
|
scalar `rho' = 0.99
|
|
di "rho >= 1, set to rho = 0.99"
|
|
}
|
|
|
|
matrix `b'[1,`rhocol'] = `rho'
|
|
scalar `s2su' = sqrt(2*`rho'/(1 - `rho'))
|
|
matrix `beta' = `b'[1,"`y':"]
|
|
|
|
quietly {
|
|
|
|
matrix score double `xb' = `beta' if `doit'
|
|
|
|
gen double `F' = . in 1
|
|
by `doit' `i': gen double `p' = cond(_n==_N,0,.) if `doit'
|
|
|
|
/* Do computation this way if only log likelihood required. */
|
|
|
|
if "`fast'" == "0" {
|
|
local m 1
|
|
while `m' <= $S_quad {
|
|
scalar `u' = `s2su'*$S_x[`m']
|
|
|
|
#delimit ;
|
|
|
|
by `doit' `i': replace `F' =
|
|
cond(_n==1,
|
|
cond(`y', normprob(`xb' + `u'),
|
|
1 - normprob(`xb' + `u')),
|
|
cond(`y', normprob(`xb' + `u'),
|
|
1 - normprob(`xb' + `u'))
|
|
*`F'[_n-1])
|
|
if `doit' ;
|
|
|
|
#delimit cr
|
|
|
|
replace `p' = `p' + $S_w[`m']*`F' if `doit'
|
|
|
|
local m = `m' + 1
|
|
}
|
|
|
|
replace `F' = sum(log(`p'/sqrt(_pi))) if `doit'
|
|
scalar `f' = `F'[_N]
|
|
|
|
exit
|
|
}
|
|
|
|
/* Do computation this way if first derivatives required. */
|
|
|
|
tempname gr
|
|
tempvar db dr
|
|
gen double `db' = 0
|
|
gen double `dr' = 0
|
|
|
|
local m 1
|
|
while `m' <= $S_quad {
|
|
scalar `u' = `s2su'*$S_x[`m']
|
|
|
|
#delimit ;
|
|
|
|
by `doit' `i': replace `F' =
|
|
cond(_n==1,
|
|
cond(`y', normprob(`xb' + `u'),
|
|
1 - normprob(`xb' + `u')),
|
|
cond(`y', normprob(`xb' + `u'),
|
|
1 - normprob(`xb' + `u'))*`F'[_n-1])
|
|
if `doit' ;
|
|
|
|
replace `p' = `p' + $S_w[`m']*`F' if `doit' ;
|
|
|
|
by `doit' `i': replace `F' =
|
|
cond(`y',1,-1)*exp(-0.5*(`xb'+`u')^2)*`F'[_N]
|
|
/cond(`y',normprob(`xb'+`u'),1-normprob(`xb'+`u'))
|
|
if `doit' ;
|
|
|
|
#delimit cr
|
|
|
|
replace `db' = `db' + $S_w[`m']*`F' if `doit'
|
|
replace `dr' = `dr' + $S_w[`m']*`u'*`F' if `doit'
|
|
|
|
local m = `m' + 1
|
|
}
|
|
|
|
/* Compute log likelihood. */
|
|
|
|
replace `F' = sum(log(`p'/sqrt(_pi))) if `doit'
|
|
scalar `f' = `F'[_N]
|
|
|
|
/* Compute first derivatives. */
|
|
|
|
by `doit' `i': replace `p' = `p'[_N] if `doit'
|
|
|
|
replace `db' = `db'/(sqrt(2*_pi)*`p') if `doit'
|
|
|
|
local nvar = colsof(`beta') - 1
|
|
if `nvar' > 0 {
|
|
matrix `beta' = `beta'[1, 1..`nvar']
|
|
local vars : colnames(`beta')
|
|
}
|
|
|
|
matrix vecaccum `g' = `db' `vars' if `doit'
|
|
|
|
replace `dr' = sum(`dr'/`p') if `doit'
|
|
|
|
matrix `gr' = (0)
|
|
matrix `gr'[1,1] = `dr'[_N]/(2*sqrt(2*_pi)*`rho'*(1-`rho'))
|
|
matrix `g' = `g' , `gr'
|
|
}
|
|
end
|