capture log close
log using lem_mata.log, replace

// 2 by 2 table, missing values on one margin.
// same example as solved by LEM
// the SV matrix (col vector) includes the 
// 4 (s,v) joint probabilities and 2 marginal (s,.)
// p is the col vector of 4 joint probabilities Pr(S=s, V=v)
// C is the row vector to multiply p; c is the element 1
// We impose the constraint C p = c, that is, the 4 probabilities add up to 1
// Cc is the row vector that is passed to Mata

mata:  // start Mata

	mata clear
	SV = (28, 45, 22, 52 , 10, 15)' 
	C = (1, 1, 1, 1)
	c = (1)
	Cc = (C,c)
	void myfun(todo, p,  SV,  lnf, g, H)
	{
	 lnf = SV[1]*log(p[1]) + SV[2]*log(p[2]) ///
		+ SV[3]*log(p[3]) + SV[4]*log(p[4]) 

	  lnf = lnf + SV[5]*log(p[1]+p[2]) ///
		+ SV[6]*log(p[3]+p[4])
	}

	S = optimize_init()

	optimize_init_evaluator(S, &myfun()) // optimize liklihood function
	optimize_init_params(S, (.25, .25, .25, .25)) // initial values of p
	optimize_init_constraints(S, Cc)  // constraints
	optimize_init_argument(S, 1, SV)                            
	optimize(S)
	optimize_result_V_oim(S)

	p = optimize(S) // estimated probabilities
	varcov_p = optimize_result_V_oim(S) // var-cov matrix of estimates
	var_p = diagonal(varcov_p) // variances of estimates
	se_p = sqrt(var_p) // se's of estimates
		
	round(p, .0001)  // print out estimated probabilities w/ 4 decimals
	round(se_p', .0001)  // print out row of estimated se's w/ 4 decimals

end // exit Mata

log close

exit