proc nlmixed data=yogurt;
  parms beta1 0 beta2 -36;
  array p[4] py pd pw ph;
  array f[4] fy fd fw fh;
  array alpha[4] alpha1-alpha4;
  alpha4=0;
  bot=0;
  do j=1 to 4;
    bot=bot+exp(beta1*f[j]+beta2*p[j]+alpha[j]);
  end;
  top=exp(beta1*f[brand]+beta2*p[brand]+alpha[brand]);
  lambda=top/bot;
  model choice ~ poisson(lambda);
  random alpha1 alpha2 alpha3
         ~normal([a1,a2,a3],[sig11,sig21,sig31,sig22,sig32,sig33]) subject=id;
run;