C++ dnorm函数代码示例

//.........这里部分代码省略.........
		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 ;
		}

/* # 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]--;
      }
//.........这里部分代码省略.........

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

}

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

//.........这里部分代码省略.........
    	    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++;

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

//.........这里部分代码省略.........
                
                /* 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];
            

/* 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];
                }
            }
        }
    }
//.........这里部分代码省略.........

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

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

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);
//.........这里部分代码省略.........

double F77_SUB(dnrm)(double *x, double *mu, double *sigma, int *give_log)
{
	return dnorm(*x, *mu, *sigma, *give_log);
}

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

/***** ***************************************************************************************** *****/
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])       ***/

//.........这里部分代码省略.........

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