void BAFT_LNsurv_update_sigSq(gsl_vector *yL,
                              gsl_vector *yU,
                              gsl_vector *yU_posinf,
                              gsl_vector *c0,
                              gsl_vector *c0_neginf,
                              gsl_matrix *X,
                              gsl_vector *y,
                              gsl_vector *beta,
                              double beta0,
                              double *sigSq,
                              double a_sigSq,
                              double b_sigSq,
                              double sigSq_prop_var,
                              int *accept_sigSq)
{
    int i, u;
    double eta, loglh, loglh_prop, logR, gamma_prop, sigSq_prop;
    double logprior, logprior_prop;
    
    int n = X -> size1;
    gsl_vector *xbeta = gsl_vector_calloc(n);
    
    loglh = 0;
    loglh_prop = 0;
    gamma_prop = rnorm(log(*sigSq), sqrt(sigSq_prop_var));
    sigSq_prop = exp(gamma_prop);
    gsl_blas_dgemv(CblasNoTrans, 1, X, beta, 0, xbeta);
    
    for(i=0;i<n;i++)
    {
        eta = beta0 + gsl_vector_get(xbeta, i);
        if(gsl_vector_get(c0_neginf, i) == 0)
        {
            loglh += dnorm(gsl_vector_get(y, i), eta, sqrt(*sigSq), 1) - pnorm(gsl_vector_get(c0, i), eta, sqrt(*sigSq), 0, 1);
            loglh_prop += dnorm(gsl_vector_get(y, i), eta, sqrt(sigSq_prop), 1) - pnorm(gsl_vector_get(c0, i), eta, sqrt(sigSq_prop), 0, 1);
        }else
        {
            loglh += dnorm(gsl_vector_get(y, i), eta, sqrt(*sigSq), 1);
            loglh_prop += dnorm(gsl_vector_get(y, i), eta, sqrt(sigSq_prop), 1);
        }        
    }
    
    logprior = (-a_sigSq-1)*log(*sigSq)-b_sigSq /(*sigSq);
    logprior_prop = (-a_sigSq-1)*log(sigSq_prop)-b_sigSq/sigSq_prop;
    
    logR = loglh_prop - loglh + logprior_prop - logprior + gamma_prop - log(*sigSq);
    
    u = log(runif(0, 1)) < logR;
    
    if(u == 1)
    {
        *sigSq = sigSq_prop;
        *accept_sigSq += 1;
    }
    
    gsl_vector_free(xbeta);
    return;
}

double Fs0_lower(double q, double a, double w, int K)
{
  double F=0;
  for(int k=K; k>=0; k--) {
    F = F - pnorm((-2*k - 2 + w)*a/sqrt(q),0,1,1,0) + pnorm((-2*k - w)*a/sqrt(q),0,1,1,0);
  }

  return 2*F;
}

  //----------------------------------------------------------------------
  std::pair<double, double> BinomialLogitCltDataImputer::impute_large_sample(
      RNG &rng, double number_of_trials, double number_of_successes,
      double linear_predictor) const {
    double information = 0.0;
    const Vector &mixing_weights(mixture_approximation.weights());
    const Vector &sigma(mixture_approximation.sigma());
    double negative_logit_support = plogis(0, linear_predictor, 1, true);
    double positive_logit_support = plogis(0, linear_predictor, 1, false);
    Vector p0 = mixing_weights / negative_logit_support;
    Vector p1 = mixing_weights / positive_logit_support;
    for (int m = 0; m < mixture_approximation.dim(); ++m) {
      p0[m] *= pnorm(0, linear_predictor, sigma[m], true);
      p1[m] *= pnorm(0, linear_predictor, sigma[m], false);
    }

    // p0 is the probability distribution over the mixture component
    // indicators for the failures.  N0 is the count of the number of
    // failures belonging to each mixture component.
    std::vector<int> N0 =
        rmultinom_mt(rng, number_of_trials - number_of_successes, p0 / sum(p0));

    // p1 is the probability distribution over the mixture component
    // indicators for the successes.  N1 is the count of the number
    // of successes in each mixture component.
    std::vector<int> N1 = rmultinom_mt(rng, number_of_successes, p1 / sum(p1));

    double simulation_mean = 0;
    double simulation_variance = 0;
    for (int m = 0; m < N0.size(); ++m) {
      int total_obs = N0[m] + N1[m];
      if (total_obs == 0) {
        continue;
      }
      double sigsq = square(sigma[m]);
      double sig4 = square(sigsq);
      information += total_obs / sigsq;
      double truncated_normal_mean;
      double truncated_normal_variance;
      double cutpoint = 0;
      if (N0[m] > 0) {
        trun_norm_moments(linear_predictor, sigma[m], cutpoint, false,
                          &truncated_normal_mean, &truncated_normal_variance);
        simulation_mean += N0[m] * truncated_normal_mean / sigsq;
        simulation_variance += N0[m] * truncated_normal_variance / sig4;
      }
      if (N1[m] > 0) {
        trun_norm_moments(linear_predictor, sigma[m], cutpoint, true,
                          &truncated_normal_mean, &truncated_normal_variance);
        simulation_mean += N1[m] * truncated_normal_mean / sigsq;
        simulation_variance += N1[m] * truncated_normal_variance / sig4;
      }
    }
    double information_weighted_sum =
        rnorm_mt(rng, simulation_mean, sqrt(simulation_variance));
    return std::make_pair(information_weighted_sum, information);
  }

gnm_float
dsnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean give_log)
{
    if (shape == 0.)
        return dnorm (x, location, scale, give_log);
    else if (give_log)
        return gnm_log (2.) + dnorm (x, location, scale, TRUE) + pnorm (shape * x, shape * location, scale, TRUE, TRUE);
    else
        return 2 * dnorm (x, location, scale, FALSE) * pnorm (shape * x, location/shape, scale, TRUE, FALSE);
}

  double pig(double x, double mu, double lambda, bool logscale){
    if(x <= 0) return logscale ? negative_infinity() :  0;
    if(mu <= 0) throw_exception<std::runtime_error>("mu <= 0 in pig");
    if(lambda <= 0) throw_exception<std::runtime_error>("lambda <= 0 in pig");

    double rlx = sqrt(lambda/x);
    double xmu = x/mu;
    double ans = pnorm(rlx * (xmu -1)) + exp(2*lambda/mu) * pnorm(-rlx*(xmu + 1));
    return logscale ? log(ans) : ans;
  }

/* Sample from a univariate truncated Normal distribution 
   (truncated both from above and below): choose either inverse cdf
   method or rejection sampling method. For rejection sampling, 
   if the range is too far from mu, it uses standard rejection
   sampling algorithm with exponential envelope function. */ 
double TruncNorm(
		 double lb,  /* lower bound */ 
		 double ub,  /* upper bound */
		 double mu,  /* mean */
		 double var, /* variance */
		 int invcdf  /* use inverse cdf method? */
		 ) {
  
  double z;
  double sigma = sqrt(var);
  double stlb = (lb-mu)/sigma;  /* standardized lower bound */
  double stub = (ub-mu)/sigma;  /* standardized upper bound */
  if(stlb > stub)
    error("TruncNorm: lower bound is greater than upper bound\n");
  if(stlb == stub) {
    warning("TruncNorm: lower bound is equal to upper bound\n");
    return(stlb*sigma + mu);
  }
  if (invcdf) {  /* inverse cdf method */
    z = qnorm(runif(pnorm(stlb, 0, 1, 1, 0), pnorm(stub, 0, 1, 1, 0)),
	      0, 1, 1, 0); 
  }
  else { /* rejection sampling method */
    double tol=2.0;
    double temp, M, u, exp_par;
    int flag=0;  /* 1 if stlb, stub <-tol */
    if(stub<=-tol){
      flag=1;
      temp=stub;
      stub=-stlb;
      stlb=-temp;
    }
    if(stlb>=tol){
      exp_par=stlb;
      while(pexp(stub,1/exp_par,1,0) - pexp(stlb,1/exp_par,1,0) < 0.000001) 
	exp_par/=2.0;
      if(dnorm(stlb,0,1,1) - dexp(stlb,1/exp_par,1) >=
	 dnorm(stub,0,1,1) - dexp(stub,1/exp_par,1)) 
	M=exp(dnorm(stlb,0,1,1) - dexp(stlb,1/exp_par,1));
      else
	M=exp(dnorm(stub,0,1,1) - dexp(stub,1/exp_par,1));
      do{ 
	u=unif_rand();
	z=-log(1-u*(pexp(stub,1/exp_par,1,0)-pexp(stlb,1/exp_par,1,0))
	       -pexp(stlb,1/exp_par,1,0))/exp_par;
      }while(unif_rand() > exp(dnorm(z,0,1,1)-dexp(z,1/exp_par,1))/M );  
      if(flag==1) z=-z;
    } 
    else{ 
      do z=norm_rand();
      while( z<stlb || z>stub ); 
    }
  }
  return(z*sigma + mu); 
}

void BAFT_LNsurv_update_beta0(gsl_vector *yL,
                              gsl_vector *yU,
                              gsl_vector *yU_posinf,
                              gsl_vector *c0,
                              gsl_vector *c0_neginf,
                              gsl_matrix *X,
                              gsl_vector *y,
                              gsl_vector *beta,
                              double *beta0,
                              double sigSq,
                              double beta0_prop_var,
                              int *accept_beta0)
{
    int i, u;
    double eta, eta_prop, loglh, loglh_prop, logR, beta0_prop, logprior, logprior_prop;
    
    int n = X -> size1;
    
    gsl_vector *xbeta = gsl_vector_calloc(n);
    
    loglh = 0;
    loglh_prop = 0;
    beta0_prop = rnorm(*beta0, sqrt(beta0_prop_var));
    gsl_blas_dgemv(CblasNoTrans, 1, X, beta, 0, xbeta);
    
    for(i=0;i<n;i++)
    {
        eta = *beta0 + gsl_vector_get(xbeta, i);
        eta_prop = beta0_prop + gsl_vector_get(xbeta, i);
        if(gsl_vector_get(c0_neginf, i) == 0)
        {
            loglh += dnorm(gsl_vector_get(y, i), eta, sqrt(sigSq), 1) - pnorm(gsl_vector_get(c0, i), eta, sqrt(sigSq), 0, 1);
            loglh_prop += dnorm(gsl_vector_get(y, i), eta_prop, sqrt(sigSq), 1) - pnorm(gsl_vector_get(c0, i), eta_prop, sqrt(sigSq), 0, 1);
        }else
        {
            loglh += dnorm(gsl_vector_get(y, i), eta, sqrt(sigSq), 1);
            loglh_prop += dnorm(gsl_vector_get(y, i), eta_prop, sqrt(sigSq), 1);
        }        
    }
    
    logprior = dnorm(*beta0, 0, pow(10,6)*sqrt(sigSq), 1);
    logprior_prop = dnorm(beta0_prop, 0, pow(10,6)*sqrt(sigSq), 1);
    
    logR = loglh_prop - loglh;
    u = log(runif(0, 1)) < logR;
    if(u == 1)
    {
        *beta0 = beta0_prop;
        *accept_beta0 += 1;
    }
    
    gsl_vector_free(xbeta);
    return;
}

void truncatedRat(double *old, double *sd, double *low, double *high, double *newvalue, double *ratio) {
  double lowlimold, upplimold, lowlimnew, upplimnew, plowold, puppold, plownew, puppnew;
  lowlimold = (*low - *old)/ *sd;
  upplimold = (*high - *old)/ *sd;
  lowlimnew = (*low - *newvalue)/ *sd;
  upplimnew = (*high - *newvalue)/ *sd;
  plowold = pnorm(lowlimold,0.0,1.0,1,0);
  puppold = pnorm(upplimold,0.0,1.0,1,0);
  plownew = pnorm(lowlimnew,0.0,1.0,1,0);
  puppnew = pnorm(upplimnew,0.0,1.0,1,0);
  *ratio = (puppold - plowold)/(puppnew - plownew);
}

gnm_float
dsnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean give_log)
{
	if (gnm_isnan (x) || gnm_isnan (shape) || gnm_isnan (location) || gnm_isnan (scale))
		return gnm_nan;

	if (shape == 0.)
		return dnorm (x, location, scale, give_log);
	else if (give_log)
		return M_LN2gnum + dnorm (x, location, scale, TRUE) + pnorm (shape * x, shape * location, scale, TRUE, TRUE);
	else
		return 2 * dnorm (x, location, scale, FALSE) * pnorm (shape * x, location/shape, scale, TRUE, FALSE);
}

double dcutpoints(const cs *liab, double *yP, int *observed, int start,int finish, double *oldcutpoints, double *newcutpoints, int stcutpoints, int ncutpoints, double sdcp, double sdl)
{
    int i,j,w;
    double llik = 0.0;

    for (j = 2 ; j < (ncutpoints-2); j++){
        llik += log(pnorm(oldcutpoints[stcutpoints+j+1]-oldcutpoints[j], 0.0, sdcp, TRUE,FALSE)-pnorm(newcutpoints[stcutpoints+j-1]-oldcutpoints[j], 0.0, sdcp, TRUE,FALSE));
        llik -= log(pnorm(newcutpoints[stcutpoints+j+1]-newcutpoints[j], 0.0, sdcp, TRUE,FALSE)-pnorm(oldcutpoints[stcutpoints+j-1]-newcutpoints[j], 0.0, sdcp, TRUE,FALSE));
    }

    llik += log(1.0-pnorm(newcutpoints[stcutpoints+ncutpoints-3]-oldcutpoints[stcutpoints+ncutpoints-2], 0.0, sdcp, TRUE,FALSE));
    llik -= log(1.0-pnorm(oldcutpoints[stcutpoints+ncutpoints-3]-newcutpoints[stcutpoints+ncutpoints-2], 0.0, sdcp, TRUE,FALSE));

    for (i = start ; i < finish; i++){
        w = yP[i];
        if(w>1 && observed[i]==1){
          if(w==(ncutpoints-1)){
            llik += log(1.0-pnorm(newcutpoints[stcutpoints+w-1], liab->x[i], sdl, TRUE,FALSE));
            llik -= log(1.0-pnorm(oldcutpoints[stcutpoints+w-1], liab->x[i], sdl, TRUE,FALSE));
          }else{
            llik += log(pnorm(newcutpoints[stcutpoints+w], liab->x[i], sdl, TRUE,FALSE)-pnorm(newcutpoints[stcutpoints+w-1], liab->x[i], sdl, TRUE,FALSE));
            llik -= log(pnorm(oldcutpoints[stcutpoints+w], liab->x[i], sdl, TRUE,FALSE)-pnorm(oldcutpoints[stcutpoints+w-1], liab->x[i], sdl, TRUE,FALSE));
          }
        }
    }
    return llik;
}

void _sir_binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
			  int *obsindex, int *stateindex, int *parindex, int *covindex,
			  int ncovars, double *covars, double t) {
  double mean, sd;
  double f;
  mean = CASE*RHO;
  sd = sqrt(CASE*RHO*(1-RHO));
  if (REPORTS > 0) {
    f = pnorm(REPORTS+0.5,mean,sd,1,0)-pnorm(REPORTS-0.5,mean,sd,1,0);
  } else {
    f = pnorm(REPORTS+0.5,mean,sd,1,0);
  }
  *lik = (give_log) ? log(f) : f;
}

gnm_float
psnorm (gnm_float x, gnm_float shape, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p)
{
	gnm_float result, h;

	if (gnm_isnan (x) || gnm_isnan (shape) ||
	    gnm_isnan (location) || gnm_isnan (scale))
		return gnm_nan;

	if (shape == 0.)
		return pnorm (x, location, scale, lower_tail, log_p);

	/* Normalize */
	h = (x - location) / scale;

	/* Flip to a lower-tail problem.  */
	if (!lower_tail) {
		h = -h;
		shape = -shape;
		lower_tail = !lower_tail;
	}

	if (gnm_abs (shape) < 10) {
		gnm_float s = pnorm (h, 0, 1, lower_tail, FALSE);
		gnm_float t = 2 * gnm_owent (h, shape);
		result = s - t;
	} else {
		/*
		 * Make use of this result for Owen's T:
		 *
		 * T(h,a) = .5N(h) + .5N(ha) - N(h)N(ha) - T(ha,1/a)
		 */
		gnm_float s = pnorm (h * shape, 0, 1, TRUE, FALSE);
		gnm_float u = gnm_erf (h / M_SQRT2gnum);
		gnm_float t = 2 * gnm_owent (h * shape, 1 / shape);
		result = s * u + t;
	}

	/*
	 * Negatives can occur due to rounding errors and hopefully for no
	 * other reason.
	 */
	result= CLAMP (result, 0.0, 1.0);

	if (log_p)
		return gnm_log (result);
	else
		return result;
}

double p_swald(double t, double alpha, double nu, double theta, int lower_tail, int log_p)
{
    double p;
    double x;

    if(log_p)
      x = exp(t);
    else
      x = t;

    p = pnorm((nu*(x-theta)-alpha) / sqrt((x-theta)), 0,1,1,0) + 
        exp(2*alpha*nu) * pnorm(-(nu*(x-theta)+alpha) / sqrt((x-theta)), 0,1,1,0);
    
    return (lower_tail ? p : 1-p);
}

// Get the angle between two vectors
double angle(const Vector& u, const Vector& w)
{
  // Get the magnitudes of the vectors
  double unorm = pnorm(u);
  double wnorm = pnorm(w);
  // Get the dot product
  double dprod = inner(u, w);
  // Use the cosine rule
  // but make sure neither is a zero vector
  double rval = 0.0;
  if(dprod > 1E-12){
    rval = std::acos(dprod/(unorm*wnorm));
  }
  return rval;
}

double exp_pnorm(double a, double b)
{
  double r=0;
  if (R_IsNaN(r) && b < -5.5) r = 1/sqrt(2) * exp(a - b*b/2) * (0.5641882/b/b/b - 1/b/sqrt(M_PI));
  else r = exp(a) * pnorm(b,0,1,1,0);
  return r;
}

void betaHyperObjectiveGr(int n, double * par, double * gr, void * ex) {
    // m is location parameter, tau is log of precision parameter
    const double m(par[0]), tau(par[1]);
    // extract objective parameters
    double * input = static_cast<double *>(ex);
    const double beta1_sum(input[0]);
    const double beta1_sqr_sum(input[1]);
    const double P(input[2]);

    const double l1(input[3]);
    const double s1(input[4]);
    const double m_bar(input[5]);
    const double nu_m(input[6]);

    double g_m, g_tau;
    const double log_inv_mills = dnorm(m * exp(.5 * tau), 0, 1, /* give_log */ 1) -
            pnorm(m * exp(.5 * tau), 0, 1, /* lower_tail */ 1, /* give_log */ 1);
    g_m = - P * exp(.5 * tau + log_inv_mills);
    g_m += exp(tau) * (beta1_sum - m * P);
    g_m += nu_m * (m_bar - m);

    g_tau = - P * m * .5 * exp(.5 * tau + log_inv_mills);
    g_tau += - exp(tau) * .5 * (beta1_sqr_sum - 2.0 * m * beta1_sum + P * m * m);
    g_tau += - exp(tau) * l1 + (s1 - 1.0) + P * .5;

    gr[0] = - g_m;
    gr[1] = - g_tau;
}

// Compute the gradient vector of the conditional log likelihood for a Gaussian-Binary model :
void Grad_Cond_Bin(double rho,double pij, double p,int *flag, double *gradcor, double *grad,
		   int *npar, double *nuis, double *thr, double u, double v)
{
  // Initialization variables:
  double dpij=0.0, dij=0.0, rvar=0.0, dpdm=0.0, f=0.0;
  double q1=0.0, q2=0.0, q3=0.0, sh=0.0, vario=0.0, z=0;
  int h=0, i=0, j=0;
  //init variables:
  z=(nuis[0]-*thr)/sqrt(nuis[2]+nuis[1]);
  rvar=nuis[2]/(nuis[2]+nuis[1]);
  //set derivatives components:
  q1=dnorm(z,0,1,0);//stand normal pdf
  q2=pnorm(z*sqrt((1-rvar*rho)/(1+rvar*rho)),0,1,1,0);// stand norm cdf
  q3=d2norm(z,z,rvar*rho);// biv stand norm pdf
  //derivatives:
  dpdm=q1/sqrt(nuis[2]+nuis[1]);/*dp/dmu*/
  dpij=2*dpdm*q2;/*dpij/dmu*/
  f=-(0.5*(nuis[0]-*thr)*dpdm)/(nuis[2]+nuis[1]);/* dp/dsill*/
  dij=2*f*q2;/* dpij/dsill*/
  vario=2*(p-pij);//variogramma binario!!!
  sh=1/(1-2*p+pij);
  // Derivative of the difference respect with the mean
  if(flag[0]==1) { grad[i]=(dpij-2*dpdm)*(1-((u+v)*nij(dpij,dpdm,pij,p)+
					   (u*v)*mij(dpij,dpdm,pij,p)))*sh+dpdm*(1-(u+v)/(2*p))/(1-p); i++; }
  // Derivative of the difference respect with the nugget
  if(flag[1]==1) { grad[i]=1; i++; }
  // Derivative of the difference respect with the sill
  if(flag[2]==1) { grad[i]=(dij-2*f)*(1-((u+v)*nij(dij,f,pij,p)+
					 (u*v)*mij(dij,f,pij,p)))*sh+f*(1-(u+v)/(2*p))/(1-p); i++; }
  // Derivatives with respect to the correlation parameters
  for(j=i;j<*npar;j++) { grad[j]=gradcor[j]*q3*rvar*(1-((u+v)*2*(p-1)/vario +
                                                (u*v)*2*(pij-2*pow(p,2)+p)/(vario*pij)))*sh; h++; }

  return;
}

SEXP pvaluecombine( SEXP RpVec, SEXP Rmethod ) {
	int k = length(RpVec);
	const char * method = CHAR(STRING_ELT(Rmethod, 0));
	
	SEXP Rcmbdpvalue = PROTECT(allocVector(REALSXP, 1));
	memset(REAL(Rcmbdpvalue), 0.0, sizeof(double));
	
	double * cmbdpvalue = REAL(Rcmbdpvalue);
	if (!strcmp(method, "fisher")) {
		for (int i=0; i<k; i++) {
			*cmbdpvalue += log(REAL(RpVec)[i]);
		}
		*cmbdpvalue = 1 - pchisq(-2 * *cmbdpvalue, 2*k, 1, 0);
	} else if (!strcmp(method, "normal") || !strcmp(method, "stouffer")) {
		for (int i=0; i<k; i++) {
			*cmbdpvalue += qnorm(REAL(RpVec)[i], 0.0, 1.0, 1, 0);
		}
		*cmbdpvalue = *cmbdpvalue / sqrt(k);
		*cmbdpvalue = pnorm(*cmbdpvalue, 0.0, 1.0, 1, 0);
	} else if (!strcmp(method, "min") || !strcmp(method, "tippett")) {
		*cmbdpvalue = REAL(RpVec)[0];
		for (int i=1; i<k; i++) {
			*cmbdpvalue = fmin2(*cmbdpvalue, REAL(RpVec)[i]);
		}
		*cmbdpvalue = 1 - pow(1-*cmbdpvalue, k);
	} else if (!strcmp(method, "max")) {
		*cmbdpvalue = REAL(RpVec)[0];
		for (int i=1; i<k; i++) {
			*cmbdpvalue = fmax2(*cmbdpvalue, REAL(RpVec)[i]);
		}
		*cmbdpvalue = pow(*cmbdpvalue, k);
	} else if (!strcmp(method, "sum")) {
		for (int i=0; i<k; i++) {
			*cmbdpvalue += REAL(RpVec)[i];
		}
		if (k <= 30) {
			*cmbdpvalue = pConvolveUniform(*cmbdpvalue, (double)k);
		} else {
			*cmbdpvalue = pnorm(*cmbdpvalue, (double)k/2.0, sqrt((double)k/12.0), 1, 0);
		}
	} else {
		*cmbdpvalue = 3.1415926;
	}
	// return
	UNPROTECT(1);
	return(Rcmbdpvalue);
}

SEXP bernoulliprobrandom(SEXP patterns, SEXP outcomex,SEXP lambdacoef, 
	SEXP gh, SEXP momentdata, SEXP probit)
{
	SEXP ans;
	int irow, outcome, index, noutcomes, nrows, ipoint, npoints, level2size, ilambda, lprobit, *rpatterns = INTEGER(patterns);
	double  *routcomex = REAL(outcomex), *rans,
		neww,newp, *rmomentdata=REAL(momentdata),
		*rgh=REAL(gh),*rlambdacoef=REAL(lambdacoef);
	double product, sum, myoutcomex, myoutcomep;
	
	lprobit = asLogical(probit);
	
	noutcomes = LENGTH(outcomex);
	nrows = LENGTH(patterns)/noutcomes;
	npoints = LENGTH(gh)/2;
	level2size=LENGTH(lambdacoef);
	
	PROTECT(ans = allocVector(REALSXP,nrows));
	
	rans = REAL(ans);
	
	
	for (irow=0; irow < nrows; irow++) {
		/* Rprintf("irow  %d\n",irow); */
		sum=0.0;
/* calculate transformed w and p */
		for (ipoint=0; ipoint < npoints; ipoint++) {
			/* Rprintf("momentdata  %f,%f\n",rmomentdata[irow],rmomentdata[nrows+irow]); */
			newp = rmomentdata[irow]+rmomentdata[nrows+irow]*rgh[ipoint];
			neww = log(rmomentdata[nrows+irow])+
				(rgh[ipoint]*rgh[ipoint])/2.0+log(rgh[npoints+ipoint])-
				newp*newp/2.0;
			/* Rprintf("newp,neww  %f,%f\n",newp,neww); */
			ilambda=0;
			product=1.0;
			for (outcome=0; outcome <noutcomes; outcome++) {
				/* calculate outcome probability for this outcome */
				myoutcomex = routcomex[outcome]+
					rlambdacoef[ilambda]*newp;
				if (lprobit)
					myoutcomep=pnorm(myoutcomex,0,1,TRUE,FALSE);
				else
					myoutcomep=1.0/(1+exp(-myoutcomex));
				ilambda=(ilambda+1) % level2size;				
				/* update likelihood for this observation */
			/*  Rprintf("myoutcomep  %f\n",myoutcomep); */
				index = irow+outcome*nrows;
				if (rpatterns[index]!=NA_INTEGER) {
				  if (rpatterns[index]==1) product = product*myoutcomep;
				  else product = product*(1-myoutcomep); 
				}
			}
			sum=sum+product*exp(neww);
		}
		rans[irow]=sum;
	}
	UNPROTECT(1);		
	return ans;
}

void C_pLausen94_all(const double *Q, double N, const double *m, int N_m, double *pval)
{
	int i,j;
	double *m1 = Calloc(N_m, double);
	double *m2 = Calloc(N_m, double);
	double *T = Calloc(N_m-1, double);
	if(N_m < 2){
		m1[0] = m[0];
		m2[0] = m[0];
		N_m = 1;
	}
	else{
		for(i = 0; i < N_m-1; i++){
			m1[i] = m[i];
			m2[i] = m[i+1];
		}
	}
	/* compute t and D */
	for(j = 0; j < N_m; j++){
		pval[j] = 0.0;
		double D = 0.0;
		for(i = 0; i < N_m-1; i++){
			T[i] = sqrt(1.0-(m1[i]*(N-m2[i]))/((N-m1[i])*m2[i]));
			D += (M_1_PI)*exp(-(pow(Q[j],2))/2)*(T[i] - ((pow(Q[j],2))/4 -1)*(pow(T[i],3))/6);
		}
		pval[j] = 2.0 * (1.0 - pnorm(Q[j], 0.0, 1.0, 1, 0)) + D;
		if(pval[j] > 1.0){
			pval[j] = 1.0;
		}
		if(pval[j] <= FLT_EPSILON){
			pval[j] = 0.0;
		}
		//*pval[j] = 1.0 - pval[j];
	}
	if(N_m - 1 < 1){
		pval[0] = 2.0 * (1.0 - pnorm(Q[0], 0.0, 1.0, 1, 0));
		if(pval[0] > 1.0){
			pval[0] = 1.0;
		}
		if(pval[0] <= FLT_EPSILON){
			pval[0] = 0.0;
		}
		//*pval[0] = 1.0 - pval[0];
	}
	Free(m1);Free(m2);Free(T);
}