st: Accept-reject procedure and -drawnorm- for simulating likelihood in -ml-

[Terence Cheng <terence.cheng@anu.edu.au>(by way of Keith Dear <keith.dear@anu.edu.au>)]

Dear Statalisters

I am writing a program that calculates the log likelihood value for 
use in -ml-. The density function that I am maximising does not have 
a closed form solution as it depends on three random coefficients. 
Simulating the density function requires drawing from the truncated 
trivariate normal distribution which I am attempting to do using the 
accept-reject method.  This first involves generating a set of three 
random numbers using -drawnorm- and dropping the observations that 
are not within the accepted range, using the accepted observations to 
simulate the density.

The problem with the accept-reject approach as described above is 
that I would need to drop the data (using -drop _all-), which also 
drops the temporary variables that Stata generates for -ml- such as 
the lnf and coefficient vectors. I attempted to use scalars to 
capture the values of the coefficient vectors as I require them for 
the intermediate calculations. I have also tried to preserve the data 
prior to dropping them and for which Stata automatically restores 
after the program ends.

There must have been something wrong with my approach because while 
implementing -ml check-, Stata returns the following message:

"The initial values are not feasible.  This may be because the 
initial values have been chosen poorly or because there is an error 
in progname and it always returns missing no matter what the 
parameter values."

I would appreciate if someone could point out where I had gone wrong. 
Is there a better way to generate draws from a truncated multivariate 
normal distribution? Thanks.

The beginning segment of my code is as follows:

****************
cap prog drop progname
program define poinorm_mlv3
version 9.0
args lnf theta1 theta2 theta3 theta4 theta5 theta6
/*Capture coefficient vectors as scalars prior to drop _all*/
tempname rho12 rho13 rho23 /*Correlation matrix*/
scalar `rho12'=[exp(2*`theta4')-1]/[exp(2*`theta4')+1]
scalar `rho13'=[exp(2*`theta5')-1]/[exp(2*`theta5')+1]
scalar `rho23'=[exp(2*`theta6')-1]/[exp(2*`theta6')+1]
if (`rho12'==. | `rho13'==. | `rho23'==.) {
qui replace `lnf'=.
exit
}

matrix Sig = (1 , `rho12' , `rho13' \ /* /*Construct corr mat*/
*/`rho12' , 1 , `rho23' \ /*
*/`rho13' , `rho23' , 1 )
matrix symeigen X V = Sig /*Test for positive df*/
if V[1,3] <= 0 {
qui replace `lnf' =.
exit
}

mat M = inv(Sig)
local repl = 5 /*Number of replications*/
local n = _N /*Sample size for density function*/
tempname xb1 xb2 xb3 y1 y2 y3
scalar `xb1' = `theta1' /*Capture xbs*/
scalar `xb2' = `theta2'
scalar `xb3' = `theta3'
scalar `y1' = $ML_y1 /*Capture depvars*/
scalar `y2' = $ML_y2
scalar `y3' = $ML_y3
preserve
qui drop _all /*drop data and variables*/
/* Using -drawnorm- to generate random numbers, accept and reject. */


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