C++ pnorm函数代码示例

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;
}