/*
**   circular.e  -  A circular regression model
**
**   The log-likelihood for circular regression contains local
**   minima where the coefficients become very large.  When
**   that happens, restart from another point.
*/


/* Simulate circular data */

  rndseed 34534567;
  nobs = 500;
  numInd = 2;

  mu = 0;
  k = 1;

  b = .5 * ones(numInd,1);
  x = rndu(nobs,numInd);
  y = rndvm(nobs,1,mu + g(x*b),k);



  b0 = 5*ones(NumInd+2,1);

 output file = circular.out reset;


 { b,fct,grd,ret } = QNewton(&lpr,b0);

  cov = invpd(hessp(&lpr,b));

  print "coefficients   standard errors";
  print;
  print b~sqrt(diag(cov));

 output off;


proc G(u);
   retp(2*atan(u));
endp;

proc lpr(b);
    local dev;
    dev = y - b[2]- G(b[3] + x * b[3:rows(b)]);
    retp(sumc(ln(mbesselei0(exp(b[1]))) - exp(b[1]) * (cos(dev) - 1)));
endp;