The estimation is done with NLMIXED.

Following the link for additional information.
http://support.sas.com/rnd/app/da/new/802ce/ets/chap4/sect13.htm



data tmp;
   do i = 1 to 100;
       x1=rannor(123); x2=rannor(123);x3=rannor(123);
       s=1+1*x1>  rannor(123);
       if s=1 then y= 1+2*x2+ 1.5*rannor(123);
       else  y= 3+1*x3+ 1*rannor(123);
       output;
    end;
run;

   proc nlmixed data=tmp;
      parms a0 a1 b0 b1 c0 c1=0 s1 s2=10  ;
      bounds s1 s2 >0;

      xbeta=a0+a1*x1;
      prob=probnorm(xbeta);

      pdf1=( 1/( s1*sqrt(2*constant('pi')) ) *exp(-0.5*( (y-(b0+b1*x2))/s1) **2) );
      pdf2=( 1/( s2*sqrt(2*constant('pi')) ) *exp(-0.5*( (y-(c0+c1*x3))/s2) **2) );

      ll= min( max(prob *  pdf1 + (1-prob)*pdf2, 1e-20), 1-1e-20);

      loglik=log(ll);

      model y ~ general(loglik);
run;