SciCodePile/SciCode-Domain-Code

1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Genesis","yandorazhang/GENESIS","src/EM_EM.cpp",".cpp","44549","1245","//------------------------------------------------3//  EM_EM.cpp4//  Goal: Estimate effect size distribution.5//6//  Author: Yan (Dora) Zhang7//  Email: yandorazhang@gmail.com8//------------------------------------------------9 10#include <cstdlib>11#include <time.h>12#include <iostream>     // std::cout13#include <algorithm>    // std::min14#include <RcppArmadillo.h>15//[[Rcpp::depends(RcppArmadillo)]]16 17// --------------------------------//--------------------------------18#include <omp.h>19//[[Rcpp::plugins(openmp)]]20// --------------------------------//--------------------------------21 22using namespace arma;23using namespace Rcpp;24using namespace std;25using namespace stats;26 27 28//--------------------------------29//--------------------------------30// [[Rcpp::export]]31vec modification_loc(vec inx_name, int K, int mx_k){32  vec inx_loc(mx_k);33  inx_loc.fill(0);34  35  for(int i = 0; i < K; i++){36    inx_loc(inx_name(i) - 1) = i+1;37  }38  return(inx_loc);39}40 41 42 43//--------------------------------44//--------------------------------45// [[Rcpp::export]]46long double cpnorm(long double x) // R: pnorm()47{48  // constants49  long double a1 =  0.254829592;50  long double a2 = -0.284496736;51  long double a3 =  1.421413741;52  long double a4 = -1.453152027;53  long double a5 =  1.061405429;54  long double p  =  0.3275911;55  56  // Save the sign of x57  int sign = 1;58  if (x < 0)59    sign = -1;60  x = fabs(x)/sqrt(2.0);61  62  // A&S formula 7.1.2663  long double t = 1.0/(1.0 + p*x);64  long double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);65  66  return 0.5*(1.0 + sign*y);67}68 69//--------------------------------70//--------------------------------71// [[Rcpp::export]]72long double sumfactorial(int n) // log(factorial(n)) = log(n!)73{74  if(n > 1)75    return log(n) + sumfactorial(n - 1);76  else77    return 0;78}79 80//--------------------------------81//--------------------------------82// [[Rcpp::export]]83long double sumfactorial_rev(int n, const int &k0) // log(n!/k0!)84{85  if(n > k0)86    return log(n) + sumfactorial_rev(n-1,k0);87  else88    return 0;89}90 91//--------------------------------92//--------------------------------93// [[Rcpp::export]]94long double cdnorm(long double x, long double mean, long double sd, bool loglog) // R: dnorm()95{96  long double res;97  res = -0.5*log((2.0*PI))  - log(sd)- pow(x-mean, 2)/(2.0*pow(sd,2));98  if(loglog){return res;}99  else{return exp(res);}100}101 102//--------------------------------103//--------------------------------104// [[Rcpp::export]]105long double cdbinom(const int &k, const int &size, long double prob, bool loglog) // R: dbinom() 106{107  long double res;108  res = sumfactorial_rev(size,size-k) - sumfactorial(k) + k*log(prob) + (size-k)*log(1-prob);109  if(loglog){return res;}110  else{return exp(res);}111}112 113//--------------------------------114//--------------------------------115// [[Rcpp::export]]116long double cdmultinom3(const int k0, const int k1, const int k2,  vec prob, bool loglog) // R: dbmultinom()117{118    long double res;119    120    res = sumfactorial_rev(k0+k1+k2,k2) - sumfactorial(k0) - sumfactorial(k1)121    + k0*log(prob(0)) + k1*log(prob(1)) + k2*log(prob(2));122    123    if(loglog){return res;}124    else{return exp(res);}125}126 127 128//--------------------------------//--------------------------------129// 2-component model130//--------------------------------//--------------------------------131 132//--------------------------------133//--------------------------------134// [[Rcpp::export]]135long double loglikelihood(const vec & par, const vec &betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int & num_threads) // loglikelihood function136{137  int K = betahat.n_elem;138  long double marginal_likelihood;139  long double pi1 = par(0);140  long double sigsq = par(1);141  long double a = par(2);142  long double res = 0;143  long double y=0;144  long double tem;145  int k,j;146  // --------------------------------//--------------------------------147  omp_set_num_threads(num_threads);148  #pragma omp parallel for shared(betahat,pi1, sigsq, a, varbetahat,ldscore,c0,Nstar) private(marginal_likelihood,y,k,j,tem) reduction(+:res)149  // --------------------------------//--------------------------------150  for(k=0; k<K; k++){151    marginal_likelihood = 0;152    y = 0;153    for(j=0; j<std::min(int(Nstar(k))+1, c0+1); j++){154      tem = cdbinom(j,Nstar(k),pi1,false);155      y += tem;156      marginal_likelihood += tem*cdnorm(betahat(k),0, sqrt( (j*sigsq)*ldscore(k)/Nstar(k) + varbetahat(k) + a), false);157    }158    res += log(marginal_likelihood/y);159  }160  161  return res;162}163 164//--------------------------------165// weights function166//--------------------------------167// [[Rcpp::export]]168mat weight(const vec & par, const vec &betahat, const vec &varbetahat, const vec &ldscore, const int & c0, const vec & Nstar, const int &num_threads) // E-step: weights function169{170  int K = betahat.n_elem;171  mat w(K,c0+1);172  w.fill(0.0);173  long double pi1 = par(0);174  long double sigsq = par(1);175  long double a  = par(2);176  vec te;177  long double tem;178  int k,j;179  180  // weights formula.181  // --------------------------------//--------------------------------182  omp_set_num_threads(num_threads);183  #pragma omp parallel for shared(sigsq,a, pi1,betahat,varbetahat,ldscore,Nstar,w,c0) private(k,j,tem)184  // --------------------------------//--------------------------------185  for( k=0; k<K; k++){186    for(j=0; j<std::min(int(Nstar(k))+1, c0+1); j++){187      tem = cdbinom(j,Nstar(k),pi1,false);188      w(k,j) = tem * cdnorm(betahat(k),0, sqrt( (j*sigsq)*ldscore(k)/Nstar(k) + varbetahat(k) + a), false);189    }190  }191  te = sum(w,1);192  te = pow(te, -1.0);193  w = (diagmat(te)) * w;194  195  return w;196}197 198 199//--------------------------------200// weights and loglikelihood function201//--------------------------------202// [[Rcpp::export]]203List weight_loglikelihood(const vec & par, const vec &betahat, const vec &varbetahat, const vec &ldscore, const int & c0, const vec & Nstar, const int &num_threads) // E-step: weight-loglikelihood function204{205  int K = betahat.n_elem;206  mat w(K,c0+1);207  w.fill(0.0);208  long double pi1 = par(0);209  long double sigsq = par(1);210  long double a  = par(2);211  vec te;212  long double tem, tem1;213  int k,j;214  long double marginal_likelihood;215  long double res = 0;216  long double y=0;217 218  // --------------------------------//--------------------------------219  omp_set_num_threads(num_threads);220  #pragma omp parallel for shared(sigsq,a, pi1,betahat,varbetahat,ldscore,Nstar,w,c0) private(marginal_likelihood,y,k,j,tem,tem1) reduction(+:res)221  // --------------------------------//--------------------------------222  for(k=0; k<K; k++){223    marginal_likelihood = 0;224    y = 0;225    for(j=0; j<std::min(int(Nstar(k))+1, c0+1); j++){226      tem = cdbinom(j,Nstar(k),pi1,false);227      tem1 = tem * cdnorm(betahat(k),0, sqrt( (j*sigsq)*ldscore(k)/Nstar(k) + varbetahat(k) + a), false);228      w(k,j) = tem1;229      y += tem;230      marginal_likelihood += tem1;231    }232    res += log(marginal_likelihood/y);233  }234  te = sum(w,1);235  te = pow(te, -1.0);236  w = (diagmat(te)) * w;237 238  return List::create(239    _[""w""] = w,240    _[""llk""] = res241  );242}243 244 245 246//--------------------------------247//--------------------------------248// [[Rcpp::export]]249long double update_pi1(const mat & w, const vec & Nstar, const int &num_threads)250{251  long double res1=0;252  long double res2=0;253  long double res;254  int c0 = w.n_cols-1;255  int K = Nstar.n_elem;256  int k,j;257  // --------------------------------//--------------------------------258  omp_set_num_threads(num_threads);259  #pragma omp parallel for shared(w,Nstar,c0) private(k,j) reduction(+:res1,res2)260  // --------------------------------//--------------------------------261  for(k=0; k<K; k++){262    for(j=0; j<std::min(int(Nstar(k))+1, c0+1); j++){263      res1 += w(k,j)*j;264      res2 += w(k,j)*Nstar(k);265    }266  }267  268  res = res1/res2;269  270  return res;271}272 273//--------------------------------274//--------------------------------275// [[Rcpp::export]]276vec onestep_varcomponent(const vec varcomponent,  const mat & w, const vec &betahat, const vec &varbetahat, const vec &ldscore, const vec & Nstar, const int &num_threads)277{278  279  int K = w.n_rows;280  int c0 = w.n_cols-1;281  long double sigsq = varcomponent(0);282  long double a = varcomponent(1);283  long double tem;284  long double s1=0, dd1=0, s2=0, dd2=0, dd12=0, det=0;285  vec result(2);286  287  int k,j;288  // --------------------------------//--------------------------------289  omp_set_num_threads(num_threads);290  #pragma omp parallel for shared(w, sigsq,a,varbetahat,ldscore,Nstar,c0) private(k,tem,j) reduction(+:s1,dd1,s2,dd2,dd12)291  // --------------------------------//--------------------------------292  for(k=0; k<K; k++)293  {294    for(j=0; j<std::min(int(Nstar(k))+1, c0+1); j++)295    {296      tem = (j*sigsq) * ldscore(k) /double(Nstar(k)) + varbetahat(k) + a;297      s1 += w(k,j)*j*ldscore(k)/(2.0*Nstar(k))*(pow(betahat(k),2)/pow(tem,2) - 1.0/(tem)  );298      s2 += w(k,j)/(2.0)*(pow(betahat(k),2)/pow(tem,2) - 1.0/(tem));299      300      dd1 += w(k, j)*j*ldscore(k)*j*ldscore(k)/(2.0*Nstar(k)*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));301      dd2 += w(k,j)/(2.0)*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));302      dd12 += w(k,j)*j*ldscore(k)/(2.0*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));303    }304  }305  306  det = (dd1*dd2 - pow(dd12,2));307  308  result(0) = sigsq - (dd2*s1 - dd12*s2)/det;309  result(1) = a - (dd1*s2 - dd12*s1)/det;310  311  return result;312}313 314//--------------------------------315//--------------------------------316// [[Rcpp::export]]317vec EM_func(const vec &par_start,318            const vec &betahat, const vec & varbetahat, const vec & ldscore,  const vec & Nstar, const int & M,319            int c0, const long double &eps1,const long double &eps2,const long double &eps3,const long double &eps, const int &Meps, const int &steps, const int &num_threads, const bool &print,320            const int &printfreq, const bool &stratification)321{322  vec par(3), prev_par(3);323  par(0) = par_start(0);324  par(1) = par_start(1);325  par(2) = par_start(2);326  long double llk = 0, error_pi, error_sigsq, error_a, increase_ll;327  long double prev_llk=0, pi1;328  vec old_sig(2);329  vec new_sig(2);330  mat w;331  vec result(10);332  List wllk;333  clock_t start_w, finish_w, finish_pi1, finish_sigsq;334 335  Rcout << ""Iteration, prev_loglikelihood, pic, sigmasq, a, Heritability, c0, Seconds(weight_llk), Seconds(pi), Seconds(variance_components)"";336  Rcout << endl;337 338  for(int i=0; i<Meps; i++){339    prev_llk = llk;340    prev_par(0) = par(0); prev_par(1) = par(1); prev_par(2) = par(2);341 342    // E-step: calculate the weights and log-likelihood343    start_w = clock();344    wllk = weight_loglikelihood(par, betahat, varbetahat, ldscore, c0, Nstar, num_threads);345    w = as<mat>(wllk[""w""]);346    llk = as<long double>(wllk[""llk""]);347    finish_w = clock();348 349    // M-step: update the proportion parameters350    pi1 = update_pi1(w,Nstar,num_threads);351    finish_pi1 = clock();352 353    // M-step: update the variance parameters354    new_sig(0) = par(1); new_sig(1) = par(2);355    if(stratification==false){new_sig(1)=0;}356    for(int j=0; j<steps; j++){357      old_sig = onestep_varcomponent(new_sig, w, betahat,varbetahat, ldscore, Nstar, num_threads);358      if(stratification==false){old_sig(1)=0;}359      if((abs(old_sig(0)-new_sig(0))<1e-20) & (abs(old_sig(1)-new_sig(1))<1e-20)) break;360      new_sig(0) = old_sig(0); new_sig(1) = old_sig(1);361      if(new_sig(0) > 1) new_sig(0) = 1e-5;362      if(new_sig(0) < 0) new_sig(0) = 1e-12;363      if(new_sig(1) > 1) new_sig(1) = 1e-6;364      if(new_sig(1) < -min(varbetahat)) new_sig(1) = -min(varbetahat)/2;365    }366    finish_sigsq = clock();367 368    // update par vector with the new parameter values369    par(0) = pi1; par(1) = new_sig(0); par(2) = new_sig(1);370 371    if(isinf(-llk) | isnan(-llk) | isnan(llk)) {par(2) = abs(par(2));}372    if((llk<prev_llk) & (c0<20)) {c0 = c0+1;}373    if((llk<prev_llk) & (c0>=20)) {break;}374 375    // output results into a vector result376    result(0) = i; result(1) = llk;377    result(2) = par(0); result(3) = par(1); result(4) = par(2);378    result(5) = par(0)*par(1)*M; result(6) = c0;379    result(7) = double(finish_w-start_w)/CLOCKS_PER_SEC;380    result(8) = double(finish_pi1-finish_w)/CLOCKS_PER_SEC;381    result(9) = double(finish_sigsq-finish_pi1)/CLOCKS_PER_SEC;382 383    if(print==true){384      if(i%printfreq==0){385        for(int r=0; r<10; r++){386          Rcout << result(r)<< "", "" ;387        }388        Rcout << endl;389      }390    }391 392    error_pi = abs(prev_par(0) - par(0)) ;393    error_sigsq = abs(prev_par(1) - par(1)) ;394    error_a = abs(prev_par(2) - par(2)) ;395    increase_ll = (llk - prev_llk)/prev_llk ;396 397    if(((error_pi< eps1) & (error_sigsq <eps2) & (error_a <eps3)) & (abs(increase_ll) <eps) ) {break;}398  }399  // calcualte the logliklihood under the new par400  llk = loglikelihood(par,betahat,varbetahat,ldscore,c0,Nstar, num_threads);401  result(1) = llk;402  return result;403}404 405//--------------------------------406//--------------------------------407// [[Rcpp::export]]408vec Sk(const vec & par, const long double &betahatk, const long double &varbetahatk, const long double &ldscorek,const int & c0, const int &Nstark) // score vector for kth SNP 409{410  411  long double pic = par(0);412  long double sigsq = par(1);413  long double a = par(2);414  long double  tem_var, ww, numerator_pic=0, denominator=0, numerator_sigsq=0, numerator_a=0;415  vec res(3);416  int j;417  418  for(j=0; j<std::min(int(Nstark)+1, c0+1); j++){419    tem_var = (j*sigsq)*ldscorek/Nstark + varbetahatk + a;420    ww = cdbinom(j,Nstark,pic,false) * cdnorm(betahatk, 0, sqrt(tem_var), false);421    denominator += ww;422    423    //pic424    numerator_pic += ww*(j/pic - (Nstark-j)/(1-pic));425    426    //sigsq427    numerator_sigsq += ww*(-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0))*j*ldscorek/Nstark ;428    429    //a430    numerator_a += ww*(-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0));431    432  }433  434  res(0) = numerator_pic/denominator;435  res(1) = numerator_sigsq/denominator;436  res(2) = numerator_a/denominator;437  438  return res;439}440 441//--------------------------------442//--------------------------------443// [[Rcpp::export]]444mat Ik(const vec & par, const long double &betahatk, const long double &varbetahatk, const long double &ldscorek,const int & c0, const int &Nstark) // Information for kth SNP, i.e., -partial(Sk)/partial(theta)445{446  447  long double pic = par(0);448  long double sigsq = par(1);449  long double a = par(2);450  451  long double  tem_var, g, Lk=0;452  vec r(3),  dLk(3);453  dLk.fill(0.0);454  mat dr(3,3), rr(3,3), ddLk(3,3), dLk2(3,3), res(3,3);455  ddLk.fill(0.0);456  dLk2.fill(0.0);457  int j;458  459  for(j=0; j<std::min(int(Nstark)+1, c0+1); j++){460    461    tem_var = (j*sigsq)*ldscorek/Nstark + varbetahatk + a;462    463    g = cdbinom(j,Nstark,pic,false) * cdnorm(betahatk, 0, sqrt(tem_var), false);464    r(0) =  (j/pic - (Nstark-j)/(1-pic));465    r(1)= (-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0))*j*ldscorek/Nstark ;466    r(2) = (-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0));467    468    Lk += g;469    470    dr(0,0) = -j/pow(pic,2.0) - (Nstark-j)/(pow(1-pic, 2.0));471    dr(0,1) = 0;472    dr(0,2) = 0;473    dr(1,0) = 0;474    dr(1,1) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*pow(j*ldscorek/Nstark,2.0) ;475    dr(1,2) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*(j*ldscorek/Nstark) ;476    dr(2,0) = 0;477    dr(2,1) = dr(1,2);478    dr(2,2) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0));479    480    rr(0,0) = r(0)*r(0);481    rr(0,1) = r(0)*r(1);482    rr(0,2) = r(0)*r(2);483    rr(1,0) = r(1)*r(0);484    rr(1,1) = r(1)*r(1);485    rr(1,2) = r(1)*r(2);486    rr(2,0) = r(2)*r(0);487    rr(2,1) = r(2)*r(1);488    rr(2,2) = r(2)*r(2);489    490    ddLk += g*dr + g*rr;491    dLk += g*r;492  }493  494  dLk2(0,0) = dLk(0)*dLk(0);495  dLk2(0,1) = dLk(0)*dLk(1);496  dLk2(0,2) = dLk(0)*dLk(2);497  dLk2(1,0) = dLk(1)*dLk(0);498  dLk2(1,1) = dLk(1)*dLk(1);499  dLk2(1,2) = dLk(1)*dLk(2);500  dLk2(2,0) = dLk(2)*dLk(0);501  dLk2(2,1) = dLk(2)*dLk(1);502  dLk2(2,2) = dLk(2)*dLk(2);503  504  res = dLk2/pow(Lk,2.0) - ddLk/Lk;505  506  return res;507}508 509//--------------------------------510//--------------------------------511// [[Rcpp::export]]512vec S(const vec & par, const vec & betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // total score for all K SNPs for each parameter513{514  int K = betahat.n_elem;515  vec res(3);516  res.fill(0.0);517  vec tem(3);518  long double res0=0, res1=0, res2=0;519  int k;520  521  // --------------------------------//--------------------------------522  omp_set_num_threads(num_threads);523  #pragma omp parallel for shared(par, betahat, varbetahat,ldscore,Nstar,c0) private(tem,k) reduction(+:res0,res1,res2)524  // --------------------------------//--------------------------------525  for(k=0; k<K; k++){526    tem = Sk(par, (betahat(k)), (varbetahat(k)), (ldscore(k)), c0, (Nstar(k)));527    res0 += tem(0);528    res1 += tem(1);529    res2 += tem(2);530  }531  res(0) = res0;532  res(1) = res1;533  res(2) = res2;534  535  return res;536}537 538 539//--------------------------------540//--------------------------------541// [[Rcpp::export]]542mat SS(const vec & par, const vec & betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // store the score function for each SNP k, summarize it as a K*3 matrix543{544  int K = betahat.n_elem;545  mat res(K,3); res.fill(0.0);546  vec tem(3);547  int k;548  549  // --------------------------------//--------------------------------550  omp_set_num_threads(num_threads);551  #pragma omp parallel for shared(par, betahat, varbetahat,ldscore,Nstar,c0,res) private(k,tem) 552  // --------------------------------//--------------------------------553  for(k=0; k<K; k++){554    tem = Sk(par,  (betahat(k)), (varbetahat(k)), (ldscore(k)), c0, (Nstar(k)));555    res(k,0) = tem(0);556    res(k,1) = tem(1);557    res(k,2) = tem(2);558  }559  560  return res;561}562 563//--------------------------------564//--------------------------------565// [[Rcpp::export]]566mat I(const vec & par, const vec & betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // 3*3 Information matrix567{568  int K = betahat.n_elem;569  mat res(3,3), tem(3,3);570  res.fill(0.0);571  long double res00=0, res01=0, res02=0, res10=0,res11=0,res12=0,res20=0,res21=0,res22=0;572  int k;573  574  // --------------------------------//--------------------------------575  omp_set_num_threads(num_threads);576  #pragma omp parallel for shared(par, betahat, varbetahat,ldscore,Nstar,c0) private(k,tem) reduction(+:res00,res01,res02,res10,res11,res12,res20,res21,res22)577  // --------------------------------//--------------------------------578  for(k=0; k<K; k++){579    tem = Ik(par,  (betahat(k)), (varbetahat(k)), (ldscore(k)), c0, (Nstar(k)));580    res00 += tem(0,0); res01 += tem(0,1); res02 += tem(0,2);581    res10 += tem(1,0); res11 += tem(1,1); res12 += tem(1,2);582    res20 += tem(2,0); res21 += tem(2,1); res22 += tem(2,2);583  }584  585  res(0,0) = res00; res(0,1) = res01; res(0,2) = res02;586  res(1,0) = res10; res(1,1) = res11; res(1,2) = res12;587  res(2,0) = res20; res(2,1) = res21; res(2,2) = res22;588  589  return res;590}591 592//--------------------------------593//--------------------------------594// [[Rcpp::export]]595List mixture_components_marginal(const vec & par, const vec &ldscore, const int & c0, const vec &Nstar, const int &num_threads)596{597  int K = ldscore.n_elem;598  599  long double pic = par(0);600  long double sigsq = par(1);601  602  int k,j;603  mat ww(K, (c0+1)); mat vv(K, (c0+1));604  ww.fill(0.0);  vv.fill(0.0);605  606  // weights formula.607  omp_set_num_threads(num_threads);608  #pragma omp parallel for shared(sigsq,pic,ldscore,Nstar,K,c0,ww,vv) private(k,j)609  for(k=0; k<K; k++){610    for(j=0; j<std::min(int(Nstar(k))+1,c0+1); j++){611      ww(k, j) = cdbinom(j,Nstar(k),pic,false);612      vv(k, j) = (j*sigsq)*ldscore(k)/Nstar(k);613    }614  }615  616  return List::create(617    _[""proportions""] = ww,618    _[""varcomponents""] = vv619  );620  621}622 623 624 625//--------------------------------//--------------------------------626// 3-component model627//--------------------------------//--------------------------------628 629//--------------------------------630//--------------------------------631// [[Rcpp::export]]632long double loglikelihood3(const vec & par, const vec &betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // loglikelihood function633{634    int K = betahat.n_elem;635    long double loginside;636    637    long double pic = par(0);638    long double p0 = par(1);639    long double sig1 = par(2);640    long double sig2 = par(3);641    long double a = par(4);642    643    long double res = 0;644    long double y=0;645    long double tem=0;646    647    long double pi1 = pic*p0;648    long double pi2 = pic*(1-p0);649    vec tem_prob(3);650    tem_prob(0) = pi1;651    tem_prob(1) = pi2;652    tem_prob(2) = 1-pi1-pi2;653    int k,j1,j2;654    655    // -------*-------*-------*-------*-------*-------*-------*656    omp_set_num_threads(num_threads);657    #pragma omp parallel for shared(a,betahat,varbetahat,ldscore,Nstar,tem_prob) private(loginside,y,k,tem,j1,j2) reduction(+:res)658    // -------*-------*-------*-------*-------*-------*-------*659    for(k=0; k<K; k++){660        loginside = 0;661        y = 0;662        for(j1=0; j1<std::min(int(Nstar(k))+1,c0+1); j1++){663            for(j2=0; j2<std::min(int(Nstar(k))+1,c0+1); j2++){664                if(Nstar(k) - j1 - j2 <0) break;665                666                tem = cdmultinom3(j1, j2, Nstar(k)- j1 - j2,tem_prob,false);667                loginside += tem*cdnorm(betahat(k),0, sqrt( (j1*sig1+ j2*sig2 )*ldscore(k)/Nstar(k) + varbetahat(k) + a), false);668                y+=tem;669            }670        }671        res += log(loginside/y);672    }673    674    return res;675}676 677//--------------------------------678// weights function679//--------------------------------680// [[Rcpp::export]]681mat weight3(const vec & par, const vec &betahat, const vec &varbetahat, const vec &ldscore, const int & c0, const vec & Nstar, const int &num_threads) // E-step: weights function682{683    int K = betahat.n_elem;684    mat w(K, (c0+1)*(c0+1));685    w.fill(0.0);686    687    long double pic = par(0);688    long double p0 = par(1);689    long double sig1 = par(2);690    long double sig2 = par(3);691    long double a = par(4) ;692    693    long double pi1 = pic*p0;694    long double pi2 = pic*(1-p0);695    vec tem_prob(3);696    tem_prob(0) = pi1;697    tem_prob(1) = pi2;698    tem_prob(2) = 1-pi1-pi2;699    700    int k,j1,j2;701    vec te;702    703    // weights formula.704    // -------*-------*-------*-------*-------*-------*-------*705    omp_set_num_threads(num_threads);706    #pragma omp parallel for  shared(a,sig1,sig2,betahat,varbetahat,ldscore,Nstar,w,tem_prob,K,c0) private(k,j1,j2)707    // -------*-------*-------*-------*-------*-------*-------*708    for(k=0; k<K; k++){709        for(j1=0; j1<std::min(int(Nstar(k))+1,c0+1); j1++){710            for(j2=0; j2<std::min(int(Nstar(k))+1,c0+1); j2++){711                if(Nstar(k) - j1 - j2 < 0) break;712                w(k, j1*(c0+1)+j2) =cdmultinom3(j1, j2, Nstar(k)- j1 - j2,tem_prob,false)*cdnorm(betahat(k),0, sqrt( (j1*sig1+ j2*sig2 )*ldscore(k)/Nstar(k) + varbetahat(k) + a), false);713            }714        }715    }716    717    te = sum(w,1);718    te = pow(te, -1.0);719    w = (diagmat(te)) * w;720    return w;721    722}723 724//--------------------------------725//--------------------------------726// [[Rcpp::export]]727List weight_loglikelihood3(const vec & par, const vec &betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // loglikelihood function728{729  int K = betahat.n_elem;730  long double loginside;731  732  long double pic = par(0);733  long double p0 = par(1);734  long double sig1 = par(2);735  long double sig2 = par(3);736  long double a = par(4);737  738  long double res = 0;739  long double y=0;740  long double tem=0,tem1=0;741  742  mat w(K, (c0+1)*(c0+1));743  w.fill(0.0);744  vec te;745  746  long double pi1 = pic*p0;747  long double pi2 = pic*(1-p0);748  vec tem_prob(3);749  tem_prob(0) = pi1;750  tem_prob(1) = pi2;751  tem_prob(2) = 1-pi1-pi2;752  int k,j1,j2;753  754  // -------*-------*-------*-------*-------*-------*-------*755  omp_set_num_threads(num_threads);756  #pragma omp parallel for shared(a,betahat,varbetahat,ldscore,Nstar,tem_prob,w,K,c0) private(tem1,loginside,y,k,tem,j1,j2) reduction(+:res)757  // -------*-------*-------*-------*-------*-------*-------*758  for(k=0; k<K; k++){759    loginside = 0;760    y = 0;761    for(j1=0; j1<std::min(int(Nstar(k))+1,c0+1); j1++){762      for(j2=0; j2<std::min(int(Nstar(k))+1,c0+1); j2++){763        if(Nstar(k) - j1 - j2 <0) break;764        tem = cdmultinom3(j1, j2, Nstar(k)- j1 - j2,tem_prob,false);765        tem1 = tem*cdnorm(betahat(k),0, sqrt( (j1*sig1+ j2*sig2 )*ldscore(k)/Nstar(k) + varbetahat(k) + a), false);766        767        w(k, j1*(c0+1)+j2) =tem1; 768        loginside += tem1;769        y+=tem;770      }771    }772    res += log(loginside/y);773  }774  775  te = sum(w,1);776  te = pow(te, -1.0);777  w = (diagmat(te)) * w;778  779  return List::create(780    _[""w""] = w,781    _[""llk""] = res782  );783}784 785 786//--------------------------------787//--------------------------------788// [[Rcpp::export]]789vec update_p3(const mat & w, const vec & Nstar, const int &num_threads)790{791    792    long double res0=0;793    long double res1=0;794    long double res2=0;795    long double tem;796    vec res(2);797    int c0 = sqrt(w.n_cols) - 1;798    int K = Nstar.n_elem;799    int k,j1,j2;800    801    // -------*-------*-------*-------*-------*-------*-------*802    omp_set_num_threads(num_threads);803    #pragma omp parallel for shared(Nstar,c0,w) private(k,tem,j1,j2) reduction(+:res0,res1,res2)804    // -------*-------*-------*-------*-------*-------*-------*805    for(k=0; k<K; k++){806        for(j1=0; j1<std::min(int(Nstar(k))+1,c0+1); j1++){807            for(j2=0; j2<std::min(int(Nstar(k))+1,c0+1); j2++){808                if(Nstar(k) - j1 - j2 <0) break;809                tem = w(k, j1*(c0+1)+j2);810                res0 += tem*(j1+j2);811                res1 += tem*Nstar(k);812                res2 += tem*j1;813            }814        }815    }816    817    res(0) = res0/res1;818    res(1) = res2/res0;819    820    return res;821}822 823//--------------------------------824//--------------------------------825// [[Rcpp::export]]826vec onestep_varcomponent3(const vec &varcomponent, const mat & w, const vec &betahat, const vec &varbetahat, const vec &ldscore, const vec & Nstar, const int &num_threads)827{828    int K = betahat.n_elem;829    int c0 = sqrt(w.n_cols) - 1;830    831    long double sig1 = varcomponent(0);832    long double sig2 = varcomponent(1);833    long double a = varcomponent(2);834    835    long double tem, det;836    long double s1 = 0, s2=0, s3 =0, dd11 = 0, dd22=0,dd33=0, dd12=0, dd13=0, dd23=0, dd123=0;837    int k,j1,j2;838    vec result(3);839    840    // -------*-------*-------*-------*-------*-------*-------*841    omp_set_num_threads(num_threads);842    #pragma omp parallel for shared(Nstar,a,sig1,sig2,ldscore,varbetahat,betahat,w) private(tem,k,j1,j2) reduction(+:s1,s2,s3,dd11,dd22,dd33,dd12,dd13,dd23,dd123)843    // -------*-------*-------*-------*-------*-------*-------*844    for(k=0; k<K; k++)845    {846        for(j1=0; j1<std::min(int(Nstar(k))+1,c0+1); j1++){847            for(j2=0; j2<std::min(int(Nstar(k))+1,c0+1); j2++){848                if(Nstar(k) - j1 - j2 <0) break;849                tem = (j1*sig1+j2*sig2) * ldscore(k) /double(Nstar(k)) + varbetahat(k) + a;850                s1 += w(k, j1*(c0+1)+j2)*j1*ldscore(k)/(2.0*Nstar(k))*(pow(betahat(k),2)/pow(tem,2) - 1.0/(tem)  );851                s2 += w(k, j1*(c0+1)+j2)*j2*ldscore(k)/(2.0*Nstar(k))*(pow(betahat(k),2)/pow(tem,2) - 1.0/(tem)  );852                s3 += w(k, j1*(c0+1)+j2)/(2.0)*(pow(betahat(k),2)/pow(tem,2) - 1.0/(tem)  );853                854                dd11 += w(k, j1*(c0+1)+j2)*j1*ldscore(k)*j1*ldscore(k)/(2.0*Nstar(k)*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));855                dd22 += w(k, j1*(c0+1)+j2)*j2*ldscore(k)*j2*ldscore(k)/(2.0*Nstar(k)*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));856                dd33 += w(k, j1*(c0+1)+j2)/(2.0)*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));857                dd12 += w(k, j1*(c0+1)+j2)*j1*ldscore(k)*j2*ldscore(k)/(2.0*Nstar(k)*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));858                dd13 += w(k, j1*(c0+1)+j2)*j1*ldscore(k)/(2.0*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));859                dd23 += w(k, j1*(c0+1)+j2)*j2*ldscore(k)/(2.0*Nstar(k))*(-2.0*pow(betahat(k),2)/pow(tem,3) + 1.0/(pow(tem,2)));860                dd123 += w(k, j1*(c0+1)+j2)*j1*ldscore(k)*j2*ldscore(k)/(2.0*Nstar(k)*Nstar(k))*(6.0*pow(betahat(k),2)/pow(tem,4) - 2.0/(pow(tem,3)));861            }862        }863    }864    865    det = dd11*(dd33*dd22 - dd23*dd23) - dd12*(dd33*dd12-dd23*dd13) + dd13*(dd23*dd12 - dd22*dd13);866    867    result(0) = sig1 - (1.0/det)*( (dd33*dd22-dd23*dd23)*s1-(dd33*dd12-dd23*dd13)*s2 + (dd23*dd12-dd22*dd13)*s3  );868    result(1) = sig2 -  (1.0/det)*( -(dd33*dd12-dd13*dd23)*s1 + (dd33*dd11-dd13*dd13)*s2 - (dd23*dd11-dd12*dd13)*s3  );869    result(2) = a - (1.0/det)*( (dd23*dd12-dd13*dd22)*s1 - (dd23*dd11-dd13*dd12)*s2 + (dd22*dd11-dd12*dd12)*s3  );870    871    return result;872}873 874//--------------------------------875//--------------------------------876// [[Rcpp::export]]877vec EM_func3(const vec &par_start, const vec & lower_pi, const vec & upper_pi,878             const vec &betahat, const vec & varbetahat, const vec & ldscore,  const vec & Nstar, const int & M,879             int c0,const long double &eps1,const long double &eps2,const long double &eps3,const long double &eps4,const long double &eps5,880             const long double &eps, const int &Meps, const int &steps, const int &num_threads, const bool &print, const int &printfreq, const bool &stratification)881{882    vec par(5), prev_par(5);883    par(0) = par_start(0);884    par(1) = par_start(1);885    par(2) = par_start(2);886    par(3) = par_start(3);887    par(4) = par_start(4);888    889    long double llk = 0,error_pi, error_p0, error_sig1, error_sig2, error_a, increase_ll;890    long double prev_llk, pic, p0, sig1, sig2, a;891    clock_t start_w, finish_w, finish_p, finish_sigsq; 892    893    List wllk;894    mat w;895    vec tem_p(2);896    vec tem_sig(3);897    vec old(3);898    vec result(12);899    900    Rcout << ""Iteration, prev_loglikelihood, pic, p1, sigmasq1, sigmasq2, a, Heritability, c0, Seconds(weight_llk), Seconds(pi), Seconds(variance_components)"";901    Rcout << endl;902    903    for(int i=0; i<Meps; i++){904        prev_llk = llk;905        prev_par(0) = par(0); prev_par(1) = par(1); prev_par(2) = par(2); prev_par(3) = par(3);prev_par(4) = par(4);906        907        pic = par(0);908        p0  = par(1);909        sig1 = par(2);910        sig2 = par(3);911        a = par(4);912        913        // update weight and log-likelihood914        start_w = clock();915        wllk = weight_loglikelihood3(par, betahat, varbetahat, ldscore,c0,Nstar,num_threads);916        w = as<mat>(wllk[""w""]);917        llk = as<long double>(wllk[""llk""]);918        finish_w = clock();919        920        // update proportion921        if ((pic>=lower_pi(0)) & (pic<=upper_pi(0)) & (p0>=lower_pi(1)) & (p0 <= upper_pi(1))) tem_p = update_p3(w,Nstar,num_threads);922        finish_p = clock(); 923        924        if (pic < lower_pi(0)) tem_p(0) = lower_pi(0);925        if (p0  < lower_pi(1)) tem_p(1) = lower_pi(1);926        if (pic > upper_pi(0)) tem_p(0) = upper_pi(0);927        if (p0  > upper_pi(1)) tem_p(1) = upper_pi(1);928        929        // update variance components930        tem_sig(0) = sig1; tem_sig(1) = sig2; tem_sig(2) = a;931        if(stratification==false){tem_sig(2)=0;}932        for(int j=0; j<steps; j++){933            old = onestep_varcomponent3(tem_sig, w, betahat,varbetahat, ldscore, Nstar,num_threads);934            if(stratification==false){old(2)=0;}935            if((abs(old(0)-tem_sig(0))<1e-20) & (abs(old(1)-tem_sig(1))<1e-20) & (abs(old(2)-tem_sig(2))<1e-20) ) break;936            tem_sig(0) = old(0); tem_sig(1) = old(1); tem_sig(2) = old(2);937        }938        if(tem_sig(0) > 1) tem_sig(0) = 1e-5;939        if(tem_sig(0) < 0) tem_sig(0) = 1e-12;940        if(tem_sig(1) > 1) tem_sig(1) = 1e-5;941        if(tem_sig(1) < 0) tem_sig(1) = 1e-12;942        if(tem_sig(2) > 1) tem_sig(2) = 1e-5;943        if(tem_sig(2) < -min(varbetahat)) tem_sig(2) = -min(varbetahat)/2;944        finish_sigsq = clock();945        946        par(0) = tem_p(0); par(1) =tem_p(1);947        par(2) = tem_sig(0); par(3) = tem_sig(1); par(4) = tem_sig(2);948        949        if(isinf(-llk) | isnan(-llk) | isnan(llk)) {par(4) = abs(par(4));}950        if((llk<prev_llk) & (c0<20)) {c0 = c0+1;}951        if((llk<prev_llk) & (c0>=20)) {break;}952        953        result(0) = i; result(1) = llk;954        result(2) = par(0); result(3) = par(1); result(4) = par(2); result(5) = par(3); result(6)= par(4);955        result(7) = M*par(0)*( par(1)*par(2) + (1-par(1))*par(3)); result(8) = c0;956        957        result(9) = double(finish_w-start_w)/CLOCKS_PER_SEC;958        result(10) = double(finish_p-finish_w)/CLOCKS_PER_SEC;959        result(11) = double(finish_sigsq-finish_p)/CLOCKS_PER_SEC;960 961        if(print==true){962          if(i%printfreq==0){963            for(int r=0; r<12; r++){964              Rcout << result(r)<< "", "" ;965            }966            Rcout << endl;967          }968        }969        970        error_pi = abs(prev_par(0) - par(0)) ;971        error_p0 = abs(prev_par(1) - par(1)) ;972        error_sig1 = abs(prev_par(2) - par(2)) ;973        error_sig2 = abs(prev_par(3) - par(3)) ;974        error_a = abs(prev_par(4) - par(4)) ;975        increase_ll = (llk - prev_llk)/prev_llk;976        977        if(((error_pi< eps1) & (error_p0<eps2) & (error_sig1 <eps3) & (error_sig2<eps4) & (error_a<eps5)) & (abs(increase_ll) <eps)){break;}978    }979    980    // update log-likelihood981    llk = loglikelihood3(par,betahat,varbetahat,ldscore,c0,Nstar,num_threads);982    result(1) = llk;983    984    return result;985    986}987 988 989//--------------------------------990//--------------------------------991// [[Rcpp::export]]992vec Sk3(const vec & par, const long double &betahatk, const long double &varbetahatk, const long double &ldscorek,const int & c0, const int &Nstark) // score vector for kth SNP993{994    995    long double pic = par(0);996    long double p1 = par(1);997    long double sig1sq = par(2);998    long double sig2sq = par(3);999    long double a = par(4);1000    1001    vec tem_prob(3);1002    tem_prob(0) = pic*p1;1003    tem_prob(1) = pic*(1-p1);1004    tem_prob(2) = 1-pic;1005    int j1,j2;1006    1007    long double ww, denominator=0,1008    numerator_pic=0,numerator_p1=0,numerator_sig1sq=0,1009    numerator_sig2sq=0, numerator_a=0;1010    vec res(5);1011    long double tem_var;1012    1013    for(j1=0; j1<std::min(Nstark+1,c0+1); j1++){1014        for(j2=0; j2<std::min(Nstark+1,c0+1); j2++){1015            if(Nstark - j1 - j2 <0) break;1016            tem_var = (j1*sig1sq+j2*sig2sq)*ldscorek/Nstark + varbetahatk + a;1017            ww = cdmultinom3(j1, j2, Nstark- j1 - j2,tem_prob,false)*cdnorm(betahatk,0, sqrt( tem_var), false);1018            denominator += ww;1019            1020            //pic1021            numerator_pic += ww*((j1+j2)/pic - (Nstark-j1-j2)/(1-pic));1022            numerator_p1 += ww*(j1/p1 - j2/(1-p1));1023            1024            //sigma1^21025            numerator_sig1sq += ww*(-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0))*j1*ldscorek/Nstark ;1026            numerator_sig2sq += ww*(-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0))*j2*ldscorek/Nstark ;1027            1028            //sigma1^21029            numerator_a += ww*(-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0));1030            1031        }1032    }1033    res(0) = numerator_pic/denominator;1034    res(1) = numerator_p1/denominator;1035    res(2) = numerator_sig1sq/denominator;1036    res(3) = numerator_sig2sq/denominator;1037    res(4) = numerator_a/denominator;1038    1039    return res;1040}1041 1042 1043 1044 1045 1046//--------------------------------1047//--------------------------------1048// [[Rcpp::export]]1049mat Ik3(const vec & par, const long double &betahatk, const long double &varbetahatk, const long double &ldscorek,const int & c0, const int &Nstark) // Information for kth SNP, i.e., -partial(Sk)/partial(theta)1050{1051    long double pic = par(0);1052    long double p1 = par(1);1053    long double sig1sq = par(2);1054    long double sig2sq = par(3);1055    long double a = par(4);1056    1057    vec tem_prob(3);1058    tem_prob(0) = pic*p1;1059    tem_prob(1) = pic*(1-p1);1060    tem_prob(2) = 1-pic;1061    int j1,j2;1062    1063    long double  tem_var, g, Lk=0;1064    vec r(5),  dLk(5);1065    dLk.fill(0.0); r.fill(0.0);1066    mat dr(5,5), rr(5,5), ddLk(5,5), dLk2(5,5), res(5,5);1067    ddLk.fill(0.0);1068    dLk2.fill(0.0);1069    dr.fill(0.0);1070    rr.fill(0.0);1071    res.fill(0.0);1072    1073    for(j1=0; j1<std::min(Nstark+1,c0+1); j1++){1074        for(j2=0; j2<std::min(Nstark+1,c0+1); j2++){1075            if(Nstark - j1 - j2 <0) break;1076            1077            tem_var = (j1*sig1sq+j2*sig2sq)*ldscorek/Nstark + varbetahatk + a;1078            g = cdmultinom3(j1, j2, Nstark- j1 - j2,tem_prob,false)*cdnorm(betahatk,0, sqrt(tem_var), false);1079            1080            r(0) = (j1+j2)/pic - (Nstark-j1-j2)/(1-pic);1081            r(1) = j1/p1 - j2/(1-p1);1082            r(2) = (-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0))*j1*ldscorek/Nstark ;1083            r(3) = (-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0))*j2*ldscorek/Nstark ;1084            r(4) = (-0.5/tem_var + 0.5*pow(betahatk,2.0)/pow(tem_var,2.0));1085            1086            Lk += g;1087            1088            dr(0,0) = -(j1+j2)/pow(pic,2.0) - (Nstark-j1-j2)/(pow(1-pic, 2.0));1089            dr(1,1) = -(j1)/pow(p1,2.0) - (j2)/(pow(1-p1, 2.0));1090            dr(2,2) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*(j1*ldscorek/Nstark)*(j1*ldscorek/Nstark);1091            dr(2,3) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*(j1*ldscorek/Nstark)*(j2*ldscorek/Nstark);1092            dr(2,4) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*(j1*ldscorek/Nstark);1093            dr(3,2) = dr(2,3);1094            dr(3,3) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*(j2*ldscorek/Nstark)*(j2*ldscorek/Nstark);1095            dr(3,4) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0))*(j2*ldscorek/Nstark);1096            dr(4,2) = dr(2,4);1097            dr(4,3) = dr(3,4);1098            dr(4,4) = (0.5/pow(tem_var,2.0) - pow(betahatk,2.0)/pow(tem_var,3.0));1099            1100            rr(0,0) = r(0)*r(0); rr(0,1) = r(0)*r(1); rr(0,2) = r(0)*r(2); rr(0,3) = r(0)*r(3); rr(0,4) = r(0)*r(4);1101            rr(1,0) = rr(0,1);   rr(1,1) = r(1)*r(1); rr(1,2) = r(1)*r(2); rr(1,3) = r(1)*r(3); rr(1,4) = r(1)*r(4);1102            rr(2,0) = rr(2,0);   rr(2,1) = rr(1,2);   rr(2,2) = r(2)*r(2); rr(2,3) = r(2)*r(3); rr(2,4) = r(2)*r(4);1103            rr(3,0) = rr(0,3);   rr(3,1) = rr(1,3);   rr(3,2) = rr(2,3);   rr(3,3) = r(3)*r(3); rr(3,4) = r(3)*r(4);1104            rr(4,0) = rr(0,4);   rr(4,1) = rr(1,4);   rr(4,2) = rr(2,4);   rr(4,3) = rr(3,4);   rr(4,4) = r(4)*r(4);1105            1106            ddLk += g*dr + g*rr;1107            dLk += g*r;1108        }1109    }1110    1111    dLk2(0,0) = dLk(0)*dLk(0);dLk2(0,1) = dLk(0)*dLk(1);dLk2(0,2) = dLk(0)*dLk(2);dLk2(0,3) = dLk(0)*dLk(3);dLk2(0,4) = dLk(0)*dLk(4);1112    dLk2(1,0) = dLk2(0,1);    dLk2(1,1) = dLk(1)*dLk(1);dLk2(1,2) = dLk(1)*dLk(2);dLk2(1,3) = dLk(1)*dLk(3);dLk2(1,4) = dLk(1)*dLk(4);1113    dLk2(2,0) = dLk2(0,2);    dLk2(2,1) = dLk2(1,2);    dLk2(2,2) = dLk(2)*dLk(2);dLk2(2,3) = dLk(2)*dLk(3);dLk2(2,4) = dLk(2)*dLk(4);1114    dLk2(3,0) = dLk2(0,3);    dLk2(3,1) = dLk2(1,3);    dLk2(3,2) = dLk2(2,3);    dLk2(3,3) = dLk(3)*dLk(3);dLk2(3,4) = dLk(3)*dLk(4);1115    dLk2(4,0) = dLk2(0,4);    dLk2(4,1) = dLk2(1,4);    dLk2(4,2) = dLk2(2,4);    dLk2(4,3) = dLk2(3,4);    dLk2(4,4) = dLk(4)*dLk(4);1116    1117    res = dLk2/pow(Lk,2.0) - ddLk/Lk;1118    1119    return res;1120}1121 1122//--------------------------------1123//--------------------------------1124// [[Rcpp::export]]1125vec S3(const vec & par, const vec & betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar) // total score for all K SNPs for each parameter1126{1127    int K = betahat.n_elem;1128    vec res(5);1129    vec tem(5);1130    res.fill(0.0);1131    1132    for(int k=0; k<K; k++){1133        tem = Sk3(par,  (betahat(k)), (varbetahat(k)), (ldscore(k)), c0, (Nstar(k)))/K;1134        res(0) += tem(0);1135        res(1) += tem(1);1136        res(2) += tem(2);1137        res(3) += tem(3);1138        res(4) += tem(4);1139    }1140    return res;1141}1142 1143//--------------------------------1144//--------------------------------1145// [[Rcpp::export]]1146mat SS3(const vec & par, const vec & betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // store the score function for each SNP k, summarize it as a K*3 matrix1147{1148    int K = betahat.n_elem;1149    mat res(K,5); res.fill(0.0);1150    vec tem(5);1151    int k;1152    1153    // --------------------------------//--------------------------------1154    omp_set_num_threads(num_threads);1155    #pragma omp parallel for shared(par, betahat, varbetahat,ldscore,Nstar,c0,res) private(k,tem)1156    // --------------------------------//--------------------------------1157    for(k=0; k<K; k++){1158        tem = Sk3(par,  (betahat(k)), (varbetahat(k)), (ldscore(k)), c0, (Nstar(k)));1159        res(k,0) = tem(0);1160        res(k,1) = tem(1);1161        res(k,2) = tem(2);1162        res(k,3) = tem(3);1163        res(k,4) = tem(4);1164    }1165    1166    return res;1167}1168 1169 1170//--------------------------------1171//--------------------------------1172// [[Rcpp::export]]1173mat I3(const vec & par, const vec & betahat, const vec &varbetahat, const vec &ldscore,const int & c0, const vec &Nstar, const int &num_threads) // 3*3 Information matrix1174{1175    int K = betahat.n_elem;1176    mat res(5,5), tem(5,5);1177    res.fill(0.0), tem.fill(0.0);1178    long double res00=0,res01=0,res02=0,res03=0,res04=0,res10=0,res11=0,res12=0,res13=0,res14=0,res20=0,res21=0,res22=0,res23=0,res24=0;1179    long double res30=0,res31=0,res32=0,res33=0,res34=0,res40=0,res41=0,res42=0,res43=0,res44=0;1180    int k;1181    1182    // --------------------------------//--------------------------------1183    omp_set_num_threads(num_threads);1184    #pragma omp parallel for shared(par, betahat, varbetahat,ldscore,Nstar,c0) private(k,tem) reduction(+:res00,res01,res02,res03,res04,res10,res11,res12,res13,res14,res20,res21,res22,res23,res24,res30,res31,res32,res33,res34,res40,res41,res42,res43,res44)1185    // --------------------------------//--------------------------------1186    for(k=0; k<K; k++){1187      tem = Ik3(par,  (betahat(k)), (varbetahat(k)), (ldscore(k)), c0, (Nstar(k)));1188      res00 += tem(0,0); res01 += tem(0,1); res02 += tem(0,2); res03 += tem(0,3); res04 += tem(0,4); 1189      res10 += tem(1,0); res11 += tem(1,1); res12 += tem(1,2); res13 += tem(1,3); res14 += tem(1,4); 1190      res20 += tem(2,0); res21 += tem(2,1); res22 += tem(2,2); res23 += tem(2,3); res24 += tem(2,4); 1191      res30 += tem(3,0); res31 += tem(3,1); res32 += tem(3,2); res33 += tem(3,3); res34 += tem(3,4); 1192      res40 += tem(4,0); res41 += tem(4,1); res42 += tem(4,2); res43 += tem(4,3); res44 += tem(4,4); 1193    }1194    1195    res(0,0) = res00; res(0,1) = res01; res(0,2) = res02; res(0,3) = res03; res(0,4) = res04;1196    res(1,0) = res10; res(1,1) = res11; res(1,2) = res12; res(1,3) = res13; res(1,4) = res14;1197    res(2,0) = res20; res(2,1) = res21; res(2,2) = res22; res(2,3) = res23; res(2,4) = res24;1198    res(3,0) = res30; res(3,1) = res31; res(3,2) = res32; res(3,3) = res33; res(3,4) = res34;1199    res(4,0) = res40; res(4,1) = res41; res(4,2) = res42; res(4,3) = res43; res(4,4) = res44;1200    

Showing the first 1,200 of 42633 lines. Download the file for the rest.