proc iml;
  start f(z) global(t,x,count);
    count=count+1;
    p=normpdf(t);
        q=normcdf(t+x*z)-normcdf(t);
        m=p*(q**(n-1));
        d=(z**(v-1))*m*exp(-0.5*v*(z**2));
        return(d);
  finish;
  start g(i) global(t,eps);
    interval=.M||.M;
        t=i;
        call quad(pm,"F",interval) eps=eps;
        return(pm);
  finish;
  start h(xx,n,v) global(x,eps);
    x=xx;
        interval=0||.M;
    call quad(p,"G",interval) eps=eps;
    l=n*(v**(v/2));
        k=gamma(v/2)*(2**(v/2-1));
        a=(l/k)*p;
        return(a);
  finish;
  count=0;
  eps=1E-11;
  b=h(0,16,5);
  print p count;
  quit;