I am trying to find the coefficients of a system of linear regression equations defined as y = b*X + ε, where ε ~ N(0,s^2). The coefficients are to be constrained so that each b_i is between 0 and 1 and sum(b_i) = 1. To do this we are trying to define w = exp(b) / sum(exp(b)) and then solve for w. Inside my linregeval() function I make this transformation, obtain the X (data) matrix from userinfo, then set the likelihood equal to lnnormalden(y, X*w’, sqrt(p2)), where p2 is the variance parameter.

When I run moptimize() it fails in “Iteration 1” for “flat or discontinuous region encountered” and the likelihood value does not change upon each sub-iteration (each time linregeval() is called within “Iteration 1”). The value of b also doesn’t change, but if I print p1 (obtained from util_xb) it appears to update on each call to linregeval(). However, I can’t obtain b from p1 since this would involve “dividing out” the X matrix.

Is there any way to solve this problem?