cap pr drop myreg4

program define myreg4

	version 8 
	args todo b lnf g negH

	tempvar theta
	tempname sigma2 
	mleval `theta' 		= `b', eq(1) 
	mleval `sigma2' 	= `b', eq(2) scalar

	mlsum `lnf' = -($ML_y1-`theta')^2/(2*`sigma2')-0.5*ln(`sigma2') 
	if (`todo' ==0 | `lnf' ==.) exit 

	tempvar fb fs

	mlvecsum `lnf' `fb' = ($ML_y1-`theta')/`sigma2',	 	eq(1) 
	mlvecsum `lnf' `fs' = ($ML_y1-`theta')^2/(2*`sigma2'^2)-1/(2*`sigma2'),	eq(2)

	matrix `g' = (`fb',`fs') 
	if (`todo'==1 | `lnf'==.) exit 

	tempvar fbb fbs fss

	mlmatsum `lnf' `fbb'	= -1/`sigma2',              						eq(1) 
	mlmatsum `lnf' `fbs'	= -($ML_y1-`theta')/(`sigma2'^2),      			    eq(1,2) 
	mlmatsum `lnf' `fss'	= -($ML_y1-`theta')^2/(`sigma2'^3)+1/(2*`sigma2'^2),	eq(2) 

	matrix `negH' = (`fbb',`fbs' \ `fbs'',`fss')

end

ml model d2 myreg4 (reg:y=x1 x2) /sigma2

ml maximize, //grad hess