本文整理汇总了C++中dnorm函数的典型用法代码示例。如果您正苦于以下问题：C++ dnorm函数的具体用法？C++ dnorm怎么用？C++ dnorm使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了dnorm函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: migC
//.........这里部分代码省略.........
for (j = 0; j < length_destinations ; j++) {
//mat_distsL[i*length_destinations + j] = fabs(space[departures2[i]-1] - space[destinations[j]-1]);//-*-
//mat_distsl[i*length_destinations + j] = fabs(space[(*space_size+departures2[i]-1)] - space[(*space_size+destinations[j]-1)]);//-*-
//mat_distsL[i*length_destinations + j] = fabs(space[departures[i]-1] - space[destinations[j]-1]); // euclidean
//mat_distsl[i*length_destinations + j] = fabs(space[(*space_size+departures[i]-1)] - space[(*space_size+destinations[j]-1)]); // euclidean
mat_distsL[i*length_destinations + j] = distkm( space[departures[i]-1], space[destinations[j]-1],
(space[(*space_size+departures[i]-1)]+space[(*space_size+destinations[j]-1)])/2, (space[(*space_size+departures[i]-1)]+space[(*space_size+destinations[j]-1)])/2 ); // km
mat_distsl[i*length_destinations + j] = distkm( (space[departures[i]-1]+space[destinations[j]-1])/2, (space[departures[i]-1]+space[destinations[j]-1])/2,
space[(*space_size+departures[i]-1)], space[(*space_size+destinations[j]-1)] ); // km
}
}
/*printf("\n\nmat_distsL\n");
for (i = 0; i < dim_mat_dists; i++) {
printf("%lf , ",mat_distsL[i]);
}
printf("\n\nmat_distsl\n");
for (i = 0; i < dim_mat_dists; i++) {
printf("%lf , ",mat_distsl[i]);
}*/
// ##### 4 ##### //
// Density matrix //
double mean = 0, sdL, sdl ;
sdL = sqrt(sigma[0]);
sdl = sqrt(sigma[1]);
int b_log = 1;
double *densitymat= NULL;
densitymat = malloc(dim_mat_dists * sizeof(double));
for (i = 0; i < length_destinations; i++) {
for (j = 0; j < length_departures; j++) {
//for (j = 0; j < length_departures2; j++) {//-*-
//densitymat[i + (length_destinations * j)] = ((dnorm(mat_distsL[i + (length_destinations * j)] , mean , sdL , b_log)) + (dnorm(mat_distsl[i + (length_destinations * j)] , mean , sdl , b_log))) - ( log((pnorm(1, mean, sdL,1,0)- pnorm(0, mean, sdL,1,0))) + log((pnorm(1,mean,sdl,1,0) - pnorm(0,mean,sdl,1,0)))) ;
densitymat[i + (length_destinations * j)] = exp( (dnorm(mat_distsL[i + (length_destinations * j)] , mean , sdL , b_log) - dnorm(0 , mean , sdL , b_log)) +
(dnorm(mat_distsl[i + (length_destinations * j)] , mean , sdl , b_log) - dnorm(0 , mean , sdl , b_log)) )
/ length_destinations ;
}
}
/*
if (sum_occupied==1){
printf("\n\nF_{i,j}\n");
for (i = 0; i < dim_mat_dists ; i++) {
printf("%e , ",densitymat[i]);
}
} */
// ##### 5 ##### //
/* printf("\n\noccupied\n");
for ( i = 0; i < (*space_size) ; i++) {
printf("%lf \t",(double)occupied[i]);
}*/
//FILE *fp;
//if (sum_occupied==1){fp = fopen("results.dat", "w");}
// Lambda //
double *L = NULL;
L = malloc(length_destinations * sizeof(double));
for (i = 0 ; i < length_destinations ; i++) {
L[i] = (double)occupied[destinations[(i)]-1];
if (L[i]==0) {
L[i] = 1;
} else if (L[i]>0) {
L[i] = *lambda ;
}
开发者ID:AdrienOliva，项目名称:Phyloland，代码行数:67，代码来源:main.c
示例2: errors
/* # code for computing a hierarchical model, with normally distributed
# level 1 errors (variance known) and level 2 follows a DP
# y[i]: observed datum for obs i
# theta[i]: level 1 mean for obs i
# phi: vector of unique values of theta (i.e., clusters)
# config[i]: cluster label / configuration indicator
####################################################################
*/
HHRESULT CGaussianMDP::sample_config
(
int *&config,
int obs,
double *sigma2,
int n,
double *y,
double *phi,
double alpha
)
{
/*
# config: vector of configuration indicators
# obs: index of observation under study
# sigma2: (known) level 1 variances(
# n: sample size
*/
int i,j,nclus,oldconfig,ind;
int sumconfig = 0;
double sumprob;
double tempphi = 0.0;
HHRESULT hr = HH_OK;
/* get number of configurations/clusters
also set up other things to check */
sumconfig = 0;
nclus = 0; /* number of configurations */
for(i=0; i<n; i++)
{
if(config[i]==config[obs]) sumconfig++;
nclus = imax2(config[i], nclus);
}
/* ## STEP 1: nothing changes if obs under study (obs) has its own
cluster w/prob */
if( (sumconfig == 1) && (runif(0.0,1.0) < (nclus-1.0)/nclus))
{
goto Cleanup;
}
// nconfig counts obs in clusters, current obs not included
for(i=0; i<nclus; i++)
{
nconfig[i] = 0;
}
for(j=0; j<n; j++)
{
nconfig[config[j]-1]++;
}
nconfig[config[obs]-1]--; /* #nclus-star */
/* STEP 2: if there are more than 1 obs in case i's cluster, then: */
if(sumconfig > 1)
{
sumprob = 0;
for(j=0; j<nclus; j++)
{
prob[j] = nconfig[j] *
dnorm(y[obs], phi[j], sqrt(sigma2[obs]), 0);
sumprob += prob[j];
}
prob[nclus] = (alpha/(nclus+1)) *
dnorm(y[obs], phi[nclus], sqrt(sigma2[obs]), 0);
sumprob+=prob[nclus];
if(sumprob==0)
{
for(j=0; j<=nclus; j++) prob[j]=1.0;
}
/* need to add in a sample-type function */
config[obs] = multinomial(nclus+1,prob);
goto Cleanup;
}
/* STEP 3: if there is just one obs in cluster but need to sample new clustr:*/
/* else s(i)=1 and need to sample new cluster */
if(sumconfig==1) /* # s/b unnec line */
{
oldconfig=config[obs];
for(i=0; i<n; i++)
{
if(config[i] > oldconfig) config[i]--;
}
//.........这里部分代码省略.........
开发者ID:gregridgeway，项目名称:hhsim，代码行数:101，代码来源:gaussianmdp.cpp
示例3: diffhfunc_v
void diffhfunc_v(double* u, double* v, int* n, double* param, int* copula, double* out)
{
int j, k=1;
double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t12, t13, t15, t16, t18, t19, t20, t21, t22, t27, t33;
double theta = param[0];
for(j=0;j<*n;j++)
{
if(*copula==0)
{
out[j]=0;
}
else if(*copula==1)
{
t1=qnorm(u[j],0.0,1.0,1,0);
t2=qnorm(v[j],0.0,1.0,1,0);
t3=t1-theta*t2;
t4=1.0-pow(theta,2);
t5=sqrt(t4);
t6=t3/t5;
t7=dnorm(t6,0.0,1.0,0);
t8=sqrt(2.0*pi);
t9=pow(t2,2);
t10=exp(-t9/2.0);
out[j]=t7*t8*(-theta)/t5/t10;
}
else if(*copula==2)
{
diffhfunc_v_tCopula_new(&u[j], &v[j], &k, param, copula, &out[j]);
}
else if(*copula==3)
{
t1 = -theta-1.0;
t2 = pow(v[j],1.0*t1);
t4 = 1/v[j];
t5 = pow(u[j],-1.0*theta);
t6 = pow(v[j],-1.0*theta);
t7 = t5+t6-1.0;
t9 = -1.0-1/theta;
t10 = pow(t7,1.0*t9);
out[j] = t10*t4*t1*t2-1/t7*t4*theta*t6*t9*t10*t2;
}
else if(*copula==4)
{
t3 = log(u[j]);
t4 = pow(-t3,1.0*theta);
t5 = log(v[j]);
t6 = pow(-t5,1.0*theta);
t7 = t4+t6;
t8 = 1/theta;
t9 = pow(t7,1.0*t8);
t10 = t6*t6;
t12 = v[j]*v[j];
t13 = 1/t12;
t15 = t5*t5;
t16 = 1/t15;
t18 = t16/t7;
t19 = exp(-t9);
t20 = t8-1.0;
t21 = pow(t7,1.0*t20);
t22 = t19*t21;
t27 = theta*t13;
t33 = t6*t13;
out[j] = t9*t10*t13*t18*t22-t22*t20*t10*t27*t18-t22*t6*t27*t16+t22*t33/t5+t22*t33*t16;
}
else if(*copula==5)
{
t1 = exp(theta);
t2 = theta*u[j];
t3 = exp(t2);
t6 = theta*v[j];
t8 = exp(t6+t2);
t10 = exp(t6+theta);
t12 = exp(t2+theta);
t13 = pow(t8-t10-t12+t1,2.0);
out[j] = t1*(t3-1.0)/t13*(theta*t8-theta*t10);
}
else if(*copula==6)
{
t2 = pow(1.0-u[j],1.0*theta);
t3 = 1.0-v[j];
t4 = pow(t3,1.0*theta);
t5 = t2*t4;
t6 = t2+t4-t5;
t8 = 1/theta-1.0;
t9 = pow(t6,1.0*t8);
t12 = 1/t3;
t19 = theta-1.0;
t20 = pow(t3,1.0*t19);
t22 = 1.0-t2;
out[j] = t9*t8*(-t4*theta*t12+t5*theta*t12)/t6*t20*t22-t9*t20*t19*t12*t22;
}
}
}
开发者ID:cran，项目名称:VineCopula，代码行数:96，代码来源:hfuncderiv.c
示例4: VB5_dmeasure
void VB5_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)
{
lik = dnorm(Lobs,L,L_sd,give_log);
}
开发者ID:claycressler，项目名称:deb_fitting，代码行数:4，代码来源:VB5.c
示例5: spMisalign
//.........这里部分代码省略.........
F77_NAME(dtrsm)(lside, lower, ntran, nUnit, &N, &P1, &one, C, &N, vU, &N);//L^{-1}[v:U] = [y:X]
F77_NAME(dgemm)(ytran, ntran, &P, &P, &N, &one, &vU[N], &N, &vU[N], &N, &zero, tmp_PP, &P); //U'U
F77_NAME(dpotrf)(lower, &P, tmp_PP, &P, &info); if(info != 0){error("c++ error: dpotrf failed\n");}
for(k = 0; k < P; k++) det += 2*log(tmp_PP[k*P+k]);
F77_NAME(dgemv)(ytran, &N, &P, &one, &vU[N], &N, vU, &incOne, &zero, tmp_P, &incOne); //U'v
F77_NAME(dtrsv)(lower, ntran, nUnit, &P, tmp_PP, &P, tmp_P, &incOne);
Q = pow(F77_NAME(dnrm2)(&N, vU, &incOne),2) - pow(F77_NAME(dnrm2)(&P, tmp_P, &incOne),2) ;
}
//
//priors, jacobian adjustments, and likelihood
//
logPostCand = 0.0;
if(KPriorName == "IW"){
logDetK = 0.0;
SKtrace = 0.0;
for(k = 0; k < m; k++){logDetK += 2*log(A[k*m+k]);}
//jacobian \sum_{i=1}^{m} (m-i+1)*log(a_ii)+log(a_ii)
for(k = 0; k < m; k++){logPostCand += (m-k)*log(A[k*m+k])+log(A[k*m+k]);}
//S*K^-1
F77_NAME(dpotri)(lower, &m, A, &m, &info); if(info != 0){error("c++ error: dpotri failed\n");}
F77_NAME(dsymm)(rside, lower, &m, &m, &one, A, &m, KIW_S, &m, &zero, tmp_mm, &m);
for(k = 0; k < m; k++){SKtrace += tmp_mm[k*m+k];}
logPostCand += -0.5*(KIW_df+m+1)*logDetK - 0.5*SKtrace;
}else{
for(k = 0; k < nLTr; k++){
logPostCand += dnorm(params[AIndx+k], ANormMu[k], sqrt(ANormC[k]), 1);
}
}
if(nugget){
for(k = 0; k < m; k++){
logPostCand += -1.0*(1.0+PsiIGa[k])*log(Psi[k])-PsiIGb[k]/Psi[k]+log(Psi[k]);
}
}
for(k = 0; k < m; k++){
logPostCand += log(phi[k] - phiUnif[k*2]) + log(phiUnif[k*2+1] - phi[k]);
if(covModel == "matern"){
logPostCand += log(nu[k] - nuUnif[k*2]) + log(nuUnif[k*2+1] - nu[k]);
}
}
logPostCand += -0.5*det-0.5*Q;
//
//MH accept/reject
//
logMHRatio = logPostCand - logPostCurrent;
if(runif(0.0,1.0) <= exp(logMHRatio)){
logPostCurrent = logPostCand;
if(amcmc){
accept[j]++;
}else{
accept[0]++;
batchAccept++;
开发者ID:cran，项目名称:spBayes，代码行数:67，代码来源:spMisalign.cpp
示例6: BAFT_LNsurv_update_beta
void BAFT_LNsurv_update_beta(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 beta_prop_var,
gsl_vector *accept_beta)
{
int i, j, u;
double eta, eta_prop, loglh, loglh_prop, logR;
int n = X -> size1;
int p = X -> size2;
gsl_vector *beta_prop = gsl_vector_calloc(p);
gsl_vector *xbeta = gsl_vector_calloc(n);
gsl_vector *xbeta_prop = gsl_vector_calloc(n);
j = (int) runif(0, p);
loglh = 0;
loglh_prop = 0;
gsl_vector_memcpy(beta_prop, beta);
gsl_vector_set(beta_prop, j, rnorm(gsl_vector_get(beta, j), sqrt(beta_prop_var)));
gsl_blas_dgemv(CblasNoTrans, 1, X, beta, 0, xbeta);
gsl_blas_dgemv(CblasNoTrans, 1, X, beta_prop, 0, xbeta_prop);
for(i=0;i<n;i++)
{
eta = beta0 + gsl_vector_get(xbeta, i);
eta_prop = beta0 + gsl_vector_get(xbeta_prop, 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);
}
}
logR = loglh_prop - loglh;
u = log(runif(0, 1)) < logR;
if(u == 1)
{
gsl_vector_memcpy(beta, beta_prop);
gsl_vector_set(accept_beta, j, gsl_vector_get(accept_beta, j) + 1);
}
gsl_vector_free(beta_prop);
gsl_vector_free(xbeta);
gsl_vector_free(xbeta_prop);
return;
}
开发者ID:cran，项目名称:SemiCompRisks，代码行数:61，代码来源:BAFT_LNsurv_Updates.c
示例7: do_constraint
//.........这里部分代码省略.........
/* The position corrections dr due to the constraints */
dsvmul(-lambda*rm*pull->grp[g].invtm, vec, dr[g]);
dsvmul( lambda*rm* pref->invtm, vec,ref_dr);
break;
case epullgPOS:
for(m=0; m<DIM; m++)
{
if (pull->dim[m])
{
lambda = r_ij[g][m] - ref[m];
/* The position corrections dr due to the constraints */
dr[g][m] = -lambda*rm*pull->grp[g].invtm;
ref_dr[m] = lambda*rm*pref->invtm;
}
else
{
dr[g][m] = 0;
ref_dr[m] = 0;
}
}
break;
}
/* DEBUG */
if (debug)
{
j = (PULL_CYL(pull) ? g : 0);
get_pullgrps_dr(pull,pbc,g,t,rinew[g],rjnew[j],-1,tmp);
get_pullgrps_dr(pull,pbc,g,t,dr[g] ,ref_dr ,-1,tmp3);
fprintf(debug,
"Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
rinew[g][0],rinew[g][1],rinew[g][2],
rjnew[j][0],rjnew[j][1],rjnew[j][2], dnorm(tmp));
if (pull->eGeom == epullgPOS)
{
fprintf(debug,
"Pull ref %8.5f %8.5f %8.5f\n",
pgrp->vec[0],pgrp->vec[1],pgrp->vec[2]);
}
else
{
fprintf(debug,
"Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f %8.5f %8.5f\n",
"","","","","","",ref[0],ref[1],ref[2]);
}
fprintf(debug,
"Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
dr[g][0],dr[g][1],dr[g][2],
ref_dr[0],ref_dr[1],ref_dr[2],
dnorm(tmp3));
fprintf(debug,
"Pull cor %10.7f %10.7f %10.7f\n",
dr[g][0],dr[g][1],dr[g][2]);
} /* END DEBUG */
/* Update the COMs with dr */
dvec_inc(rinew[g], dr[g]);
dvec_inc(rjnew[PULL_CYL(pull) ? g : 0],ref_dr);
}
/* Check if all constraints are fullfilled now */
for(g=1; g<1+pull->ngrp; g++)
{
pgrp = &pull->grp[g];
开发者ID:martinhoefling，项目名称:gromacs，代码行数:66，代码来源:pull.c
示例8: do_pull_pot
/* Pulling with a harmonic umbrella potential or constant force */
static void do_pull_pot(int ePull,
t_pull *pull, t_pbc *pbc, double t, real lambda,
real *V, tensor vir, real *dVdl)
{
int g,j,m;
dvec dev;
double ndr,invdr;
real k,dkdl;
t_pullgrp *pgrp;
/* loop over the groups that are being pulled */
*V = 0;
*dVdl = 0;
for(g=1; g<1+pull->ngrp; g++)
{
pgrp = &pull->grp[g];
get_pullgrp_distance(pull,pbc,g,t,pgrp->dr,dev);
k = (1.0 - lambda)*pgrp->k + lambda*pgrp->kB;
dkdl = pgrp->kB - pgrp->k;
switch (pull->eGeom)
{
case epullgDIST:
ndr = dnorm(pgrp->dr);
invdr = 1/ndr;
if (ePull == epullUMBRELLA)
{
pgrp->f_scal = -k*dev[0];
*V += 0.5* k*dsqr(dev[0]);
*dVdl += 0.5*dkdl*dsqr(dev[0]);
}
else
{
pgrp->f_scal = -k;
*V += k*ndr;
*dVdl += dkdl*ndr;
}
for(m=0; m<DIM; m++)
{
pgrp->f[m] = pgrp->f_scal*pgrp->dr[m]*invdr;
}
break;
case epullgDIR:
case epullgDIRPBC:
case epullgCYL:
if (ePull == epullUMBRELLA)
{
pgrp->f_scal = -k*dev[0];
*V += 0.5* k*dsqr(dev[0]);
*dVdl += 0.5*dkdl*dsqr(dev[0]);
}
else
{
ndr = 0;
for(m=0; m<DIM; m++)
{
ndr += pgrp->vec[m]*pgrp->dr[m];
}
pgrp->f_scal = -k;
*V += k*ndr;
*dVdl += dkdl*ndr;
}
for(m=0; m<DIM; m++)
{
pgrp->f[m] = pgrp->f_scal*pgrp->vec[m];
}
break;
case epullgPOS:
for(m=0; m<DIM; m++)
{
if (ePull == epullUMBRELLA)
{
pgrp->f[m] = -k*dev[m];
*V += 0.5* k*dsqr(dev[m]);
*dVdl += 0.5*dkdl*dsqr(dev[m]);
}
else
{
pgrp->f[m] = -k*pull->dim[m];
*V += k*pgrp->dr[m]*pull->dim[m];
*dVdl += dkdl*pgrp->dr[m]*pull->dim[m];
}
}
break;
}
if (vir)
{
/* Add the pull contribution to the virial */
for(j=0; j<DIM; j++)
{
for(m=0;m<DIM;m++)
{
vir[j][m] -= 0.5*pgrp->f[j]*pgrp->dr[m];
}
}
}
}
//.........这里部分代码省略.........
开发者ID:martinhoefling，项目名称:gromacs，代码行数:101，代码来源:pull.c
示例9: dlognorm
Type dlognorm(Type x, Type meanlog, Type sdlog, int give_log=0){
Type logres = dnorm( log(x), meanlog, sdlog, true) - log(x);
if(give_log) return logres; else return exp(logres);
}
开发者ID:merrillrudd，项目名称:2016_Spatio-temporal_models，代码行数:4，代码来源:glm_hw.cpp
示例10: get_pullgrp_distance
void get_pullgrp_distance(t_pull *pull,t_pbc *pbc,int g,double t,
dvec dr,dvec dev)
{
static gmx_bool bWarned=FALSE; /* TODO: this should be fixed for thread-safety,
but is fairly benign */
t_pullgrp *pgrp;
int m;
dvec ref;
double drs,inpr;
pgrp = &pull->grp[g];
get_pullgrp_dr(pull,pbc,g,t,dr);
if (pull->eGeom == epullgPOS)
{
for(m=0; m<DIM; m++)
{
ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
}
}
else
{
ref[0] = pgrp->init[0] + pgrp->rate*t;
}
switch (pull->eGeom)
{
case epullgDIST:
/* Pull along the vector between the com's */
if (ref[0] < 0 && !bWarned)
{
fprintf(stderr,"\nPull reference distance for group %d is negative (%f)\n",g,ref[0]);
bWarned = TRUE;
}
drs = dnorm(dr);
if (drs == 0)
{
/* With no vector we can not determine the direction for the force,
* so we set the force to zero.
*/
dev[0] = 0;
}
else
{
/* Determine the deviation */
dev[0] = drs - ref[0];
}
break;
case epullgDIR:
case epullgDIRPBC:
case epullgCYL:
/* Pull along vec */
inpr = 0;
for(m=0; m<DIM; m++)
{
inpr += pgrp->vec[m]*dr[m];
}
dev[0] = inpr - ref[0];
break;
case epullgPOS:
/* Determine the difference of dr and ref along each dimension */
for(m=0; m<DIM; m++)
{
dev[m] = (dr[m] - ref[m])*pull->dim[m];
}
break;
}
}
开发者ID:martinhoefling，项目名称:gromacs，代码行数:69，代码来源:pull.c
示例11: DATA_STRING
Type objective_function<Type>::operator() ()
{
DATA_STRING(distr);
DATA_INTEGER(n);
Type ans = 0;
if (distr == "norm") {
PARAMETER(mu);
PARAMETER(sd);
vector<Type> x = rnorm(n, mu, sd);
ans -= dnorm(x, mu, sd, true).sum();
}
else if (distr == "gamma") {
PARAMETER(shape);
PARAMETER(scale);
vector<Type> x = rgamma(n, shape, scale);
ans -= dgamma(x, shape, scale, true).sum();
}
else if (distr == "pois") {
PARAMETER(lambda);
vector<Type> x = rpois(n, lambda);
ans -= dpois(x, lambda, true).sum();
}
else if (distr == "compois") {
PARAMETER(mode);
PARAMETER(nu);
vector<Type> x = rcompois(n, mode, nu);
ans -= dcompois(x, mode, nu, true).sum();
}
else if (distr == "compois2") {
PARAMETER(mean);
PARAMETER(nu);
vector<Type> x = rcompois2(n, mean, nu);
ans -= dcompois2(x, mean, nu, true).sum();
}
else if (distr == "nbinom") {
PARAMETER(size);
PARAMETER(prob);
vector<Type> x = rnbinom(n, size, prob);
ans -= dnbinom(x, size, prob, true).sum();
}
else if (distr == "nbinom2") {
PARAMETER(mu);
PARAMETER(var);
vector<Type> x = rnbinom2(n, mu, var);
ans -= dnbinom2(x, mu, var, true).sum();
}
else if (distr == "exp") {
PARAMETER(rate);
vector<Type> x = rexp(n, rate);
ans -= dexp(x, rate, true).sum();
}
else if (distr == "beta") {
PARAMETER(shape1);
PARAMETER(shape2);
vector<Type> x = rbeta(n, shape1, shape2);
ans -= dbeta(x, shape1, shape2, true).sum();
}
else if (distr == "f") {
PARAMETER(df1);
PARAMETER(df2);
vector<Type> x = rf(n, df1, df2);
ans -= df(x, df1, df2, true).sum();
}
else if (distr == "logis") {
PARAMETER(location);
PARAMETER(scale);
vector<Type> x = rlogis(n, location, scale);
ans -= dlogis(x, location, scale, true).sum();
}
else if (distr == "t") {
PARAMETER(df);
vector<Type> x = rt(n, df);
ans -= dt(x, df, true).sum();
}
else if (distr == "weibull") {
PARAMETER(shape);
PARAMETER(scale);
vector<Type> x = rweibull(n, shape, scale);
ans -= dweibull(x, shape, scale, true).sum();
}
else if (distr == "AR1") {
PARAMETER(phi);
vector<Type> x(n);
density::AR1(phi).simulate(x);
ans += density::AR1(phi)(x);
}
else if (distr == "ARk") {
PARAMETER_VECTOR(phi);
vector<Type> x(n);
density::ARk(phi).simulate(x);
ans += density::ARk(phi)(x);
}
else if (distr == "MVNORM") {
PARAMETER(phi);
matrix<Type> Sigma(5,5);
for(int i=0; i<Sigma.rows(); i++)
for(int j=0; j<Sigma.rows(); j++)
Sigma(i,j) = exp( -phi * abs(i - j) );
density::MVNORM_t<Type> nldens = density::MVNORM(Sigma);
//.........这里部分代码省略.........
开发者ID:kaskr，项目名称:adcomp，代码行数:101，代码来源:check_simulations.cpp
示例12: F77_SUB
double F77_SUB(dnrm)(double *x, double *mu, double *sigma, int *give_log)
{
return dnorm(*x, *mu, *sigma, *give_log);
}
开发者ID:Syed-Arshad，项目名称:DPpackage，代码行数:4，代码来源:ToolsRfun.c
示例13: DATA_MATRIX
Type objective_function<Type>::operator() () {
// data:
DATA_MATRIX(x_ij);
DATA_VECTOR(y_i);
DATA_IVECTOR(k_i); // vector of IDs
DATA_INTEGER(n_k); // number of IDs
DATA_INTEGER(n_j); // number of IDs
DATA_VECTOR(b1_cov_re_i); // predictor data for random slope
DATA_VECTOR(sigma1_cov_re_i); // predictor data for random slope
//DATA_VECTOR(sigma2_cov_re_i); // predictor data for random slope
// parameters:
PARAMETER_VECTOR(b_j)
PARAMETER_VECTOR(sigma_j);
PARAMETER(log_b0_sigma);
PARAMETER_VECTOR(b0_k);
PARAMETER(log_b1_sigma);
PARAMETER_VECTOR(b1_k);
PARAMETER(log_sigma0_sigma);
PARAMETER(log_sigma1_sigma);
PARAMETER_VECTOR(sigma0_k);
PARAMETER_VECTOR(sigma1_k);
int n_data = y_i.size(); // get number of data points to loop over
// Linear predictor
vector<Type> linear_predictor_i(n_data);
vector<Type> linear_predictor_sigma_i(n_data);
linear_predictor_i = x_ij*b_j;
linear_predictor_sigma_i = x_ij*sigma_j;
Type nll = 0.0; // initialize negative log likelihood
for(int i = 0; i < n_data; i++){
nll -= dnorm(
y_i(i),
b0_k(k_i(i)) + b1_k(k_i(i)) * b1_cov_re_i(i) +
linear_predictor_i(i),
sqrt(exp(
sigma0_k(k_i(i)) +
sigma1_k(k_i(i)) * sigma1_cov_re_i(i) +
linear_predictor_sigma_i(i))),
true);
}
for(int k = 0; k < n_k; k++){
nll -= dnorm(b0_k(k), Type(0.0), exp(log_b0_sigma), true);
nll -= dnorm(b1_k(k), Type(0.0), exp(log_b1_sigma), true);
nll -= dnorm(sigma0_k(k), Type(0.0), exp(log_sigma0_sigma), true);
nll -= dnorm(sigma1_k(k), Type(0.0), exp(log_sigma1_sigma), true);
//nll -= dnorm(sigma2_k(k), Type(0.0), exp(log_sigma2_sigma), true);
}
// Reporting
Type b0_sigma = exp(log_b0_sigma);
Type b1_sigma = exp(log_b1_sigma);
Type sigma0_sigma = exp(log_sigma0_sigma);
Type sigma1_sigma = exp(log_sigma1_sigma);
//Type sigma2_sigma = exp(log_sigma2_sigma);
vector<Type> b1_b1_k(n_k);
vector<Type> sigma1_sigma1_k(n_k);
for(int k = 0; k < n_k; k++){
// these are fixed-effect slopes + random-effect slopes
b1_b1_k(k) = b_j(n_j) + b1_k(k);
sigma1_sigma1_k(k) = sigma_j(n_j) + sigma1_k(k);
}
REPORT( b0_k );
REPORT( b1_k );
REPORT( b_j );
REPORT( sigma0_k );
REPORT( sigma1_k );
//REPORT( sigma2_k );
REPORT(b0_sigma);
REPORT(b1_sigma);
REPORT(sigma0_sigma);
REPORT(sigma1_sigma);
//REPORT(sigma2_sigma);
REPORT(b1_b1_k);
REPORT(sigma1_sigma1_k);
//ADREPORT( b0_k );
//ADREPORT( b1_k );
//ADREPORT( b_j );
//ADREPORT( sigma0_k );
//ADREPORT( sigma1_k );
//ADREPORT( sigma2_k );
//ADREPORT(b0_sigma);
//ADREPORT(b1_sigma);
//ADREPORT(sigma0_sigma);
//ADREPORT(sigma1_sigma);
//ADREPORT(sigma2_sigma);
//ADREPORT(b1_b1_k);
//ADREPORT(sigma1_sigma1_k);
return nll;
}
开发者ID:NCEAS，项目名称:pfx-commercial，代码行数:100，代码来源:revenue_diff0_nooffset.cpp
示例14: NMix_PredCondDensCDFMarg
/***** ***************************************************************************************** *****/
void
NMix_PredCondDensCDFMarg(double* dens,
double* qdens,
int* err,
const int* calc_dens,
const int* nquant,
const double* qprob,
const int* icond,
const double* y,
const int* p,
const int* n,
const int* chK,
const double* chw,
const double* chmu,
const double* chLi,
const int* M)
{
const char *fname = "NMix_PredCondDensCDFMarg";
*err = 0;
if (*p <= 1){
*err = 1;
error("%s: Dimension must be at least 2.\n", fname);
}
if (*icond < 0 || *icond >= *p){
*err = 1;
error("%s: Incorrect index of the margin by which we condition,\n", fname);
}
/***** Variables which will (repeatedly) be used *****/
/***** ========================================= *****/
int m0, i0, i1, t, i, j;
double dtmp;
double csigma; /* to keep std. deviation of the margin by which we condition */
double cov_m0_icond; /* to keep covariance between two margins */
double mu_cond; /* to keep conditional mean when computing conditional cdf */
double sigma_cond; /* to keep conditional std. deviation when computing conditional cdf */
double *densP;
double *dP;
double y2[2]; /* to keep 2-component vector of grid values */
double mu2[2]; /* to keep 2-component vector of means */
double Li2[3]; /* to keep lower triangle of 2x2 matrix */
double * Li2P;
const int *n0;
const int *K;
const double *w, *mu, *Li;
const double *ycP, *y0P, *y0start;
const double *wP = NULL;
const double *muP = NULL;
const double *LiP = NULL;
const int LTp = (*p * (*p + 1))/2; /** length of lower triangles of covariance matrices **/
const int icdiag = (*icond * (2 * (*p) - (*icond) + 1))/2; /** index of diagonal element for icond margin **/
const int TWO = 2;
double log_dets[2];
log_dets[1] = -TWO * M_LN_SQRT_2PI;
/***** lgrid: Total length of the marginal grids (except the grid by which we condition) *****/
/***** lcgrid: Length of the grid of values by which we condition *****/
/***** ycond: Pointer to the first value by which we condition *****/
/***** ldens: Length of the array dens *****/
/***** ========================================================================================= *****/
int lgrid = 0;
int lcgrid;
const double *ycond;
ycond = y;
n0 = n;
for (m0 = 0; m0 < *icond; m0++){
lgrid += *n0;
ycond += *n0;
n0++;
}
lcgrid = *n0;
n0++;
for (m0 = *icond + 1; m0 < *p; m0++){
lgrid += *n0;
n0++;
}
int ldens = (lgrid + 1) * lcgrid;
/***** Working array *****/
/***** ============= *****/
double *dwork = Calloc(2 + LTp + lcgrid * (2 + lgrid), double);
double *dwork_dMVN, *Sigma, *dens_denom, *w_fycond, *dens_numer;
double *SigmaP, *dens_denomP, *w_fycondP, *dens_numerP, *cSigma;
dwork_dMVN = dwork;
Sigma = dwork + 2; /** space to store Sigma_j (LT(p)) **/
dens_denom = Sigma + LTp; /** space to store denominator when computing conditional densities **/
/** = {sum_{k<K} w_k*f(ycond[i]): i < lcgrid} **/
w_fycond = dens_denom + lcgrid; /** space to store {w_k*f(ycond[i]: i < lcgrid)} for fixed k **/
/** * this is needed when computing conditional cdf's **/
dens_numer = w_fycond + lcgrid; /** space to store numerator when computing conditional densities **/
/*** REMARK: dens_numer will be sorted in this way: ***/
/*** f(y0|ycond=ycond[0]), ..., f(y0|ycond[last]), ..., f(y[p-1]|ycond[0]), ..., f(y[p-1]|ycond[last]) ***/
//.........这里部分代码省略.........
开发者ID:cran，项目名称:mixAK，代码行数:101，代码来源:NMix_PredCondDensCDFMarg.cpp
示例15: rnorm_truncated
void
rnorm_truncated (double *sample, int *n, double *mu,
double *sigma, double *lower, double *upper)
{
int k;
int change;
double a, b;
double logt1 = log(0.150), logt2 = log(2.18), t3 = 0.725, t4 = 0.45;
double z, tmp, lograt;
GetRNGstate();
for (k=0; k<(*n); k++)
{
change = 0;
a = (lower[k] - mu[k])/sigma[k];
b = (upper[k] - mu[k])/sigma[k];
// First scenario
if( (a == R_NegInf) || (b == R_PosInf))
{
if(a == R_NegInf)
{
change = 1;
a = -b;
b = R_PosInf;
}
// The two possibilities for this scenario
if(a <= 0.45) z = norm_rs(a, b);
else z = exp_rs(a, b);
if(change) z = -z;
}
// Second scenario
else if((a * b) <= 0.0)
{
// The two possibilities for this scenario
if((dnorm(a, 0.0, 1.0, 1) <= logt1) || (dnorm(b, 0.0, 1.0, 1) <= logt1))
{
z = norm_rs(a, b);
}
else z = unif_rs(a,b);
}
// Third scenario
else
{
if(b < 0)
{
tmp = b; b = -a; a = -tmp; change = 1;
}
lograt = dnorm(a, 0.0, 1.0, 1) - dnorm(b, 0.0, 1.0, 1);
if(lograt <= logt2) z = unif_rs(a,b);
else if((lograt > logt1) && (a < t3)) z = half_norm_rs(a,b);
else z = exp_rs(a,b);
if(change) z = -z;
}
sample[k] = sigma[k]*z + mu[k];
}
PutRNGstate();
}
开发者ID:YiranChen，项目名称:truncatedNormals，代码行数:63，代码来源:trunc_norm.c
示例16: init_pull_coord
static void init_pull_coord(t_pull_coord *pcrd, int coord_index_for_output,
char *dim_buf,
const char *origin_buf, const char *vec_buf,
warninp_t wi)
{
int m;
dvec origin, vec;
char buf[STRLEN];
if (pcrd->eType == epullCONSTRAINT && (pcrd->eGeom == epullgCYL ||
pcrd->eGeom == epullgDIRRELATIVE ||
pcrd->eGeom == epullgANGLE ||
pcrd->eGeom == epullgANGLEAXIS ||
pcrd->eGeom == epullgDIHEDRAL))
{
gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
epull_names[pcrd->eType],
epullg_names[pcrd->eGeom],
epull_names[epullUMBRELLA]);
}
if (pcrd->eType == epullEXTERNAL)
{
if (pcrd->externalPotentialProvider[0] == '\0')
{
sprintf(buf, "The use of pull type '%s' for pull coordinate %d requires that the name of the module providing the potential external is set with the option %s%d%s",
epull_names[pcrd->eType], coord_index_for_output,
"pull-coord", coord_index_for_output, "-potential-provider");
warning_error(wi, buf);
}
if (pcrd->rate != 0)
{
sprintf(buf, "The use of pull type '%s' for pull coordinate %d requires that the pull rate is zero",
epull_names[pcrd->eType], coord_index_for_output);
warning_error(wi, buf);
}
if (pcrd->eGeom == epullgCYL)
{
/* Warn the user of a PBC restriction, caused by the fact that
* there is no reference value with an external pull potential.
*/
sprintf(buf, "With pull type '%s' and geometry '%s', the distance component along the cylinder axis between atoms in the cylinder group and the COM of the pull group should be smaller than half the box length",
epull_names[pcrd->eType], epullg_names[pcrd->eGeom]);
warning_note(wi, buf);
}
}
process_pull_dim(dim_buf, pcrd->dim, pcrd);
string2dvec(origin_buf, origin);
if (pcrd->group[0] != 0 && dnorm(origin) > 0)
{
gmx_fatal(FARGS, "The pull origin can only be set with an absolute reference");
}
/* Check the given initial reference value and warn for dangerous values */
if (pcrd->eGeom == epullgDIST)
{
if (pcrd->bStart && pcrd->init < 0)
{
sprintf(buf, "The initial reference distance set by pull-coord-init is set to a negative value (%g) with geometry %s while distances need to be non-negative. "
"This may work, since you have set pull-coord-start to 'yes' which modifies this value, but only for certain starting distances. "
"If this is a mistake you may want to use geometry %s instead.",
pcrd->init, EPULLGEOM(pcrd->eGeom), EPULLGEOM(epullgDIR));
warning(wi, buf);
}
}
else if (pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgANGLEAXIS)
{
if (pcrd->bStart && (pcrd->init < 0 || pcrd->init > 180))
{
/* This value of pcrd->init may be ok depending on pcrd->bStart which modifies pcrd->init later on */
sprintf(buf, "The initial reference angle set by pull-coord-init (%g) is outside of the allowed range [0, 180] degrees for geometry (%s). "
"This may work, since you have set pull-coord-start to 'yes' which modifies this value, but only for certain starting angles.",
pcrd->init, EPULLGEOM(pcrd->eGeom));
warning(wi, buf);
}
}
else if (pcrd->eGeom == epullgDIHEDRAL)
{
if (pcrd->bStart && (pcrd->init < -180 || pcrd->init > 180))
{
sprintf(buf, "The initial reference angle set by pull-coord-init (%g) is outside of the allowed range [-180, 180] degrees for geometry (%s). "
"This may work, since you have set pull-coord-start to 'yes' which modifies this value, but only for certain starting angles.",
pcrd->init, EPULLGEOM(pcrd->eGeom));
warning(wi, buf);
}
}
/* Check and set the pull vector */
clear_dvec(vec);
string2dvec(vec_buf, vec);
if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgCYL || pcrd->eGeom == epullgDIRPBC || pcrd->eGeom == epullgANGLEAXIS)
{
if (dnorm2(vec) == 0)
{
gmx_fatal(FARGS, "With pull geometry %s the pull vector can not be 0,0,0",
//.........这里部分代码省略.........
开发者ID:friforever，项目名称:gromacs，代码行数:101，代码来源:readpull.cpp
示例17: do_constraint
//.........这里部分代码省略.........
{
gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
}
{
double q, c_a, c_b, c_c;
c_a = diprod(r_ij[c], r_ij[c]);
c_b = diprod(unc_ij, r_ij[c])*2;
c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
if (c_b < 0)
{
q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
lambda = -q/c_a;
}
else
{
q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
lambda = -c_c/q;
}
if (debug)
{
fprintf(debug,
"Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
c_a, c_b, c_c, lambda);
}
}
/* The position corrections dr due to the constraints */
dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
dr_tot[c] += -lambda*dnorm(r_ij[c]);
break;
case epullgDIR:
case epullgDIRPBC:
case epullgCYL:
/* A 1-dimensional constraint along a vector */
a = 0;
for (m = 0; m < DIM; m++)
{
vec[m] = pcrd->vec[m];
a += unc_ij[m]*vec[m];
}
/* Select only the component along the vector */
dsvmul(a, vec, unc_ij);
lambda = a - ref;
if (debug)
{
fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
}
|
客服电话
APP下载
官方微信
请发表评论