C++ rnorm函数代码示例

int test_spd_orderings(taucs_ccs_matrix* A,
		       double* x, double* y, double* b, double* z)
{
  int rc;
  char* metis[] = {"taucs.factor.ordering=metis", NULL};
  char* genmmd[] = {"taucs.factor.ordering=genmmd", NULL};
  char* colamd[] = {"taucs.factor.ordering=colamd", NULL};
  char* amd[] = {"taucs.factor.ordering=amd", NULL};
  void* opt_arg[] = { NULL };
  
  rc = taucs_factor_solve(A,NULL,1, y,b,metis,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  rc = taucs_factor_solve(A,NULL,1, y,b,genmmd,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  rc = taucs_factor_solve(A,NULL,1, y,b,amd,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  /* colamd should fail on symmetric matrices */
  rc = taucs_factor_solve(A,NULL,1, y,b,colamd,opt_arg);
  if (rc == TAUCS_SUCCESS) return TAUCS_ERROR; 

  printf("TESING SYMMETRIC ORDERINGS SUCCEDDED\n");

  return TAUCS_SUCCESS;
}

int test_spd_factorizations(taucs_ccs_matrix* A, 
			    double* x, double* y, double* b, double* z)
{
  int rc;
  char* mf[]   = {"taucs.factor.mf=true", NULL};
  char* ll[]   = {"taucs.factor.ll=true", NULL};
  char* mfmd[] = {"taucs.factor.mf=true", "taucs.maxdepth=5", NULL};
  char* llmd[] = {"taucs.factor.ll=true", "taucs.maxdepth=5", NULL};
  char* ooc[]  = {"taucs.ooc=true", "taucs.ooc.basename=/tmp/taucs-test", NULL};
  void* opt_arg[] = { NULL };
  
  rc = taucs_factor_solve(A,NULL,1, y,b,ooc,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  rc = taucs_factor_solve(A,NULL,1, y,b,mf,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  rc = taucs_factor_solve(A,NULL,1, y,b,ll,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  /* low depth should fail  */
  rc = taucs_factor_solve(A,NULL,1, y,b,mfmd,opt_arg);
  if (rc == TAUCS_SUCCESS) return TAUCS_ERROR; 

  rc = taucs_factor_solve(A,NULL,1, y,b,llmd,opt_arg);
  if (rc == TAUCS_SUCCESS) return TAUCS_ERROR; 

  printf("TESING SPD FACTORIZATIONS SUCCEDDED\n");

  return TAUCS_SUCCESS;
}

void factor_model_row(double *x, unsigned int row, unsigned int n,
		      unsigned int p, unsigned int n_factors, double sigma)
{
	register unsigned int j, k, l;

	double factor = 0.0;

	l = row;
	for (j = 0; j < p; j++) {
		x[l] = 0.0;
		l += n;
	}

	for (k = 0; k < n_factors; k++) {

		factor = rnorm(0.0, sigma);
		l = row;
		for (j = 0; j < p; j++) {
			x[l] += factor * rnorm(0.0, sigma);
			l += n;
		}

	}

	l = row;

	for (j = 0; j < p; j++) {
		x[l] += rnorm(0.0, sigma);
		l += n;
	}

}

// simple 2D Ornstein-Uhlenbeck process simulation
static void sim_ou2 (double *x1, double *x2,
		     double alpha1, double alpha2, double alpha3, double alpha4, 
		     double sigma1, double sigma2, double sigma3)
{
  double eps[2], xnew[2];

  if (!(R_FINITE(*x1))) return;
  if (!(R_FINITE(*x2))) return;
  if (!(R_FINITE(alpha1))) return;
  if (!(R_FINITE(alpha2))) return;
  if (!(R_FINITE(alpha3))) return;
  if (!(R_FINITE(alpha4))) return;
  if (!(R_FINITE(sigma1))) return;
  if (!(R_FINITE(sigma2))) return;
  if (!(R_FINITE(sigma3))) return;

  eps[0] = rnorm(0,1);
  eps[1] = rnorm(0,1);

  xnew[0] = alpha1*(*x1)+alpha3*(*x2)+sigma1*eps[0];
  xnew[1] = alpha2*(*x1)+alpha4*(*x2)+sigma2*eps[0]+sigma3*eps[1];

  *x1 = xnew[0];
  *x2 = xnew[1];
}

// bivariate normal measurement error simulator
void ou2_rmeasure (double *y, double *x, double *p, 
		   int *obsindex, int *stateindex, int *parindex, int *covindex,
		   int ncovar, double *covar, 
		   double t) 
{
  double sd = fabs(TAU);
  Y1 = rnorm(x[X1],sd);
  Y2 = rnorm(x[X2],sd);
}

 void LocalLinearTrendModule::SimulateData(int time_dimension) {
   trend_.resize(time_dimension);
   double level = initial_level_;
   double slope = initial_slope_;
   for (int i = 0; i < time_dimension; ++i) {
     trend_[i] = level;
     level += slope + rnorm(0, level_sd_);
     slope += rnorm(0, slope_sd_);
   }
 }

// simple 2D Ornstein-Uhlenbeck process simulation
static void sim_ou2 (double *x1, double *x2,
                     double alpha1, double alpha2, double alpha3, double alpha4,
                     double sigma1, double sigma2, double sigma3)
{
    double eps[2], xnew[2];

    eps[0] = rnorm(0,1);
    eps[1] = rnorm(0,1);

    xnew[0] = alpha1*(*x1)+alpha3*(*x2)+sigma1*eps[0];
    xnew[1] = alpha2*(*x1)+alpha4*(*x2)+sigma2*eps[0]+sigma3*eps[1];

    *x1 = xnew[0];
    *x2 = xnew[1];
}

void leftTruncNorm(double *mu, double *sigma2, double *x){
 int check1, check2;
 double alphaStar, u, muMinus, z;
 muMinus = -*mu/sqrt(*sigma2);
 if (muMinus <= 0.0){
  check1 = FALSE;
  while(check1 == FALSE){
   GetRNGstate();
   z = rnorm(0.0,1.0);
   PutRNGstate();
   check1 = (z > muMinus);
  }
 } else {
  alphaStar = 0.5 * (muMinus + sqrt(muMinus * muMinus + 4.0));
  check2 = FALSE;
  while(check2 == FALSE){
   GetRNGstate();
   z = muMinus + rexp(1/alphaStar);
   PutRNGstate();
   GetRNGstate();
   u = runif(0.0,1.0);
   PutRNGstate();
   check2 = (u <= exp(-0.5*(z-alphaStar) * (z-alphaStar)));
  }
 }
 *x = *mu + z * sqrt(*sigma2);
}

 Vec ArModel::simulate(int n, const Vec &y0) const {
   if(y0.size() != number_of_lags()){
     ostringstream err;
     err << "Error in ArModel::simulate." << endl
         << "Initial state value y0 was size " << y0.size()
         << ", but the model has " << number_of_lags() << " lags."
         << endl;
     report_error(err.str());
   }
   const Vec &phi(this->phi());
   std::deque<double> lags(y0.rbegin(), y0.rend());
   Vec ans;
   ans.reserve(n);
   for(int i = 0; i < n; ++i) {
     double mu = 0;
     for(int lag = 0; lag < number_of_lags(); ++lag) {
       mu += phi[lag] * lags[lag];
     }
     double y = rnorm(mu, sigma());
     lags.push_front(y);
     lags.pop_back();
     ans.push_back(y);
   }
   return ans;
 }

static void xLz(float c,float x,float y,float d){
	if(T==MT)
		glTriangle(x,y,x+cos(d+M_PI/64)*c,y+sin(d+M_PI/64)*c,x+cos(d-M_PI/64)*c,y+sin(d-M_PI/64)*c);
	c*=c;
	for(int i=0;i<2;i++)
		if(dst2(x,y,Px[i],Py[i])<c&&fabsf(rnorm(d-dir(x,y,Px[i],Py[i])))<M_PI/64)Ph[i]--;
}

SEXP mutate_constants_normal(SEXP sexp, double p, double mu, double sigma) {
  SEXP c;
  switch (TYPEOF(sexp)) { // switch for speed
  case NILSXP:
    return sexp; // do nothing with nils
  case REALSXP:
    if (unif_rand() < p) { // mutate constant with probability p
      PROTECT(c = allocVector(REALSXP, 1));
      REAL(c)[0] = REAL(sexp)[0] + rnorm(mu, sigma);
      UNPROTECT(1);
      return c;
    } else {
      return sexp;
    }
  case LANGSXP: {
    int function_arity = 0;
    SEXP tail_e, e;
    PROTECT(tail_e = R_NilValue);
    for (SEXP iterator = CDR(sexp); !isNull(iterator); iterator = CDR(iterator)) { // recurse on actual parameters
      function_arity++; // determine arity on the fly
      SEXP mutated_parameter;
      PROTECT(mutated_parameter = mutate_constants_normal(CAR(iterator), p, mu, sigma));
      PROTECT(tail_e = CONS(mutated_parameter, tail_e));
    }
    PROTECT(e = LCONS(CAR(sexp), tail_e));
    UNPROTECT(2 * function_arity + 2);
    return e;
  }
  case LISTSXP:
    error("mutate_constants_normal: unexpected LISTSXP");
  default: // base case
    return sexp; // do nothing
  }
}

int test_spd_factorsolve(taucs_ccs_matrix* A, 
			 double* x, double* y, double* b, double* z)
{
  int rc;
  void* F = NULL;
  char* factor[] = {"taucs.solve=false", NULL};
  char* solve [] = {"taucs.factor=false", NULL};
  void* opt_arg[] = { NULL };
  
  /* solve without a factorization should fail */
  rc = taucs_factor_solve(A,NULL,1, y,b,solve,opt_arg);
  if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;

  /* solve without a factorization should fail */
  rc = taucs_factor_solve(A,&F,1, y,b,solve,opt_arg);
  if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;

  rc = taucs_factor_solve(A,&F,1, y,b,factor,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;

  rc = taucs_factor_solve(A,&F,1, y,b,solve,opt_arg);
  if (rc != TAUCS_SUCCESS) return rc;
  if (rnorm(A,x,y,z)) return TAUCS_ERROR;

  printf("TESING SPD FACTORSOLVE SUCCEDDED\n");

  return TAUCS_SUCCESS;
}

 Vec MVTR::simulate_fake_x()const{
   uint p = xdim();
   Vec x(p);
   x[0] = 1.0;
   for(uint i=0; i<p; ++i) x[i] = rnorm();
   return x;
 }

/**
 * Simulate beta using the naive Gibbs update
 *
 * @param da an SEXP struct
 *
 */
static void sim_beta(SEXP da){
  int *dm = DIMS_SLOT(da), *k = K_SLOT(da);
  int nB = dm[nB_POS];
  double *beta = FIXEF_SLOT(da), *mh_sd = MHSD_SLOT(da), *l = CLLIK_SLOT(da), 
    *pm = PBM_SLOT(da), *pv = PBV_SLOT(da), *acc = ACC_SLOT(da);
  double xo, xn, l1, l2, A;

  /* initialize llik_mu*/
  *l = llik_mu(da);
  for (int j = 0; j < nB; j++){
    *k = j;
    xo = beta[j];
    xn = rnorm(xo, mh_sd[j]);
    l1 = *l;
    l2 = post_betak(xn, da);
    A = exp(l2 - l1 + 0.5 * (xo - pm[j]) * (xo - pm[j]) / pv[j]);
    /* determine whether to accept the sample */
    if (A < 1 && runif(0, 1) >= A){ /* not accepted */
      *l = l1;       /* revert the likelihood (this is updated in post_betak) */
    }
    else {
      beta[j] = xn;
      acc[j]++;    
    }
  }                  /* update the mean using the new beta */                    
  if (dm[nU_POS]) cpglmm_fitted(beta, 1, da);
  else cpglm_fitted(beta, da);  
}

double CGaussianMDP::sample_phi0
(
   int n,
   double *ysamp,
   double *sigmasamp,
   double s2,
   double m
)
{
   int i;
   double s=0.0;
   double ys=0.0;
   double var;
   double mn;
   
   for(i=0; i<n; i++)
   {
      s += 1.0/sigmasamp[i];
      ys += ysamp[i]/sigmasamp[i];
   }
   
   var = 1.0/(1.0/s2 + s);
   mn = (m/s2 + ys)*var;
   return rnorm(mn,sqrt(var));
}