@@ -48,29 +48,30 @@

   alpha_u = array(0, dim = c(N, K));

   alpha_s = array(0, dim = c(N, K));

-  varp_s1 = array(sqrt(3), dim = c(K, 1));

+  varp_s1 = array(sqrt(5), dim = c(K, 1)); # variance of the proposal distribution 1 for alpha_s

   # sqrt(var)

-  varp_s2 = array(sqrt(10), dim = c(K, 1));

+  varp_s2 = array(sqrt(15), dim = c(K, 1)); #variance of the proposal distribution 2 for alpha_s

   # sqrt(var)

-  varp_s2[K] = sqrt(20);

-  varp_s1[K] = sqrt(10);

+  varp_s2[K] = sqrt(25);

+  varp_s1[K] = sqrt(15);

+  # alpha_s is proposed using a mixture distribution

-  pvar_s = array(0.8, dim = c(K, 1));

+  pvar_s = array(0.8, dim = c(K, 1)); # mixing for the mixture proposal distribution

   pvar_s[K] = 0.6;

-  varp_u = array(sqrt(8), dim = c(K, 1));

+  varp_u = array(sqrt(10), dim = c(K, 1)); #variance of the proposal distribution for alpha_u

   pp = array(0.65, dim = c(I, 1))

   pb1 <- clamp(1.5 / median(indi[, K]), 0, 0.9)

   pb2 <- clamp(5.0 / median(indi[, K]), 0, 0.9)

-  lambda_s = rep(0, K);

+  lambda_s = rep(0, K);  

   lambda_s[1:K1] = (10 ^ -2) * max(N_s, N_u)

   lambda_s[K] = max(N_s, N_u) - sum(lambda_s[1:K1])

   lambda_u = lambda_s

-  alpha_u[1, 1:(K - 1)] = 10

+  alpha_u[1, 1:(K - 1)] = 10 #initializaion 

   alpha_u[1, K] = 150

-  alpha_s[1, 1:(K - 1)] = 10

+  alpha_s[1, 1:(K - 1)] = 10 #initialization 

   alpha_s[1, K] = 100

   #################### acceptance rate ###########################

@@ -84,7 +85,8 @@

     if (tt %% 1000 == 0) vmessage("Iteration ", tt, " of ", N, ".")

     # update alphau

-    res2 <- .Call(C_updatealphau_noPu_Exp, alphaut = alpha_u[tt - 1,], n_s = n_s, n_u = n_u, I = I, K = K, lambda_u = lambda_u, var_p = varp_u, ttt = ttt, gammat = gamma[,, tt - 1])

+    #res2 <- .Call(C_updatealphau_noPu_Exp, alphaut = alpha_u[tt - 1,], n_s = n_s, n_u = n_u, I = I, K = K, lambda_u = lambda_u, var_p = varp_u, ttt = ttt, gammat = gamma[,, tt - 1])

+    res2 <- .Call(C_updatealphau_noPu_Exp_MH, alphaut = alpha_u[tt - 1,], n_s = n_s, n_u = n_u, I = I, K = K, lambda_u = lambda_u, var_p = varp_u, gammat = gamma[,, tt - 1])

     if (length(alpha_u[tt, ]) != length(res2$alphau_tt)) {

       vmessage("res2 alphau_tt length:", length(res2$alphau_tt), "\n")

       vmessage("alpha_u[tt,] length: ", length(alpha_u[tt,]), "\n")

@@ -95,6 +97,7 @@

     #update gamma

     res1 <- .Call(C_updategammak_noPu, n_s = n_s, n_u = n_u, gammat = gamma[,, tt - 1], I = I, K = K, SS = SS, alphau = alpha_u[tt,], alphas = alpha_s[tt - 1,], alpha = 1, mk = mk, Istar = Istar,

       mKstar = mKstar, pp = pp, pb1 = pb1, pb2 = pb2, indi = indi)

+    

     gamma[,, tt] = res1$gamma_tt;

     if (length(A_gm[, tt]) != length(res1$Ag)) {

       vmessage("res1 Ag length: ", length(res1$Ag), "\n")

@@ -107,7 +110,8 @@

     mKstar = res1$mKstar;

     # update alphas

-    res3 <- .Call(C_updatealphas_Exp, alphast = alpha_s[tt - 1,], n_s = n_s, K = K, I = I, lambda_s = lambda_s, gammat = gamma[,, tt], var_1 = varp_s1, var_2 = varp_s2, p_var = pvar_s, ttt = ttt)

+    # res3 <- .Call(C_updatealphas_Exp, alphast = alpha_s[tt - 1,], n_s = n_s, K = K, I = I, lambda_s = lambda_s, gammat = gamma[,, tt], var_1 = varp_s1, var_2 = varp_s2, p_var = pvar_s, ttt = ttt)

+    res3 <- .Call(C_updatealphas_Exp_MH, alphast = alpha_s[tt - 1,], n_s = n_s, K = K, I = I, lambda_s = lambda_s, gammat = gamma[,, tt], var_1 = varp_s1, var_2 = varp_s2, p_var = pvar_s)

     if (length(alpha_s[tt, ]) != length(res3$alphas_tt)) {

       vmessage("res3 alphas_tt length:", length(res3$alphas_tt), "\n")

       vmessage("alpha_s[tt,] length: ", length(alpha_s[tt,]), "\n")

@@ -1,4 +1,3 @@

-// Copyright [2014] <Fred Hutchinson Cancer Research Center>

 // This file was automatically generated by Kmisc::registerFunctions()

 #include <R.h>

@@ -34,6 +33,15 @@ SEXP updatealphas_Exp(SEXP alphast,

                       SEXP var_2,

                       SEXP p_var,

                       SEXP ttt);

+SEXP updatealphas_Exp_MH(SEXP alphast,

+                         SEXP n_s,

+                         SEXP K,

+                         SEXP I,

+                         SEXP lambda_s,

+                         SEXP gammat,

+                         SEXP var_1,

+                         SEXP var_2,

+                         SEXP p_var);

 SEXP updatealphau_noPu_Exp(SEXP alphaut,

                            SEXP n_s,

                            SEXP n_u,

@@ -43,6 +51,14 @@ SEXP updatealphau_noPu_Exp(SEXP alphaut,

                            SEXP var_p,

                            SEXP ttt,

                            SEXP gammat);

+SEXP updatealphau_noPu_Exp_MH(SEXP alphaut,

+                              SEXP n_s,

+                              SEXP n_u,

+                              SEXP I,

+                              SEXP K,

+                              SEXP lambda_u,

+                              SEXP var_p,

+                              SEXP gammat);

 SEXP updategammak_noPu(SEXP n_s,

                        SEXP n_u,

                        SEXP gammat,

@@ -69,7 +85,9 @@ R_CallMethodDef callMethods[] = {

     {"C_COMPASS_CellCounts", (DL_FUNC) & _COMPASS_CellCounts, 2},

     {"C_samplePuPs", (DL_FUNC) & samplePuPs, 9},

     {"C_updatealphas_Exp", (DL_FUNC) & updatealphas_Exp, 10},

+    {"C_updatealphas_Exp_MH", (DL_FUNC) & updatealphas_Exp_MH, 9},

     {"C_updatealphau_noPu_Exp", (DL_FUNC) & updatealphau_noPu_Exp, 9},

+    {"C_updatealphau_noPu_Exp_MH", (DL_FUNC) & updatealphau_noPu_Exp_MH, 8},

     {"C_updategammak_noPu", (DL_FUNC) & updategammak_noPu, 16},

     {NULL, NULL, 0}};

@@ -1,5 +1,3 @@

-// Copyright [2014] <Fred Hutchinson Cancer Research Center>

-

 // Generated by using Rcpp::compileAttributes() -> do not edit by hand

 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

@@ -10,25 +8,24 @@ using namespace Rcpp;

 // CellCounts

 IntegerMatrix CellCounts(List x, List combos);

 RcppExport SEXP _COMPASS_CellCounts(SEXP xSEXP, SEXP combosSEXP) {

-  BEGIN_RCPP

-  Rcpp::RObject rcpp_result_gen;

-  Rcpp::RNGScope rcpp_rngScope_gen;

-  Rcpp::traits::input_parameter<List>::type x(xSEXP);

-  Rcpp::traits::input_parameter<List>::type combos(combosSEXP);

-  rcpp_result_gen = Rcpp::wrap(CellCounts(x, combos));

-  return rcpp_result_gen;

-  END_RCPP

+BEGIN_RCPP

+    Rcpp::RObject rcpp_result_gen;

+    Rcpp::RNGScope rcpp_rngScope_gen;

+    Rcpp::traits::input_parameter< List >::type x(xSEXP);

+    Rcpp::traits::input_parameter< List >::type combos(combosSEXP);

+    rcpp_result_gen = Rcpp::wrap(CellCounts(x, combos));

+    return rcpp_result_gen;

+END_RCPP

 }

 // CellCounts_character

 IntegerMatrix CellCounts_character(List data, List combinations);

-RcppExport SEXP _COMPASS_CellCounts_character(SEXP dataSEXP,

-                                              SEXP combinationsSEXP) {

-  BEGIN_RCPP

-  Rcpp::RObject rcpp_result_gen;

-  Rcpp::RNGScope rcpp_rngScope_gen;

-  Rcpp::traits::input_parameter<List>::type data(dataSEXP);

-  Rcpp::traits::input_parameter<List>::type combinations(combinationsSEXP);

-  rcpp_result_gen = Rcpp::wrap(CellCounts_character(data, combinations));

-  return rcpp_result_gen;

-  END_RCPP

+RcppExport SEXP _COMPASS_CellCounts_character(SEXP dataSEXP, SEXP combinationsSEXP) {

+BEGIN_RCPP

+    Rcpp::RObject rcpp_result_gen;

+    Rcpp::RNGScope rcpp_rngScope_gen;

+    Rcpp::traits::input_parameter< List >::type data(dataSEXP);

+    Rcpp::traits::input_parameter< List >::type combinations(combinationsSEXP);

+    rcpp_result_gen = Rcpp::wrap(CellCounts_character(data, combinations));

+    return rcpp_result_gen;

+END_RCPP

 }

new file mode 100644

@@ -0,0 +1,183 @@

+#include <Rcpp.h>

+

+extern "C" double digamma(double);

+

+// [[register]]

+RcppExport SEXP updatealphas_Exp_MH(SEXP alphast, SEXP n_s, SEXP K, SEXP I,

+                                    SEXP lambda_s, SEXP gammat, SEXP var_1,

+                                    SEXP var_2, SEXP p_var) {

+  BEGIN_RCPP

+  Rcpp::NumericVector xalphast(alphast);

+

+  Rcpp::IntegerMatrix xn_s(n_s);

+  Rcpp::IntegerMatrix xgammat(gammat);

+  int xI = Rcpp::as<int>(I);

+  int xK = Rcpp::as<int>(K);

+  Rcpp::NumericVector sqrt_var1(var_1);

+  Rcpp::NumericVector sqrt_var2(var_2);

+  // int xtt = Rcpp::as<int>(ttt);

+  Rcpp::NumericVector xlambda_s(lambda_s);

+  Rcpp::IntegerVector xAalphas(xK);

+

+  Rcpp::RNGScope scope;

+  Rcpp::NumericVector xp_var(p_var);  // proposal mixture

+

+  // double delF = 0.0;

+  double psik = 0.;

+  double log1 = 0.0;

+  double log2 = 0.0;

+  double sums = 0.;

+  double sum_alp_ns = 0.0;

+  double sum_alp = 0.0;

+  double sum_gl_alp = 0.0;

+  double sum_gl_alp_ns = 0.0;

+  int flag1 = 0;

+  int flagkk = 0;

+  int lp1 = 0;

+  for (int kk = 0; kk < xK; kk++) {

+    // delF = 0.0;

+    psik = digamma(xalphast[kk]);

+    log1 = 0.0;

+    log2 = 0.0;

+    for (int i = 0; i < xI; i++) {

+      lp1 = 0;

+      for (int k = 0; k < xK; k++) {

+        if (xgammat(i, k) == 1) {

+          lp1 += 1;

+        }

+      }

+

+      std::vector<int> p1(lp1);

+      flag1 = 0;

+      flagkk = 0;

+      for (int k = 0; k < xK; k++) {

+        if (xgammat(i, k) == 1) {

+          p1[flag1] = k;

+          flag1 += 1;

+          if (k == kk) {

+            flagkk = 1;

+          }

+        }

+      }

+      sum_alp_ns = 0.0;

+      sum_alp = 0.0;

+      sum_gl_alp = 0.0;

+      sum_gl_alp_ns = 0.0;

+      for (int k = 0; k < lp1; k++) {

+        sums = xalphast[p1[k]] + xn_s(i, p1[k]);

+        sum_alp_ns += sums;

+        sum_alp += xalphast[p1[k]];

+        sum_gl_alp += lgamma(xalphast[p1[k]]);

+        sum_gl_alp_ns += lgamma(sums);

+      }

+      // if (flagkk > 0) {

+      // delF += digamma(xn_s(i, kk) + xalphast[kk]) - psik -

+      //  digamma(sum_alp_ns) + digamma(sum_alp);

+      //}

+      if (lp1 > 0) {

+        log2 += -(sum_gl_alp - lgamma(sum_alp)) +

+                (sum_gl_alp_ns - lgamma(sum_alp_ns));

+      }

+    }

+    // double mean_p = std::max(0.01, xalphast[kk] + delF / xtt);

+    Rcpp::NumericVector alpha_s_p = Rcpp::rnorm(1, xalphast[kk], sqrt_var1[kk]);

+

+    if (Rcpp::as<double>(Rcpp::rbinom(1, 1, xp_var[kk])) == 1) {

+      alpha_s_p = Rcpp::rnorm(1, xalphast[kk], sqrt_var1[kk]);

+    } else {

+      alpha_s_p = Rcpp::rnorm(1, xalphast[kk], sqrt_var2[kk]);

+    }

+

+    if (alpha_s_p[0] > 0.0) {

+      std::vector<double> alp(xK);

+

+      for (int i = 0; i < xK; i++) {

+        alp[i] = xalphast[i];

+      }

+      alp[kk] = alpha_s_p[0];

+      // log2 += log(xp_var[kk]*gsl_ran_gaussian_pdf(alp[kk]-mean_p,

+      // sqrt_var1[kk])+(1-xp_var[kk])*gsl_ran_gaussian_pdf(alp[kk]-mean_p,

+      // sqrt_var2[kk]));

+      log2 +=

+          log(xp_var[kk] * Rf_dnorm4(alp[kk], xalphast[kk], sqrt_var1[kk], 0) +

+              (1 - xp_var[kk]) *

+                  Rf_dnorm4(alp[kk], xalphast[kk], sqrt_var2[kk], 0));

+      // delF = 0.0;

+      psik = digamma(alp[kk]);

+      for (int i = 0; i < xI; i++) {

+        lp1 = 0;

+        for (int k = 0; k < xK; k++) {

+          if (xgammat(i, k) == 1) {

+            lp1 += 1;

+          }

+        }

+

+        std::vector<int> p1(lp1);

+        flag1 = 0;

+        flagkk = 0;

+        for (int k = 0; k < xK; k++) {

+          if (xgammat(i, k) == 1) {

+            p1[flag1] = k;

+            flag1 += 1;

+            if (k == kk) {

+              flagkk = 1;

+            }

+          }

+        }

+

+        sum_alp_ns = 0.0;

+        sum_alp = 0.0;

+        sum_gl_alp = 0.0;

+        sum_gl_alp_ns = 0.0;

+        for (int k = 0; k < lp1; k++) {

+          sums = alp[p1[k]] + xn_s(i, p1[k]);

+          sum_alp_ns += sums;

+          sum_alp += alp[p1[k]];

+          sum_gl_alp += lgamma(alp[p1[k]]);

+          sum_gl_alp_ns += lgamma(sums);

+        }

+        // if (flagkk > 0) {

+        //   delF += digamma(xn_s(i, kk) + xalphast[kk]) - psik -

+        //     digamma(sum_alp_ns) + digamma(sum_alp);

+        // }

+        if (lp1 > 0) {

+          log1 += -(sum_gl_alp - lgamma(sum_alp)) +

+                  (sum_gl_alp_ns - lgamma(sum_alp_ns));

+        }

+      }

+      // mean_p = std::max(0.01, alp[kk] + delF / xtt);

+      // log1 +=log(xp_var[kk]*gsl_ran_gaussian_pdf(xalphast[kk]-mean_p,

+      // sqrt_var1[kk])+(1-xp_var[kk])*gsl_ran_gaussian_pdf(xalphast[kk]-mean_p,

+      // sqrt_var2[kk]));

+      log1 +=

+          log(xp_var[kk] * Rf_dnorm4(xalphast[kk], alp[kk], sqrt_var1[kk], 0) +

+              (1 - xp_var[kk]) *

+                  Rf_dnorm4(xalphast[kk], alp[kk], sqrt_var2[kk], 0));

+

+      // log1 += log(gsl_ran_exponential_pdf(alp[kk],xlambda_s[kk]));

+      // //exponential prior

+      log1 += Rf_dexp(alp[kk], xlambda_s[kk], 1);

+

+      // log2 +=

+      // log(gsl_ran_exponential_pdf(xalphast[kk],xlambda_s[kk]));//exponential

+      // prior

+      log2 += Rf_dexp(xalphast[kk], xlambda_s[kk], 1);

+

+      // if (alp[kk]<0 || alp[kk]>xlambda_s[kk]) {log1+=log(0);} //Uniform prior

+      // if (xalphast[kk]<0 || xalphast[kk]>xlambda_s[kk]) {log2+=log(0);}

+      // //Uniform prior

+

+      if (log(Rcpp::as<double>(Rcpp::runif(1))) <= (log1 - log2)) {

+        xalphast[kk] = alp[kk];

+        xAalphas[kk] = 1;

+      } else {

+        xAalphas[kk] = 0;

+      }

+    }

+  }

+

+  return Rcpp::List::create(Rcpp::Named("alphas_tt") = xalphast,

+                            Rcpp::Named("Aalphas") = xAalphas);

+

+  END_RCPP

+}

new file mode 100644

@@ -0,0 +1,209 @@

+#include <Rcpp.h>

+extern "C" double digamma(double);

+

+// Standard MH algorithm for updating alpha_u

+

+// [[register]]

+RcppExport SEXP updatealphau_noPu_Exp_MH(SEXP alphaut, SEXP n_s, SEXP n_u,

+                                         SEXP I, SEXP K, SEXP lambda_u,

+                                         SEXP var_p, SEXP gammat) {

+  BEGIN_RCPP

+  Rcpp::IntegerMatrix xgammat(gammat);

+  Rcpp::NumericVector xalphaut(alphaut);

+  Rcpp::IntegerMatrix xn_s(n_s);

+  Rcpp::IntegerMatrix xn_u(n_u);

+  int xI = Rcpp::as<int>(I);

+  int xK = Rcpp::as<int>(K);

+  Rcpp::NumericVector sqrt_var(var_p);

+  // int xtt = Rcpp::as<int>(ttt);

+  Rcpp::NumericVector xlambda_u(lambda_u);

+  Rcpp::IntegerVector xAalphau(xK);

+  Rcpp::RNGScope scope;

+

+  // double delF = 0.0;

+  double log1 = 0.0;

+  double log2 = 0.0;

+  double sum_alphau = 0.0;

+  int flag1 = 0;

+  int flag0 = 0;

+  int flagkk = 0;

+  int lp0 = 0;

+  int lp1 = 0;

+  double sum_nusalphau = 0.0;

+  double sum_nualphau = 0.0;

+  double sums = 0.;

+  for (int kk = 0; kk < xK; kk++) {

+    // delF = 0.0;

+    log1 = 0.0;

+    log2 = 0.0;

+    sum_alphau = 0.0;

+    for (int s = 0; s < xK; s++) {

+      sum_alphau += xalphaut[s];

+    }

+    log2 -= xI * lgamma(xalphaut[kk]);

+    // delF += xI * (digamma(sum_alphau) - digamma(xalphaut[kk]));

+    log2 += xI * lgamma(sum_alphau);

+    for (int i = 0; i < xI; i++) {

+      lp1 = 0;

+      for (int k = 0; k < xK; k++) {

+        if (xgammat(i, k) == 1) {

+          lp1 += 1;

+        }

+      }

+      lp0 = xK - lp1;

+      std::vector<int> p1(lp1);

+      flag1 = 0;

+      std::vector<int> p0(lp0);

+      flag0 = 0;

+      flagkk = 0;  // whether gamma_k = 1

+

+      for (int k = 0; k < xK; k++) {

+        if (xgammat(i, k) == 1) {

+          p1[flag1] = k;

+          flag1 += 1;

+          if (k == kk) {

+            flagkk = 1;

+          }

+        } else {

+          p0[flag0] = k;

+          flag0 += 1;

+        }

+      }

+      if (flagkk == 1) {

+        log2 += lgamma(xn_u(i, kk) + xalphaut[kk]);

+        // delF += digamma(xn_u(i, kk) + xalphaut[kk]);

+        sum_nualphau = 0.0;

+        sum_nusalphau = 0.0;

+        for (int k = 0; k < lp1; k++) {

+          sums = xn_u(i, p1[k]) + xalphaut[p1[k]];

+          sum_nualphau += sums;

+          sum_nusalphau += (sums + xn_s(i, p1[k]));

+        }

+        log2 -= lgamma(sum_nualphau);

+        log2 += lgamma(sum_nusalphau + 1);

+        // delF -= digamma(sum_nualphau);

+        // delF += digamma(sum_nusalphau + 1);

+

+        for (int k = 0; k < lp0; k++) {

+          sum_nusalphau += (xn_u(i, p0[k]) + xalphaut[p0[k]] + xn_s(i, p0[k]));

+        }

+        // delF -= digamma(sum_nusalphau + 1);

+        log2 -= lgamma(sum_nusalphau + 1);

+      } else {

+        log2 += lgamma(xn_u(i, kk) + xalphaut[kk] + xn_s(i, kk));

+        // delF += digamma(xn_u(i, kk) + xalphaut[kk] + xn_s(i, kk));

+        sum_nusalphau = 0.0;

+        for (int k = 0; k < xK; k++) {

+          sum_nusalphau += xn_u(i, k) + xalphaut[k] + xn_s(i, k);

+        }

+        log2 -= lgamma(sum_nusalphau + 1);

+        // delF -= digamma(sum_nusalphau + 1);

+      }

+    }

+    // double mean_p = std::max(0.01, xalphaut[kk] + delF / xtt);

+    Rcpp::NumericVector alpha_u_p = Rcpp::rnorm(1, xalphaut[kk], sqrt_var[kk]);

+    if (alpha_u_p[0] > 0.0) {

+      std::vector<double> alp(xK);

+      for (int i = 0; i < xK; i++) {

+        alp[i] = xalphaut[i];

+      }

+      alp[kk] = alpha_u_p[0];

+

+      // log2 += log(gsl_ran_gaussian_pdf(alp[kk]-mean_p, sqrt_var[kk]));

+      log2 += Rf_dnorm4(alp[kk], xalphaut[kk], sqrt_var[kk], 1);

+

+      // delF = 0.0;

+      sum_alphau = 0.0;

+      for (int s = 0; s < xK; s++) {

+        sum_alphau += alp[s];

+      }

+      log1 -= xI * lgamma(alp[kk]);

+      // delF += xI * (digamma(sum_alphau) - digamma(alp[kk]));

+      log1 += xI * lgamma(sum_alphau);

+      for (int i = 0; i < xI; i++) {

+        lp1 = 0;

+        for (int k = 0; k < xK; k++) {

+          if (xgammat(i, k) == 1) {

+            lp1 += 1;

+          }

+        }

+        lp0 = xK - lp1;

+        std::vector<int> p1(lp1);

+        flag1 = 0;

+        std::vector<int> p0(lp0);

+        flag0 = 0;

+        flagkk = 0;  // whether gamma_k = 1

+

+        for (int k = 0; k < xK; k++) {

+          if (xgammat(i, k) == 1) {

+            p1[flag1] = k;

+            flag1 += 1;

+            if (k == kk) {

+              flagkk = 1;

+            }

+          } else {

+            p0[flag0] = k;

+            flag0 += 1;

+          }

+        }

+        if (flagkk == 1) {

+          log1 += lgamma(xn_u(i, kk) + alp[kk]);

+          // delF += digamma(xn_u(i, kk) + alp[kk]);

+          sum_nualphau = 0.0;

+          sum_nusalphau = 0.0;

+          for (int k = 0; k < lp1; k++) {

+            sums = xn_u(i, p1[k]) + alp[p1[k]];

+            sum_nualphau += sums;

+            sum_nusalphau += (sums + xn_s(i, p1[k]));

+          }

+          log1 -= lgamma(sum_nualphau);

+          log1 += lgamma(sum_nusalphau + 1);

+          // delF -= digamma(sum_nualphau);

+          // delF += digamma(sum_nusalphau + 1);

+

+          for (int k = 0; k < lp0; k++) {

+            sum_nusalphau += (xn_u(i, p0[k]) + alp[p0[k]] + xn_s(i, p0[k]));

+          }

+          // delF -= digamma(sum_nusalphau + 1);

+          log1 -= lgamma(sum_nusalphau + 1);

+        } else {

+          log1 += lgamma(xn_u(i, kk) + alp[kk] + xn_s(i, kk));

+          // delF += digamma(xn_u(i, kk) + alp[kk] + xn_s(i, kk));

+          sum_nusalphau = 0.0;

+          for (int k = 0; k < xK; k++) {

+            sum_nusalphau += xn_u(i, k) + alp[k] + xn_s(i, k);

+          }

+          log1 -= lgamma(sum_nusalphau + 1);

+          // delF -= digamma(sum_nusalphau + 1);

+        }

+      }

+      // mean_p = std::max(0.01, alp[kk] + delF / xtt);

+

+      // log1 +=log(gsl_ran_gaussian_pdf(xalphaut[kk]-mean_p, sqrt_var[kk]));

+      log1 += Rf_dnorm4(xalphaut[kk], alp[kk], sqrt_var[kk], 1);

+

+      // log1 += log(gsl_ran_exponential_pdf(alp[kk],xlambda_u[kk]));

+      // //exponential prior

+      log1 += Rf_dexp(alp[kk], xlambda_u[kk], 1);

+

+      // log2 += log(gsl_ran_exponential_pdf(xalphaut[kk],xlambda_u[kk]));

+      // //exponential prior

+      log2 += Rf_dexp(xalphaut[kk], xlambda_u[kk], 1);

+

+      // if (alp[kk]<0 || alp[kk]>xlambda_u[kk]) {log1+=log(0);} //Uniform prior

+      // if (xalphaut[kk]<0 || xalphaut[kk]>xlambda_u[kk]) {log2+=log(0);}

+      // //Uniform prior

+

+      if (log(Rcpp::as<double>(Rcpp::runif(1))) <= (log1 - log2)) {

+        xalphaut[kk] = alp[kk];

+        xAalphau[kk] = 1;

+      } else {

+        xAalphau[kk] = 0;

+      }

+    }

+  }

+  return Rcpp::List::create(Rcpp::Named("alphau_tt") = xalphaut,

+                            Rcpp::Named("Aalphau") = xAalphau);

+

+  END_RCPP

+}