本文整理汇总了C++中PutRNGstate函数的典型用法代码示例。如果您正苦于以下问题：C++ PutRNGstate函数的具体用法？C++ PutRNGstate怎么用？C++ PutRNGstate使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了PutRNGstate函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: twosample_ks_cm_ad
//.........这里部分代码省略.........
freqg[j] = 0;
for (j=0; j<*n; j++) { /* n times randomly select an index and increase its frequency */
r = (unif_rand()**n);
freqg[r] += 1;
}
/* now take corresponding number of observations from the data */
r = 0;
for (j=0; j<*n; j++) {
for (k=0; k<freqg[j]; k++) { /* put the j-th observation freqg[j]-times to the bootstrap sample */
eventg[r] = event[j];
++r;
}
}
for (first_event_g=0; (event[first_event_g]==0)&&(first_event_g<*n) ; first_event_g++)
;
for (j=0; j<*n; j++) { /* each observation is randomly assigned to one of the groups; group from input is reused */
if (unif_rand() < ((double) y1/(y1+y2))) {
group[j] = 1;
--y1;
} else {
group[j] = 2;
--y2;
}
}
y1 = y1bak;
y2 = y2bak;
ks_cm_ad_process(eventg, group, n, &y1, &y2, rho, gamma, &first_event_g, du_temp, dsigma, ks_b, cm_w, cm_b, ad_w, ad_b);
/* store this bootstrapped process for plotting */
if (i<*nboot_plot) {
for (j=0; j<*n; j++)
du_boot_plot[j+i**n] = du_temp[j];
}
y1 = y1bak;
y2 = y2bak;
ks_cm_ad_stat(n,du_temp,g,ks_b,cm_w,cm_b,ad_w,ad_b,&temp1,&temp2,&temp3,&temp4,&temp5,&temp6);
*pval_ks_w_boot += (temp1>*stat_ks_w);
*pval_ks_b_boot += (temp2>*stat_ks_b);
*pval_cm_w_boot += (temp3>*stat_cm_w);
*pval_cm_b_boot += (temp4>*stat_cm_b);
*pval_ad_w_boot += (temp5>*stat_ad_w);
*pval_ad_b_boot += (temp6>*stat_ad_b);
}
*pval_ks_w_boot /= (double) *nboot;
*pval_ks_b_boot /= (double) *nboot;
*pval_cm_w_boot /= (double) *nboot;
*pval_cm_b_boot /= (double) *nboot;
*pval_ad_w_boot /= (double) *nboot;
*pval_ad_b_boot /= (double) *nboot;
}
/* PERMUTATION tests (g is 1) */
*pval_ks_w_perm = *pval_ks_b_perm = *pval_cm_w_perm = *pval_cm_b_perm = *pval_ad_w_perm = *pval_ad_b_perm = 0.;
for (i=0; i<*n; i++)
g[i] = 1.;
if (*nperm>0) {
for (i=0; i<*nperm; i++) {
for (j=0; j<*n; j++) { /* each observation is randomly assigned to one of the groups; group from input is reused (destroyed) */
if (unif_rand() < ((double) y1/(y1+y2))) {
group[j] = 1;
--y1;
} else {
group[j] = 2;
--y2;
}
}
y1 = y1bak;
y2 = y2bak;
ks_cm_ad_process(event, group, n, &y1, &y2, rho, gamma, &first_event, du_temp, dsigma, ks_b, cm_w, cm_b, ad_w, ad_b);
/* store this permutation process for plotting */
if (i<*nperm_plot) {
for (j=0; j<*n; j++)
du_perm_plot[j+i**n] = du_temp[j];
}
y1 = y1bak;
y2 = y2bak;
ks_cm_ad_stat(n,du_temp,g,ks_b,cm_w,cm_b,ad_w,ad_b,&temp1,&temp2,&temp3,&temp4,&temp5,&temp6);
*pval_ks_w_perm += (temp1>*stat_ks_w);
*pval_ks_b_perm += (temp2>*stat_ks_b);
*pval_cm_w_perm += (temp3>*stat_cm_w);
*pval_cm_b_perm += (temp4>*stat_cm_b);
*pval_ad_w_perm += (temp5>*stat_ad_w);
*pval_ad_b_perm += (temp6>*stat_ad_b);
}
*pval_ks_w_perm /= (double) *nperm;
*pval_ks_b_perm /= (double) *nperm;
*pval_cm_w_perm /= (double) *nperm;
*pval_cm_b_perm /= (double) *nperm;
*pval_ad_w_perm /= (double) *nperm;
*pval_ad_b_perm /= (double) *nperm;
}
/* ASYMPTOTIC p-values of KS */
*pval_ks_w_asympt = 1. - psupw(*stat_ks_w,50);
*pval_ks_b_asympt = 1. - psupb(*stat_ks_b,.5,50);
PutRNGstate();
}示例2: SimOneMVN_mxIW
void SimOneMVN_mxIW(double *nu, double *Lbdinvhlf, double *f1f2, int *pd,
int *pnreps, int *pN, double *es, double *YY)
{
int i, j, k, l, d, d2, N, nreps, mxnreps, Jrand;
int *lbuff;
double xd, sm, tstnu, nu_1, nu_2, zz, lambda, f_1, f_2;
double *df, *pW, *SgmHlf, *Y, *xbuff, *Sigma, *sig, *SigInv;
N = *pN;
d = *pd;
xd = (double) d;
d2 = d*d;
mxnreps=0;
for(l=0;l<N;l++) if(mxnreps < *(pnreps+l)) mxnreps = *(pnreps+l);
lbuff = (int *)S_alloc( 1,sizeof(int));
df = (double *)S_alloc( 1, sizeof(double));
pW = (double *)S_alloc( d2, sizeof(double));
xbuff = (double *)S_alloc( d, sizeof(double));
SgmHlf = (double *)S_alloc( d2, sizeof(double));
Y = (double *)S_alloc(mxnreps*d, sizeof(double));
Sigma = (double *)S_alloc( d2, sizeof(double));
SigInv = (double *)S_alloc( d2, sizeof(double));
sig = (double *)S_alloc( d, sizeof(double));
f_1 = *f1f2;
f_2 = *(f1f2+1);
lambda = *nu/(2 * xd + 2.0);
nu_1 = (2.0 * xd + 2.0)*(1.0 + (lambda - 1.0) *(1.0 - f_1)/(1-f_1/f_2));
nu_2 = (2.0 * xd + 2.0)*(1.0 + (lambda - 1.0) * f_2 );
tstnu = f_1/(nu_2 - (2.0*xd+2.0)) + (1.0-f_1)/(nu_1 - (2.0*xd+2.0));
tstnu = 1.0/tstnu + 2.0*xd + 2.0;
Rprintf("nu_1=%g, nu_2=%g, nu ?= %g", nu_1, nu_2, tstnu);
GetRNGstate();
/* NOTE: */
/* this block computes the average std dev over genes from the model */
/* its diagonal elements, passed to the pointer, sig (of size 3) */
/* are used for the purposes of assigning mean value to Y's under the alternative */
/* */
for(i=0;i<d;i++)
for(j=0;j<d;j++){
sm = 0.0;
for(k=0;k<d;k++)
sm += *(Lbdinvhlf + d*i + k) * (*(Lbdinvhlf + d*j + k));
*(SigInv + d*j + i) = sm;
}
matinv(SigInv, Sigma, pd);
for(i=0;i<d;i++) *(sig + i) = pow((*(Sigma + d*i + i))/(*nu - 2.0*xd - 2.0), 0.5);
for(l=0;l<N;l++){
/* Pick J = 1 w.p. 1-f_1, J= 2 w.p. f_1. */
/* Then draw an InvWish_d(nu_J, Lambda) matrix. This is done */
/* using the result: if Sigma^(-1) ~ Wish_d(nu-d-1, Lambda^(-1)) then */
/* Sigma ~ InvWish_d(nu, Lambda). I simulate N i.i.d. Wish_d(nu-d-1,Lambda^(-1)) */
/* matrices and then invert to get Sigma. One more catch, my rwishart routine */
/* uses the cholesky square root of the parameter matrix, Lambda instead of Lambda */
/* itself. Since I want the parameter matrix in the Wishart to be Lambda^(-1) then */
/* I should pass its cholesky square root which is Lbdinvhlf, e.g. the cholesky */
/* square root of Lambda inverse. That is computed in the calling R script and */
/* passed in. Notice the need to check that Lambda is nonsingular and that */
/* nu > 2*d + 2 (required so that the expected value of the inverse wishart */
/* is finite.) */
/* */
zz = unif_rand();
Jrand = 1*(zz < f_1);
*df = (1-Jrand)*nu_1 + Jrand*nu_2 - xd - 1.0;
rwishart1(df, pd, Lbdinvhlf, pW);
matinv(pW, Sigma, pd);
/* */
/* Sigma ~ (1-f_1)*InvWish_d(nu_1, Lambda) + f_1*InvWish_d(nu_2, Lambda) */
/* */
/* Next, use Sigma to simulate i.i.d. N(0_d, Sigma)'s */
nreps = *(pnreps + l);
*lbuff = nreps*d;
rnormn(lbuff, Y);
chol(Sigma, SgmHlf, pd);
for(i=0;i<nreps;i++){
for(j=0;j<d;j++){
sm = 0.0;
for(k=0;k<d;k++) sm += (*(SgmHlf +d*j +k))*(*(Y +d*i +k));
*(xbuff+j) = sm;
}
for(j=0;j<d;j++) *(Y + d*i + j) = *(xbuff + j) + *(es + l)*(*(sig + j));
}
for(i=0;i<(nreps*d);i++) *(YY + mxnreps*d*l + i) = *(Y+i);
}
PutRNGstate();
}示例3: circemb
//.........这里部分代码省略.........
double u, v;
int card = (i != 0) * (i != halfM) + 2 * (j != 0) * (j != halfM);
switch (card){
case 3:
//B(1) = {1}, B^c(1) = {2}
j1 = (m - i) + m * j;
j2 = i + m * (m - j);
u = norm_rand();
v = norm_rand();
a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
a[j1] *= u;
ia[j1] *= v;
a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
a[j2] *= u;
ia[j2] *= -v;
//B(2) = {1,2}, B^c(2) = {0}
j1 = (m - i) + m * (m - j);
j2 = i + m * j;
u = norm_rand();
v = norm_rand();
a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
a[j1]*= u;
ia[j1] *= v;
a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
a[j2]*= u;
ia[j2] *= -v;
break;
case 1:
//B(1) = 0, B^c(1) = {1}
j1 = i + m * j;
j2 = m - i + m * j;
u = norm_rand();
v = norm_rand();
a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
a[j1] *= u;
ia[j1] *= v;
a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
a[j2] *= u;
ia[j2] *= -v;
break;
case 2:
//B(1) = 0, B^c(1) = {2}
j1 = i + m * j;
j2 = i + m * (m - j);
u = norm_rand();
v = norm_rand();
a[j1] = ia[j1] = M_SQRT1_2 * rho[j1];
a[j1] *= u;
ia[j1] *= v;
a[j2] = ia[j2] = M_SQRT1_2 * rho[j2];
a[j2] *= u;
ia[j2] *= -v;
break;
case 0:
j1 = i + m * j;
a[j1] = rho[j1] * norm_rand();
ia[j1] = 0;
break;
}
}
/* ---------- Computation of Q \Lambda^1/2 Q* Z ------------ */
int maxf, maxp, *iwork;
double *work;
/* The next lines is only valid for 2d random fields. I need to
change if m_1 \neq m_2 as here I suppose that m_1 = m_2 = m */
fft_factor(m, &maxf, &maxp);
work = (double *)R_alloc(4 * maxf, sizeof(double));
iwork = (int *)R_alloc(maxp, sizeof(int));
fft_work(a, ia, m, m, 1, -1, work, iwork);
fft_factor(m, &maxf, &maxp);
work = (double *)R_alloc(4 * maxf, sizeof(double));
iwork = (int *)R_alloc(maxp, sizeof(int));
fft_work(a, ia, 1, m, m, -1, work, iwork);
for (i=nbar;i--;)
ans[i + k * nbar] = isqrtMbar * a[i % *ngrid + m * (i / *ngrid)];
}
PutRNGstate();
if (*nugget > 0){
int dummy = *nsim * nbar;
double sqrtNugget = sqrt(*nugget);
GetRNGstate();
for (i=dummy;i--;)
ans[i] += sqrtNugget * norm_rand();
PutRNGstate();
}
free(a); free(ia);
return;
}示例4: MCMC_wrapper
/*****************
void MCMC_wrapper
Wrapper for a call from R.
and don't forget that tail -> head
*****************/
void MCMC_wrapper(int *dnumnets, int *nedges,
int *tails, int *heads,
int *dn, int *dflag, int *bipartite,
int *nterms, char **funnames,
char **sonames,
char **MHproposaltype, char **MHproposalpackage,
double *inputs, double *theta0, int *samplesize,
double *sample, int *burnin, int *interval,
int *newnetworktails,
int *newnetworkheads,
int *fVerbose,
int *attribs, int *maxout, int *maxin, int *minout,
int *minin, int *condAllDegExact, int *attriblength,
int *maxedges,
int *status){
int directed_flag;
Vertex n_nodes, nmax, bip;
/* Edge n_networks; */
Network nw[1];
Model *m;
MHproposal MH;
n_nodes = (Vertex)*dn;
/* n_networks = (Edge)*dnumnets; */
nmax = (Edge)abs(*maxedges);
bip = (Vertex)*bipartite;
GetRNGstate(); /* R function enabling uniform RNG */
directed_flag = *dflag;
m=ModelInitialize(*funnames, *sonames, &inputs, *nterms);
/* Form the network */
nw[0]=NetworkInitialize(tails, heads, nedges[0],
n_nodes, directed_flag, bip, 0, 0, NULL);
MH_init(&MH,
*MHproposaltype, *MHproposalpackage,
inputs,
*fVerbose,
nw, attribs, maxout, maxin, minout, minin,
*condAllDegExact, *attriblength);
*status = MCMCSample(&MH,
theta0, sample, *samplesize,
*burnin, *interval,
*fVerbose, nmax, nw, m);
MH_free(&MH);
/* Rprintf("Back! %d %d\n",nw[0].nedges, nmax); */
/* record new generated network to pass back to R */
if(*status == MCMC_OK && *maxedges>0 && newnetworktails && newnetworkheads)
newnetworktails[0]=newnetworkheads[0]=EdgeTree2EdgeList(newnetworktails+1,newnetworkheads+1,nw,nmax-1);
ModelDestroy(m);
NetworkDestroy(nw);
PutRNGstate(); /* Disable RNG before returning */
}示例5: updateIndicesX
//.........这里部分代码省略.........
error("'den' must be of type 'double'.");
int *neighbors = INTEGER(sneighbors);
int nneigh = asInteger(snneigh);
if (nneigh <= 0)
error("The number of neighbors must be positive.");
int k = asInteger(sk);
if (k <= 0)
error("The number of components must be positive.");
int *Z = INTEGER(sZ);
int ldZ = LENGTH(sZ) / k;
if (ldZ <= 0)
error("The leading dimension of 'Z' must be positive.");
int ldN = LENGTH(sneighbors) / nneigh;
if (ldN <= 0)
error("The leading dimension of 'neighbors' must be positive.");
int ldD = LENGTH(sden) / k;
if (ldD <= 0)
error("The leading dimension of 'den' must be positive.");
if (ldZ - 1 != ldN || ldZ - 1 != ldD || ldN != ldD)
error("The leading dimension of 'Z', 'neighbors' and 'den' do not match.");
int ldC = LENGTH(scheck) / k;
double *den = REAL(sden);
double *check = REAL(scheck);
int nblocks = LENGTH(blocks);
int b;
GetRNGstate();
for (b = 0; b < nblocks; b++) {
SEXP spoints = VECTOR_ELT(blocks, b);
int n = LENGTH(spoints);
int *points = INTEGER(spoints);
double *U = (double *) R_alloc(n, sizeof(double));
/*double *prob = (double *) R_alloc(k, sizeof(double));*/
/* int *Ni = (int *) R_alloc(k, sizeof(int));*/
int i;
for (i = 0; i < n; i++)
U[i] = unif_rand();
#pragma omp parallel for firstprivate(k, ldD, ldN, ldZ, ldC, nneigh, neighbors, \
Z, den, check, points, U)
for (i = 0; i < n; i++) {
int j, m;
double s = 0.0;
int number = 0;
double prob[10];
int Ni[10];
/* compute the posterior weights for the different classes */
for (j = 0; j < k; j++) {
Ni[j] = 0;
for (m = 0; m < nneigh; m++) {
int mm = neighbors[points[i] - 1 + m * ldN] - 1;
Ni[j] += Z[mm + j * ldZ];
}
}
/*compute the index to use the check table*/
for (j = 0; j < k; j++) {
number += Ni[j] * pow(nneigh+1, j);
}
for (j = 0; j < k; j++) {
/* use tabled value check[or + j*ldC] = exp(beta * Ni[j]) */
prob[j] = den[points[i] - 1 + j * ldD] * check[number + j * ldC];
s += prob[j];
}
/* fix up the weights and convert to probabilities */
if (s != 0.0 && R_FINITE(s))
for (j = 0; j < k; j++)
prob[j] = prob[j] / s;
else
for (j = 0; j < k; j++)
prob[j] = 1.0 / k;
/* generate new Z[points[i] - 1,] entries */
for (j = 0; j < k; j++)
Z[points[i] - 1 + j * ldZ] = 0;
for (m = 0, s = 0.0, j = 0; j < k - 1; j++) {
s += prob[j];
if (U[i] > s) m++;
}
Z[points[i] - 1 + m * ldZ] = 1;
}
}
PutRNGstate();
return sZ;
}示例6: rgeomdirect
void rgeomdirect(double *coord, int *nObs, int *nSite, int *dim,
int *covmod, int *grid, double *sigma2, double *nugget,
double *range, double *smooth, double *uBound,
double *ans){
/* This function generates random fields for the geometric model
coord: the coordinates of the locations
nObs: the number of observations to be generated
nSite: the number of locations
dim: the random field is generated in R^dim
covmod: the covariance model
grid: Does coord specifies a grid?
sigma2: the variance of the geometric gaussian process
nugget: the nugget parameter
range: the range parameter
smooth: the smooth parameter
ans: the generated random field */
int i, j, neffSite, lagi = 1, lagj = 1, oneInt = 1;
const double loguBound = log(*uBound), halfSigma2 = 0.5 * *sigma2;
double sigma = sqrt(*sigma2), sill = 1 - *nugget;
if (*grid){
neffSite = R_pow_di(*nSite, *dim);
lagi = neffSite;
}
else{
neffSite = *nSite;
lagj = *nObs;
}
double *covmat = malloc(neffSite * neffSite * sizeof(double)),
*gp = malloc(neffSite * sizeof(double));
buildcovmat(nSite, grid, covmod, coord, dim, nugget, &sill, range,
smooth, covmat);
/* Compute the Cholesky decomposition of the covariance matrix */
int info = 0;
F77_CALL(dpotrf)("U", &neffSite, covmat, &neffSite, &info);
if (info != 0)
error("error code %d from Lapack routine '%s'", info, "dpotrf");
GetRNGstate();
for (i=*nObs;i--;){
double poisson = 0;
int nKO = neffSite;
while (nKO) {
/* The stopping rule is reached when nKO = 0 i.e. when each site
satisfies the condition in Eq. (8) of Schlather (2002) */
poisson += exp_rand();
double ipoisson = -log(poisson), thresh = loguBound + ipoisson;
/* We simulate one realisation of a gaussian random field with
the required covariance function */
for (j=neffSite;j--;)
gp[j] = norm_rand();
F77_CALL(dtrmv)("U", "T", "N", &neffSite, covmat, &neffSite, gp, &oneInt);
nKO = neffSite;
double ipoissonMinusHalfSigma2 = ipoisson - halfSigma2;
for (j=neffSite;j--;){
ans[j * lagj + i * lagi] = fmax2(sigma * gp[j] + ipoissonMinusHalfSigma2,
ans[j * lagj + i * lagi]);
nKO -= (thresh <= ans[j * lagj + i * lagi]);
}
}
}
PutRNGstate();
/* So fare we generate a max-stable process with standard Gumbel
margins. Switch to unit Frechet ones */
for (i=*nObs * neffSite;i--;)
ans[i] = exp(ans[i]);
free(covmat); free(gp);
return;
}示例7: rgeomcirc
//.........这里部分代码省略.........
case 1:
whittleMatern(dist, mbar, zero, sill, *range, *smooth, rho);
break;
case 2:
cauchy(dist, mbar, zero, sill, *range, *smooth, rho);
break;
case 3:
powerExp(dist, mbar, zero, sill, *range, *smooth, rho);
break;
case 4:
bessel(dist, mbar, *dim, zero, sill, *range, *smooth, rho);
break;
}
/* Compute the eigen values to check if the circulant embbeding
matrix is positive definite */
/* Note : The next lines is only valid for 2d random fields. I
need to change if there are m_1 \neq m_2 as I suppose that m_1
= m_2 = m */
int maxf, maxp;
fft_factor(m, &maxf, &maxp);
double *work = (double *)R_alloc(4 * maxf, sizeof(double));
int *iwork = (int *)R_alloc(maxp, sizeof(int));
fft_work(rho, irho, m, m, 1, -1, work, iwork);
fft_factor(m, &maxf, &maxp);
work = (double *)R_alloc(4 * maxf, sizeof(double));
iwork = (int *)R_alloc(maxp, sizeof(int));
fft_work(rho, irho, 1, m, m, -1, work, iwork);
//Check if the eigenvalues are all positive
for (i=mbar;i--;){
notPosDef |= (rho[i] <= 0) || (fabs(irho[i]) > 0.001);
}
if (notPosDef){
k++;
m = HCN[k];
halfM = m / 2;
mbar = m * m;
}
if (k > 30)
error("Impossible to embbed the covariance matrix");
} while (notPosDef);
/* --------- end of the embedding stage --------- */
/* Computation of the square root of the eigenvalues */
for (i=mbar;i--;){
rho[i] = sqrt(rho[i]);
irho[i] = 0;//No imaginary part
}
int mdag = m / 2 + 1, mdagbar = mdag * mdag;
double isqrtMbar = 1 / sqrt(mbar);
double *a = (double *)R_alloc(mbar, sizeof(double));
double *ia = (double *)R_alloc(mbar, sizeof(double));
GetRNGstate();
for (i=*nObs;i--;){
int nKO = nbar;
double poisson = 0;
while (nKO) {
/* The stopping rule is reached when nKO = 0 i.e. when each site
satisfies the condition in Eq. (8) of Schlather (2002) */
int j;
double *gp = (double *)R_alloc(nbar, sizeof(double));
poisson += exp_rand();
double ipoisson = -log(poisson), thresh = loguBound + ipoisson;
/* We simulate one realisation of a gaussian random field with
the required covariance function */
circcore(rho, a, ia, m, halfM, mdag, mdagbar, *ngrid, nbar, isqrtMbar, *nugget, gp);
nKO = nbar;
double ipoissonMinusHalfSigma2 = ipoisson - halfSigma2;
for (j=nbar;j--;){
ans[j + i * nbar] = fmax2(sigma * gp[j] + ipoissonMinusHalfSigma2,
ans[j + i * nbar]);
nKO -= (thresh <= ans[j + i * nbar]);
}
}
}
PutRNGstate();
/* So fare we generate a max-stable process with standard Gumbel
margins. Switch to unit Frechet ones */
for (i=*nObs * nbar;i--;)
ans[i] = exp(ans[i]);
return;
}示例8: exitRNGScope
// [[Rcpp::register]]
unsigned long exitRNGScope() {
RNGScopeCounter--;
if (RNGScopeCounter == 0) PutRNGstate();
return RNGScopeCounter ;
}示例9: count_permus_with_at_least_k_unfixed_points
void count_permus_with_at_least_k_unfixed_points ( int*n, int *k, double * res){
GetRNGstate();
Generic gen;
*res = (double)gen.count_permus_with_at_least_k_unfixed_points(*n, *k);
PutRNGstate();
}示例10: sampler
//.........这里部分代码省略.........
// update beta
for(j = 0; j < *q; j++){
beta[j] = alphaLong[j]*ksi[j];
}
//update eta, eta+offset
F77_CALL(dgemm)("N", "N", n, &oneInt, q, one, X, n, beta, q, zero, eta, n);
for(int i=0; i<*n; i++) etaOffset[i] = eta[i] + offset[i];
//update sigma_eps
if(*family == 0){
//resid = y - eta - offset
F77_CALL(dcopy)(n, y, &nrv, resid, &nrv); //resid <- y
F77_CALL(daxpy)(n, minusOne, etaOffset, &nrv, resid, &nrv); //resid <- resid - eta - offset
//rss = resid'resid
rss = F77_CALL(ddot)(n, resid, &oneInt, resid, &oneInt);
//update sigma2
invSigma2 = rgamma(*n/2 + *b1, 1/(rss/2 + *b2));
sqrtInvSigma2 = R_pow(invSigma2, 0.5);
scale[0] = sqrtInvSigma2;
*sigma2 = 1 / invSigma2;
}
if(i >= *burnin){
/* report progress */
if(*verbose){
for(j=0; j<9; j++){
if(i == pcts[j]){
Rprintf(".");
#ifdef Win32
R_FlushConsole();
#endif
break;
}
}
}
/* save samples*/
if(i == keep){
for(j = 0; j < *q; j++){
(betaMatR)[save + j*nSamp] = beta[j];
(ksiMatR)[save + j*nSamp] = ksi[j];
}
for(j=0; j < *p; j++){
(alphaMatR)[save + j*nSamp] = alpha[j];
}
for(j=0; j < *pPen; j++){
(tau2MatR)[save + j*nSamp] = tau2[j];
(gammaMatR)[save + j*nSamp] = gamma[j];
(probV1MatR)[save + j*nSamp] = p1[j];
}
(wMatR)[save] = *w;
(sigma2MatR)[save] = *sigma2;
likMatR[save] = logLik(y, etaOffset, *family, scale, *n);
(logPostMatR)[save] = updateLogPost(y, alpha, varAlpha,
ksi, varKsi, scale, *b1, *b2, gamma, *w, *alphaW, *betaW,
tau2, *a1, *a2, *n, *q, *p, *pPen, pIncluded, qKsiNoUpdate, priorMeanKsi, *family, likMatR[save]);
keep = keep + *thin;
save ++;
R_CheckUserInterrupt();
}
} else {
if(*verbose){
if(i == (*burnin-1)){
Rprintf("b");
#ifdef Win32
R_FlushConsole();
#endif
}
}
}
} /* end for i*/
PutRNGstate();
if(*verbose) Rprintf(".");
if(*family > 0) {
acceptance = 0.0;
for(j=0; j<*blocksAlpha; j++) acceptance += acceptAlpha[j];
acceptance = 0.0;
if(qKsiNoUpdate < *q){
for(j=0; j<*blocksKsi; j++) acceptance += acceptKsi[j];
}
}
Free(etaOffset); Free(XKsiUpdate); Free(XAlpha); Free(resid); Free(eta);
Free(offsetKsi); Free(modeKsi); Free(priorMeanKsi); Free(ksiUpdate);
Free(offsetAlpha);
Free(modeAlpha);
Free(priorMeanAlpha);
Free(varAlpha);
Free(alphaLong);
Free(penAlphaSq);
freeXBlockQR(AlphaBlocks, *blocksAlpha);
if(qKsiNoUpdate < *q) freeXBlockQR(KsiBlocks, *blocksKsi);
Free(p1);
return(R_NilValue);
}/* end sampler ()*/示例11: BweibCorSurvmcmc
//.........这里部分代码省略.........
move = 70; */
if(move == 3)
{
BweibSurv_updateSC(beta, &alpha, &kappa, V, survTime, survEvent, cluster, survCov, c, d);
}
/* updating cluster-specific random effect: V
move = 130;*/
if(move == 4)
{
BweibSurv_updateCP(beta, alpha, kappa, V, zeta, survTime, survEvent, cluster, survCov, n_j, accept_V, mhProp_V_var);
}
/* updating precition parameter of the random effect distribution: zeta
move = 160;*/
if(move == 5)
{
BweibSurv_updateVP(V, &zeta, rho1, rho2);
}
/* */
/* Storing posterior samples */
if( ( (MM+1) % *thin ) == 0 && (MM+1) > (*numReps * *burninPerc))
{
StoreInx = (MM+1)/(*thin)- (*numReps * *burninPerc)/(*thin);
samples_alpha[StoreInx - 1] = alpha;
samples_kappa[StoreInx - 1] = kappa;
if(*p >0)
{
for(j = 0; j < *p; j++) samples_beta[(StoreInx - 1) * (*p) + j] = gsl_vector_get(beta, j);
}
for(j = 0; j < *J; j++) samples_V[(StoreInx - 1) * (*J) + j] = gsl_vector_get(V, j);
samples_zeta[StoreInx - 1] = zeta;
if(MM == (*numReps - 1))
{
/* */
if(*p >0)
{
for(j = 0; j < *p; j++) samples_misc[j] = (int) gsl_vector_get(accept_beta, j);
}
/* */
samples_misc[*p] = accept_alpha;
/* */
for(i = 0; i < *J; i++) samples_misc[*p + 1 + i] = (int) gsl_vector_get(accept_V, i);
}
}
if( ( (MM+1) % 10000 ) == 0)
{
time(&now);
Rprintf("iteration: %d: %s\n", MM+1, ctime(&now));
R_FlushConsole();
R_ProcessEvents();
}
}
PutRNGstate();
return;
}示例12: euler_model_simulator
//.........这里部分代码省略.........
case 2: // fixed step
dt = dtt;
nstep = num_map_steps(t,time[step],dt);
break;
default:
errorcall(R_NilValue,"unrecognized 'method'"); // # nocov
break;
}
for (k = 0; k < nstep; k++) { // loop over Euler steps
// interpolate the covar functions for the covariates
table_lookup(&covariate_table,t,cp);
for (j = 0, pm = ps, xm = xt; j < nreps; j++, pm += npars, xm += nvars) { // loop over replicates
switch (mode) {
case Rfun: // R function
{
double *xp = REAL(xvec);
double *pp = REAL(pvec);
double *tp = REAL(tvec);
double *dtp = REAL(dtvec);
double *ap;
*tp = t;
*dtp = dt;
memcpy(xp,xm,nvars*sizeof(double));
memcpy(pp,pm,npars*sizeof(double));
if (first) {
PROTECT(ans = eval(fcall,rho)); nprotect++; // evaluate the call
if (LENGTH(ans) != nvars) {
errorcall(R_NilValue,"user 'step.fun' returns a vector of %d state variables but %d are expected: compare initial conditions?",
LENGTH(ans),nvars);
}
PROTECT(nm = GET_NAMES(ans)); nprotect++;
use_names = !isNull(nm);
if (use_names) {
posn = INTEGER(PROTECT(matchnames(Snames,nm,"state variables"))); nprotect++;
}
ap = REAL(AS_NUMERIC(ans));
first = 0;
} else {
ap = REAL(AS_NUMERIC(eval(fcall,rho)));
}
if (use_names) {
for (i = 0; i < nvars; i++) xm[posn[i]] = ap[i];
} else {
for (i = 0; i < nvars; i++) xm[i] = ap[i];
}
}
break;
case native: // native code
(*ff)(xm,pm,sidx,pidx,cidx,ncovars,cp,t,dt);
break;
default:
errorcall(R_NilValue,"unrecognized 'mode' %d",mode); // # nocov
break;
}
}
t += dt;
if ((meth == 0) && (k == nstep-2)) { // penultimate step
dt = time[step]-t;
t = time[step]-dt;
}
}
}
}
if (mode==1) {
PutRNGstate();
unset_pomp_userdata();
}
UNPROTECT(nprotect);
return X;
}示例13: classRF
//.........这里部分代码省略.........
(nright[n] - nrightimp[n])) / nout[n];
}
}
if (noutall > 0) {
imprt[m + nclass*mdim] +=
((double)(nrightall - nrightimpall)) / noutall;
impsd[m + nclass*mdim] +=
((double) (nrightall - nrightimpall) *
(nrightall - nrightimpall)) / noutall;
}
}
}
}
/* DO PROXIMITIES */
if (iprox) {
computeProximity(prox, oobprox, nodex, jin, oobpair, near);
/* proximity for test data */
if (*testdat) {
computeProximity(proxts, 0, nodexts, jin, oobpair, ntest);
/* Compute proximity between testset and training set. */
for (n = 0; n < ntest; ++n) {
for (k = 0; k < near; ++k) {
if (nodexts[n] == nodex[k])
proxts[n + ntest * (k+ntest)] += 1.0;
}
}
}
}
if (keepf) idxByNnode += *nrnodes;
if (keepInbag) idxByNsample += nsample0;
}
PutRNGstate();
/* Final processing of variable importance. */
for (m = 0; m < mdim; m++) tgini[m] /= Ntree;
if (imp) {
for (m = 0; m < mdim; ++m) {
if (localImp) { /* casewise measures */
for (n = 0; n < nsample; ++n) impmat[m + n*mdim] /= out[n];
}
/* class-specific measures */
for (k = 0; k < nclass; ++k) {
av = imprt[m + k*mdim] / Ntree;
impsd[m + k*mdim] =
sqrt(((impsd[m + k*mdim] / Ntree) - av*av) / Ntree);
imprt[m + k*mdim] = av;
/* imprt[m + k*mdim] = (se <= 0.0) ? -1000.0 - av : av / se; */
}
/* overall measures */
av = imprt[m + nclass*mdim] / Ntree;
impsd[m + nclass*mdim] =
sqrt(((impsd[m + nclass*mdim] / Ntree) - av*av) / Ntree);
imprt[m + nclass*mdim] = av;
imprt[m + (nclass+1)*mdim] = tgini[m];
}
} else {
for (m = 0; m < mdim; ++m) imprt[m] = tgini[m];
}
/* PROXIMITY DATA ++++++++++++++++++++++++++++++++*/
if (iprox) {
for (n = 0; n < near; ++n) {示例14: SummStats
void SummStats(Edge n_edges, Vertex n_nodes, Vertex *tails, Vertex *heads,
Network *nwp, Model *m, double *stats, Network *y0){
Vertex *nodes; /*may need to change*/
nodes = (Vertex *)malloc( 1 * sizeof(Vertex));
GetRNGstate(); /* R function enabling uniform RNG */
ShuffleEdges(tails,heads,n_edges); /* Shuffle edgelist. */
for (unsigned int termi=0; termi < m->n_terms; termi++)
m->termarray[termi].dstats = m->workspace;
/* Doing this one toggle at a time saves a lot of toggles... */
for(Edge e=0; e<n_edges; e++){
ModelTerm *mtp = m->termarray;
double *statspos=stats;
nodes[0] = 0;
for (unsigned int termi=0; termi < m->n_terms; termi++, mtp++){
if(!mtp->s_func){
(*(mtp->d_func))(1, tails+e, heads+e, nodes,
mtp, nwp, y0); /* Call d_??? function */
for (unsigned int i=0; i < mtp->nstats; i++,statspos++)
*statspos += mtp->dstats[i];
}else statspos += mtp->nstats;
}
ToggleEdge(tails[e],heads[e],nwp);
}
/* Doing this one toggle at a time saves a lot of toggles... */
for(int v=0; v<n_nodes; v++){
if(m->nodalstatus[v]== 2.0){
nodes[0] = (int) (v+1);
ModelTerm *mtp = m->termarray;
double *statspos=stats;
for (unsigned int termi=0; termi < m->n_terms; termi++, mtp++){
if(!mtp->s_func){
(*(mtp->d_func))(1, tails, heads, nodes,
mtp, nwp, y0); /* Call d_??? function */
for (unsigned int i=0; i < mtp->nstats; i++,statspos++)
*statspos += mtp->dstats[i];
}else statspos += mtp->nstats;
}
ToggleNode(v+1,nwp);
}
}
ModelTerm *mtp = m->termarray;
double *dstats = m->workspace;
double *statspos=stats;
for (unsigned int termi=0; termi < m->n_terms; termi++, dstats+=mtp->nstats, mtp++ ){
if(mtp->s_func){
(*(mtp->s_func))(mtp, nwp); /* Call s_??? function */
for (unsigned int i=0; i < mtp->nstats; i++,statspos++)
*statspos = mtp->dstats[i];
}else statspos += mtp->nstats;
}
PutRNGstate();
}示例15: simStahl
void simStahl(int *n_sim, double *nu, double *p, double *L,
int *nxo, double *loc, int *max_nxo,
int *n_bins4start)
{
double **Loc, scale;
double curloc=0.0, u;
double *startprob, step;
int i, j, n_nixo;
/* re-organize loc as a doubly index array */
Loc = (double **)R_alloc(*n_sim, sizeof(double *));
Loc[0] = loc;
for(i=1; i < *n_sim; i++)
Loc[i] = Loc[i-1] + *max_nxo;
GetRNGstate();
if(fabs(*nu - 1.0) < 1e-8) { /* looks like a Poisson model */
for(i=0; i< *n_sim; i++) {
R_CheckUserInterrupt(); /* check for ^C */
nxo[i] = rpois(*L);
if(nxo[i] > *max_nxo)
error("Exceeded maximum number of crossovers.");
for(j=0; j < nxo[i]; j++)
Loc[i][j] = runif(0.0, *L);
}
}
else {
scale = 1.0 / (2.0 * *nu * (1.0 - *p));
/* set up starting distribution */
startprob = (double *)R_alloc(*n_bins4start, sizeof(double));
step = *L/(double)*n_bins4start;
startprob[0] = 2.0*(1.0 - *p)*pgamma(0.5*step, *nu, scale, 0, 0)*step;
for(i=1; i< *n_bins4start; i++) {
R_CheckUserInterrupt(); /* check for ^C */
startprob[i] = startprob[i-1] +
2.0*(1.0 - *p)*pgamma(((double)i+0.5)*step, *nu, scale, 0, 0)*step;
}
for(i=0; i< *n_sim; i++) {
R_CheckUserInterrupt(); /* check for ^C */
nxo[i] = 0;
/* locations of chiasmata from the gamma model */
/* shape = nu, rate = 2*nu*(1-p) [scale = 1/{2*nu*(1-p)}] */
u = unif_rand();
if( u > startprob[*n_bins4start-1] )
curloc = *L+1;
else {
for(j=0; j< *n_bins4start; j++) {
if(u <= startprob[j]) {
curloc = ((double)j+0.5)*step;
if(unif_rand() < 0.5) {
nxo[i] = 1;
Loc[i][0] = curloc;
}
break;
}
}
}
if(curloc < *L) {
while(curloc < *L) {
curloc += rgamma(*nu, scale);
if(curloc < *L && unif_rand() < 0.5) {
if(nxo[i] > *max_nxo)
error("Exceeded maximum number of crossovers.");
Loc[i][nxo[i]] = curloc;
(nxo[i])++;
}
}
}
/* locations of crossovers from the no interference mechanism */
if(*p > 0) {
n_nixo = rpois(*L * *p);
if(n_nixo + nxo[i] > *max_nxo)
error("Exceeded maximum number of crossovers.");
for(j=0; j < n_nixo; j++)
Loc[i][nxo[i]+j] = runif(0.0, *L);
nxo[i] += n_nixo;
}
}
}
/* sort the results */
for(i=0; i< *n_sim; i++)
R_rsort(Loc[i], nxo[i]);
PutRNGstate();
}