#include <nufamily.h>
Definition at line 76 of file nufamily.h.
|
|
Definition at line 63 of file nufamily.cpp. References in_connector, kernal, myfather, mymother, myoffspring, numcut, numoffs, out_connector, pop, seq_id, terminal, tn_genotype, and tn_qtl.
00064 {
00065 seq_id = 0;
00066 kernal = 0;
00067 terminal = 0;
00068 numoffs = 0;
00069 numcut = 0;
00070 myfather = 0;
00071 mymother = 0;
00072 myoffspring = 0;
00073 tn_genotype = 0;
00074 tn_qtl = 4;
00075 pop = 0;
00076 in_connector = 0;
00077 out_connector = 0;
00078 }
|
|
|
Definition at line 93 of file nufamily.h. References release().
00093 {release();}
|
|
||||||||||||||||
|
Definition at line 388 of file nufamily.cpp. References matvec::Individual::anterior, matvec::Matrix< double >::begin(), father_indx, fullsibs_prob(), matvec::Individual::get_penetrance(), get_posterior(), mother_indx, myfather, mymother, matvec::Individual::nspouse(), numoffs, and tn_genotype. Referenced by terminal_peeling().
00389 {
00390 /////////////////////////////////////////////////
00391 // I must be one of offspring in the NuFamily
00392 //////////////////////////////////////////////////
00393 if (!trans) throw exception("NuFamily::anterior(arg1,arg2,arg3): null arg2");
00394 Vector<double> post_vec(tn_genotype);
00395 Vector<double> pen_vec(tn_genotype);
00396 Vector<double> apm(tn_genotype);
00397 Vector<double> apf(tn_genotype);
00398
00399 unsigned um,uf,ui;
00400 double ai,val,**trme,scale = 0.0;
00401 if ( numoffs > 1 ) {
00402 scale = fullsibs_prob(I,trans,WSpace);
00403 }
00404 scale += mymother->anterior[tn_genotype];
00405 for (um=0; um<tn_genotype; um++) apm[um] = mymother->anterior[um];
00406 if (mymother->nspouse() > 1) {
00407 //new version
00408 scale += get_posterior(mymother,father_indx,post_vec);
00409 //old version
00410 // scale += get_posterior(mymother,mates_indx[1],post_vec);
00411 for (um=0; um<tn_genotype; um++) apm[um] *= post_vec[um];
00412 }
00413
00414 scale += myfather->anterior[tn_genotype];
00415 for (uf=0; uf<tn_genotype; uf++) apf[uf] = myfather->anterior[uf];
00416 if (myfather->nspouse() > 1) {
00417 scale += get_posterior(myfather,mother_indx,post_vec); //new version
00418
00419 // scale += get_posterior(myfather,mates_indx[0],post_vec); //old version
00420 for (uf=0; uf<tn_genotype; uf++) apf[uf] *= post_vec[uf];
00421 }
00422
00423 scale += I->get_penetrance(pen_vec);
00424 if (numoffs > 1) {
00425 for (val=0.0,ui=0; ui<tn_genotype; ui++) {
00426 trme = trans[ui].begin();
00427 for (ai=0.0,um=0; um<tn_genotype; um++) {
00428 for(uf=0; uf<tn_genotype; uf++) {
00429 ai += apm[um]*trme[um][uf]*WSpace[um][uf]*apf[uf];
00430 }
00431 }
00432 val += I->anterior[ui] = ai*pen_vec[ui];
00433 }
00434 }
00435 else {
00436 for (val=0.0,ui=0; ui<tn_genotype; ui++) {
00437 trme = trans[ui].begin();
00438 for (ai=0.0,um=0; um<tn_genotype; um++) {
00439 for(uf=0; uf<tn_genotype; uf++) ai += apm[um]*trme[um][uf]*apf[uf];
00440 }
00441 val += I->anterior[ui] = ai*pen_vec[ui];
00442 }
00443 }
00444 for (ui=0; ui<tn_genotype; ui++) I->anterior[ui] /= val;
00445 scale += std::log(val);
00446 I->anterior[tn_genotype] = scale;
00447 }
|
|
|
Definition at line 241 of file nufamily.cpp. References connectors, father_gid(), matvec::Individual::gid(), mother_gid(), myfather, mymother, myoffspring, and numoffs. Referenced by matvec::Population::peeling_sequence().
00242 {
00243 connectors.clear();
00244 if (father_gid() > 1) connectors.push_back(myfather);
00245 if (mother_gid() > 1) connectors.push_back(mymother);
00246 for (unsigned j=0; j<numoffs; j++) {
00247 if (myoffspring[j]->gid() > 1) connectors.push_back(myoffspring[j]);
00248 }
00249 return connectors.size();
00250 }
|
|
|
Definition at line 263 of file nufamily.cpp. References matvec::Population::analysis_type, matvec::Individual::anterior_iw, matvec::Vector< T >::begin(), matvec::Individual::connect_keeper, connectors, father_indx, matvec::Population::get_genotype_freq(), matvec::Individual::group_id, matvec::Individual::id(), matvec::Individual::initial_anterior(), matvec::Individual::initial_multi_anterior(), matvec::Individual::initial_multi_m_anterior(), matvec::Individual::initial_multi_m_posterior(), matvec::Individual::initial_multi_posterior(), matvec::Individual::initial_posterior(), matvec::Individual::loop_connector, mother_indx, matvec::Individual::myfather, myfather, matvec::Individual::mymother, mymother, myoffspring, matvec::Individual::myoffspring, numcut, matvec::Individual::numoffs, numoffs, matvec::Individual::numoffs_spouse, matvec::Individual::numspouse, matvec::Individual::offs_tree, matvec::Population::P_ndim, matvec::Individual::p_origin, pop, matvec::Individual::posterior, matvec::Individual::posterior_iw, matvec::Individual::related_family, matvec::Individual::reset_id(), matvec::Population::size(), matvec::Individual::spouselist, tn_genotype, and tn_qtl. Referenced by matvec::Population::break_loop().
00263 {
00264 unsigned damid = mymother->id();
00265 unsigned sireid = myfather->id();
00266 if (damid > pop->size()) damid -= pop->size();
00267 if (sireid > pop->size()) sireid -= pop->size();
00268 //std::cout << " " << pop->ind_name(unwelcome->id())
00269 // << " " << pop->ind_name(damid) << " " << pop->ind_name(sireid) << std::endl;
00270 //Vector<double> genotype_freq = pop->get_genotype_freq();
00271 Individual *I = new Individual(*unwelcome);
00272 std::string analysis_type=pop->analysis_type; //BRS
00273 unwelcome->group_id -= 1;
00274 unwelcome->connect_keeper = 1;
00275 I->reset_id(pop->size()+unwelcome->id());
00276 I->group_id = 1;
00277 I->loop_connector = 0;
00278 I->connect_keeper = 0;
00279 I->related_family.clear();
00280 I->offs_tree.clear();
00281
00282 std::list<Individual *>::iterator pos;
00283 pos=connectors.begin(); //BRS flag
00284 while (pos != connectors.end()) {
00285 if ((*pos) == unwelcome) { pos = connectors.erase(pos);}
00286 else { pos++; } //BRS
00287
00288 }
00289 if (I->numoffs_spouse) {
00290 delete [] I->numoffs_spouse;
00291 I->numoffs_spouse = 0;
00292 }
00293 if (I->spouselist) {
00294 delete [] I->spouselist;
00295 I->spouselist = 0;
00296 I->numspouse = 0;
00297 }
00298 if (I->myoffspring) {
00299 delete [] I->myoffspring;
00300 I->myoffspring = 0;
00301 // I->numoffs = numoffs; // however, I->myoffspring is not needed
00302 }
00303 I->anterior_iw = 0;
00304 if (analysis_type == "multipoint") { //BRS added this if statement
00305 doubleMatrix pen(pop->P_ndim,4);
00306 I->initial_multi_anterior(pen);
00307 if (I->posterior) {
00308 delete [] I->posterior;
00309 I->posterior = 0;
00310 I->initial_multi_posterior(-1);
00311 //need to delete m_posterior here also but not sure how.
00312 }
00313 }
00314 else if (analysis_type == "multipoint_m") { //BRS added this if statement
00315 doubleMatrix pen(pop->P_ndim,4);
00316 I->initial_multi_m_anterior(tn_qtl);
00317 if (I->posterior) {
00318 delete [] I->posterior;
00319 I->posterior = 0;
00320 I->initial_multi_m_posterior(tn_qtl,-1);
00321 //need to delete m_posterior here also but not sure how.
00322 }
00323 }
00324 else {
00325 Vector<double> genotype_freq = pop->get_genotype_freq(); // BRS moved this here
00326 I->initial_anterior(genotype_freq.begin(),tn_genotype);
00327 if (I->posterior) {
00328 delete [] I->posterior;
00329 I->posterior = 0;
00330 }
00331 }
00332 if (myfather == unwelcome) {
00333 I->myfather = 0;
00334 I->mymother = 0;
00335 I->numspouse = 1;
00336 I->spouselist = new Individual*[1];
00337 I->spouselist[0] = mymother;
00338 I->p_origin = 0; // phenotype is extended (not original)
00339 if (analysis_type == "multipoint") { //BRS added this if statement
00340 I->initial_multi_posterior(-1);
00341 }
00342 if (analysis_type == "multipoint_m") { //BRS added this if statement
00343 I->initial_multi_m_posterior(tn_qtl,-1);
00344 }
00345 else {
00346 I->initial_posterior(tn_genotype);
00347 }
00348 myfather = I;
00349 unwelcome->posterior_iw[mother_indx] = 3; // unwelcome has extended ped
00350 mother_indx = 0;
00351 return;
00352 }
00353 if (mymother == unwelcome) {
00354 I->myfather = 0;
00355 I->mymother = 0;
00356 I->numspouse = 1;
00357 I->spouselist = new Individual*[1];
00358 I->spouselist[0] = myfather;
00359 I->p_origin = 0; // phenotype is extended (not original)
00360 if (analysis_type == "multipoint") { //BRS added this if statement
00361 I->initial_multi_posterior(-1);
00362 }
00363 if (analysis_type == "multipoint_m") { //BRS added this if statement
00364 I->initial_multi_m_posterior(tn_qtl,-1);
00365 }
00366 else {
00367 I->initial_posterior(tn_genotype);
00368 }
00369 mymother = I;
00370 unwelcome->posterior_iw[father_indx] = 3; // unwelcome has extended ped
00371 father_indx = 0;
00372 return;
00373 }
00374 for (unsigned i=0; i<numoffs; i++) {
00375 if (myoffspring[i] == unwelcome) {
00376 numcut++;
00377 I->numoffs = 0;
00378 I->p_origin = 1; // new offspring keeps original phenotype
00379 unwelcome->p_origin = 0; // the connector's phenotype => extended
00380 myoffspring[i] = I;
00381 unwelcome->anterior_iw = 3; // unwelcome has extended ped
00382 return;
00383 }
00384 }
00385 throw exception("NuFamily::cutting(unwelcome): you have probably found a bug!");
00386 }
|
|
|
Definition at line 553 of file nufamily.cpp. References father_id(), matvec::Population::ind_name(), mother_id(), pop, and matvec::Population::size().
00553 {
00554
00555 int ident;
00556 int popsize=pop->size();
00557
00558 std::cout << " Father " << pop->ind_name(father_id());
00559 std::cout << " Mother " << pop->ind_name(mother_id());
00560 /*
00561 for (unsigned i=0; i<numoffs; i++) {
00562 ident=myoffspring[i]->id();
00563 std::cout << " child" << i+1;
00564 if (ident > popsize) {
00565 // this is a cut individual so need to get the original one
00566 ident -= popsize; // correct id;
00567 std::cout << " (cut) " << pop->ind_name(ident);
00568 }
00569 else
00570 std::cout << " " << pop->ind_name(ident);
00571 // std::cout << " or " << pop->ind_name(myoffspring[i]->id());
00572 }
00573 */
00574 std::cout << std::endl;
00575 }
|
|
|
|
|
|
Definition at line 2319 of file nufamily.cpp. References matvec::Individual::genotNodeVector, matvec::MaternalPaternalAlleles::maternal, myfather, matvec::Individual::myid, mymother, myoffspring, numoffs, and matvec::MaternalPaternalAlleles::paternal. Referenced by matvec::Population::LGGenotypeElimination().
02319 {
02320 set<MaternalPaternalAlleles> momWorkSet, momSaveSet, dadWorkSet, dadSaveSet;
02321 set<MaternalPaternalAlleles> loopSet;
02322 vector<set<MaternalPaternalAlleles> > offspringWorkVector, offspringSaveVector;
02323 vector<set<MaternalPaternalAlleles> > intersectionVector;
02324 bool sampled = mymother->genotNodeVector[locus].sampled;
02325 unsigned genotypeState = mymother->genotNodeVector[locus].genotypeState;
02326 if(!sampled){
02327 copy(mymother->genotNodeVector[locus].genotypeVector.begin(),
02328 mymother->genotNodeVector[locus].genotypeVector.end(),
02329 inserter(momWorkSet,momWorkSet.begin()) );
02330 }
02331 else {
02332 momWorkSet.insert(mymother->genotNodeVector[locus].genotypeVector[genotypeState]);
02333 mymother->genotNodeVector[locus].genotypeState = 0;
02334 }
02335 sampled = myfather->genotNodeVector[locus].sampled;
02336 genotypeState = myfather->genotNodeVector[locus].genotypeState;
02337 if(!sampled){
02338 copy(myfather->genotNodeVector[locus].genotypeVector.begin(),
02339 myfather->genotNodeVector[locus].genotypeVector.end(),
02340 inserter(dadWorkSet,dadWorkSet.begin()) );
02341 }
02342 else{
02343 dadWorkSet.insert(myfather->genotNodeVector[locus].genotypeVector[genotypeState]);
02344 myfather->genotNodeVector[locus].genotypeState = 0;
02345 }
02346 offspringWorkVector.resize(numoffs);
02347 offspringSaveVector.resize(numoffs);
02348 intersectionVector.resize(numoffs);
02349 for (unsigned i=0;i<numoffs;i++){
02350 Individual *offsPointer = myoffspring[i];
02351 sampled = offsPointer->genotNodeVector[locus].sampled;
02352 genotypeState = offsPointer->genotNodeVector[locus].genotypeState;
02353 if(!sampled){
02354 copy(offsPointer->genotNodeVector[locus].genotypeVector.begin(),
02355 offsPointer->genotNodeVector[locus].genotypeVector.end(),
02356 inserter(offspringWorkVector[i], offspringWorkVector[i].begin()) );
02357 }
02358 else {
02359 offspringWorkVector[i].insert(offsPointer->genotNodeVector[locus].genotypeVector[genotypeState]);
02360 offsPointer->genotNodeVector[locus].genotypeState = 0;
02361 }
02362 }
02363 set<MaternalPaternalAlleles>::iterator itMom, itDad;
02364 for (itMom = momWorkSet.begin(); itMom!=momWorkSet.end(); itMom++){
02365 for (itDad = dadWorkSet.begin(); itDad!=dadWorkSet.end(); itDad++){
02366 loopSet.clear();
02367 MaternalPaternalAlleles matPat;
02368 matPat.maternal = itMom->maternal;
02369 matPat.paternal = itDad->maternal;
02370 loopSet.insert(matPat);
02371 matPat.maternal = itMom->maternal;
02372 matPat.paternal = itDad->paternal;
02373 loopSet.insert(matPat);
02374 matPat.maternal = itMom->paternal;
02375 matPat.paternal = itDad->maternal;
02376 loopSet.insert(matPat);
02377 matPat.maternal = itMom->paternal;
02378 matPat.paternal = itDad->paternal;
02379 loopSet.insert(matPat);
02380 bool saveNot = false;
02381 for (unsigned i=0;i<numoffs;i++){
02382 intersectionVector[i].clear();
02383 set_intersection(loopSet.begin(),loopSet.end(),
02384 offspringWorkVector[i].begin(),
02385 offspringWorkVector[i].end(),
02386 inserter(intersectionVector[i],intersectionVector[i].begin()) );
02387 if(intersectionVector[i].size()==0) {
02388 saveNot = true;
02389 break;
02390 }
02391 }
02392 if (saveNot) continue;
02393 momSaveSet.insert(*itMom);
02394 dadSaveSet.insert(*itDad);
02395 for (unsigned i=0;i<numoffs;i++){
02396 copy(intersectionVector[i].begin(), intersectionVector[i].end(),
02397 inserter(offspringSaveVector[i],offspringSaveVector[i].begin()) );
02398 }
02399 }
02400 }
02401 mymother->genotNodeVector[locus].genotypeVector.clear();
02402 if (momSaveSet.size()==0){
02403 std::cout <<"Incompatible genotypes in nuclear family\n";
02404 std::cout <<"Locus: " << locus << endl;
02405 std::cout <<"Father: " << myfather->myid << endl;
02406 std::cout <<"Mother: " << mymother->myid << endl;
02407 throw exception ("NuFamily::eliminateGenotypes: Incompatible genotypes in pedigree");
02408 }
02409 copy(momSaveSet.begin(),momSaveSet.end(), inserter(mymother->genotNodeVector[locus].genotypeVector,
02410 mymother->genotNodeVector[locus].genotypeVector.begin()) );
02411 myfather->genotNodeVector[locus].genotypeVector.clear();
02412 if (dadSaveSet.size()==0){
02413 std::cout <<"Incompatible genotypes in nuclear family\n";
02414 std::cout <<"Locus: " << locus << endl;
02415 std::cout <<"Father: " << myfather->myid << endl;
02416 std::cout <<"Mother: " << mymother->myid << endl;
02417 throw exception ("NuFamily::eliminateGenotypes: Incompatible genotypes in pedigree");
02418 }
02419 copy(dadSaveSet.begin(),dadSaveSet.end(), inserter(myfather->genotNodeVector[locus].genotypeVector,
02420 myfather->genotNodeVector[locus].genotypeVector.begin()) );
02421 for (unsigned i=0;i<numoffs;i++){
02422 Individual *offsPointer = myoffspring[i];
02423 offsPointer->genotNodeVector[locus].genotypeVector.clear();
02424 if (offspringSaveVector[i].size()==0){
02425 std::cout <<"Incompatible genotypes in nuclear family\n";
02426 std::cout <<"Locus: " << locus << endl;
02427 std::cout <<"Father: " << myfather->myid << endl;
02428 std::cout <<"Mother: " << mymother->myid << endl;
02429 std::cout <<"Offspring: " << offsPointer->myid << endl;
02430 throw exception ("NuFamily::eliminateGenotypes: Incompatible genotypes in pedigree");
02431 }
02432 copy(offspringSaveVector[i].begin(),offspringSaveVector[i].end(),
02433 inserter(offsPointer->genotNodeVector[locus].genotypeVector,
02434 offsPointer->genotNodeVector[locus].genotypeVector.begin()) );
02435 }
02436 }
|
|
|
Definition at line 120 of file nufamily.h.
00120 {return myfather;}
|
|
|
Definition at line 96 of file nufamily.h. Referenced by matvec::Population::build_nufamily(), matvec::Population::maxant_maxpost(), and matvec::Population::maxant_maxpost_old().
00096 { myfather = sire;}
|
|
|
Definition at line 99 of file nufamily.cpp. References matvec::Individual::group_id, and myfather. Referenced by build_connectors().
00100 {
00101 if (myfather) {
00102 return myfather->group_id;
00103 }
00104 else {
00105 return 0;
00106 }
00107 }
|
|
|
Definition at line 80 of file nufamily.cpp. References matvec::Individual::id(), and myfather. Referenced by display().
00081 {
00082 if (myfather) {
00083 return myfather->id();
00084 }
00085 else {
00086 return 0;
00087 }
00088 }
|
|
||||||||||||||||
|
Definition at line 192 of file nufamily.cpp. References matvec::Matrix< double >::assign(), matvec::Individual::get_penetrance(), get_posterior(), myoffspring, matvec::Individual::nspouse(), numoffs, and tn_genotype. Referenced by anterior(), and posterior().
00193 {
00194 if (!trans) throw exception("NuFamily::fullsibs_prob(arg1,arg2,arg3): null arg2");
00195 WSpace.assign(1.0);
00196 unsigned j,um,uf,uj,nsp;
00197 double val, scale;
00198 Vector<double> pen_vec(tn_genotype);
00199 Vector<double> post_vec(tn_genotype);
00200 Individual *J;
00201
00202 for (scale = 0.0, j=0; j<numoffs; j++) {
00203 J = myoffspring[j];
00204 if (J != excludeI) {
00205 nsp = J->nspouse();
00206 if (nsp) scale += get_posterior(J,-1,post_vec);
00207 scale += J->get_penetrance(pen_vec);
00208 for (um=0; um<tn_genotype; um++) {
00209 for (uf=0; uf<tn_genotype; uf++) {
00210 val = 0.0;
00211 if (nsp) {
00212 for (uj=0; uj<tn_genotype; uj++) {
00213 val += trans[uj][um][uf]*pen_vec[uj]*post_vec[uj];
00214 }
00215 }
00216 else {
00217 for (uj=0; uj<tn_genotype; uj++) {
00218 val += trans[uj][um][uf]*pen_vec[uj];
00219 }
00220 }
00221 WSpace[um][uf] *= val;
00222 }
00223 }
00224 }
00225 }
00226 // std::cout << "scale from fullsibs " << scale << std::endl;
00227 return scale;
00228 }
|
|
||||||||||||||||
|
Definition at line 919 of file nufamily.cpp. References matvec::Individual::m_posterior, matvec::Individual::m_posterior_scale, matvec::Individual::n_switches(), matvec::Individual::nspouse(), matvec::Population::P_ndim, and pop. Referenced by multi_anterior(), multi_fullsibs_prob(), multi_llh(), and multi_posterior().
00920 {
00921
00922 /////////////////////////////////////////////////////////
00923 // if excJ < 0, say -1, means taking all posteriors
00924 /////////////////////////////////////////////////////////
00925
00926 // make sure post_mat is initialized properly !!!!!!!!!!!!!!!!
00927
00928 unsigned n_switches = I->n_switches();
00929 int spouse,i,j,k,nsp = I->nspouse();
00930 double retval=0.0;
00931
00932 for (spouse=0; spouse<nsp; spouse++) {
00933 if (spouse != excJ) {
00934 for (k=0;k<pop->P_ndim;k++){
00935 for (i=0;i<n_switches; i++){
00936 for (j=0;j<4;j++){
00937 post_mat[k][i][j] *= (I->m_posterior[spouse])[k][i][j];
00938 }
00939 }
00940 }
00941 // std::cout << "spouse " << spouse << " scale " << (I->m_posterior_scale[spouse]) << std::endl;
00942 retval += I->m_posterior_scale[spouse];
00943 }
00944 // std::cout << I->id() << " return value for spouse " << spouse << " is " << retval << std::endl;
00945
00946 }
00947 return retval;
00948 }
|
|
||||||||||||||||
|
Definition at line 173 of file nufamily.cpp. References matvec::Vector< double >::begin(), matvec::Individual::posterior, and tn_genotype. Referenced by anterior(), fullsibs_prob(), log_likelihood(), and posterior().
00174 {
00175 /////////////////////////////////////////////////////////
00176 // if excJ < 0, say -1, means taking all posteriors
00177 /////////////////////////////////////////////////////////
00178 int i,j,nsp = I->nspouse();
00179 double retval,*ve;
00180 for (i=0; i<tn_genotype; i++) vec[i] = 1.0;
00181 for (retval=0.0,i=0; i<nsp; i++) {
00182 if (i != excJ) {
00183 ve = I->posterior[i].begin();
00184 for (j=0; j<tn_genotype; j++) vec[j] *= *ve++;
00185 retval += *ve;
00186 }
00187 }
00188 // std::cout << "retval from get posterior " << retval << std::endl;
00189 return retval;
00190 }
|
|
|
Definition at line 541 of file nufamily.cpp. References matvec::Population::anterior_posterior(), myfather, myoffspring, numoffs, and pop.
00542 {
00543 //////////////////////////////////////////////////////////////
00544 // REQUIREMENTS: initilization for anterior and posterior
00545 //////////////////////////////////////////////////////////////
00546 pop->anterior_posterior(mymother,WSpace);
00547 pop->anterior_posterior(myfather,WSpace);
00548 for (unsigned j=0; j<numoffs; j++) {
00549 pop->anterior_posterior(myoffspring[j],WSpace);
00550 }
00551 }
|
|
||||||||||||
|
Definition at line 485 of file nufamily.cpp. References matvec::Individual::anterior, get_posterior(), matvec::Vector< double >::inner_product(), kernal, myfather, mymother, posterior(), and tn_genotype.
00486 {
00487 if (!kernal) throw exception(" NuFamily::log_likelihood(): must be kernal nufamily");
00488 Vector<double> post_vec(tn_genotype);
00489 //Note the change in order!!!!!!!
00490 double llh = mymother->anterior[tn_genotype];
00491 posterior(mymother,myfather,trans,WSpace);
00492 llh += get_posterior(mymother,-1,post_vec);
00493 double scale = mymother->anterior.inner_product(post_vec);
00494 llh += std::log(scale);
00495 return llh;
00496 }
|
|
|
Definition at line 121 of file nufamily.h.
00121 {return mymother;}
|
|
|
Definition at line 97 of file nufamily.h. Referenced by matvec::Population::build_nufamily(), matvec::Population::maxant_maxpost(), and matvec::Population::maxant_maxpost_old().
00097 { mymother = dam;}
|
|
|
Definition at line 109 of file nufamily.cpp. References matvec::Individual::group_id, and mymother. Referenced by build_connectors().
00110 {
00111 if (mymother) {
00112 return mymother->group_id;
00113 }
00114 else {
00115 return 0;
00116 }
00117 }
|
|
|
Definition at line 90 of file nufamily.cpp. References matvec::Individual::id(), and mymother. Referenced by display().
00091 {
00092 if (mymother) {
00093 return mymother->id();
00094 }
00095 else {
00096 return 0;
00097 }
00098 }
|
|
||||||||||||
|
Definition at line 1045 of file nufamily.cpp. References matvec::Individual::anterior_iw, matvec::Individual::base(), father_indx, matvec::Individual::id(), matvec::Individual::loop_connector, mother_indx, multi_anterior(), multi_posterior(), myfather, mymother, myoffspring, numoffs, pop, matvec::Population::popmember, matvec::Individual::posterior_iw, and matvec::Population::size(). Referenced by matvec::Population::multi_geno_dist_peeling(), and matvec::Population::multipoint_likelihood().
01045 {
01046
01047 // mother_indx or indx[0] is the index of the mother in the father's spouse list
01048 // father_indx or indx[1] is the index of the father in the mother's spouse list
01049
01050 int spouse_index,i,i_peel=1, tag;
01051 int popsize=pop->size();
01052 Individual* child;
01053
01054 // calculate posterior for father through mother
01055
01056 spouse_index = mother_indx;
01057 if ((!(myfather->base()) && (myfather->loop_connector)) || (myfather->id() > popsize)) {
01058 if (!(myfather->posterior_iw[spouse_index] == 2 ||
01059 myfather->posterior_iw[spouse_index] == 3)) {
01060 // std::cout << "Extending myfather " << myfather->id() << std::endl;
01061 multi_posterior(myfather ,mymother, post_mat_f,i_peel);
01062 }
01063 }
01064
01065 // calculate posterior for mother through father
01066
01067 spouse_index = father_indx;
01068 if ((!(mymother->base()) && (mymother->loop_connector)) || (mymother->id() > popsize)) {
01069 if (!(mymother->posterior_iw[spouse_index] == 2 ||
01070 mymother->posterior_iw[spouse_index] == 3)) {
01071 multi_posterior(mymother ,myfather, post_mat_f,i_peel);
01072 }
01073 }
01074
01075 // calculate anteriors for children
01076
01077 for (i=0;i<numoffs; i++){
01078 child = myoffspring[i];
01079 tag=child->id();
01080 // std::cout << "Extending child original " << tag << " ";
01081 if (tag > popsize) {
01082 child=pop->popmember[tag-popsize-1];
01083 tag -= popsize;
01084 }
01085 // std::cout << "Extending child " << tag << std::endl;
01086 if (child->loop_connector) {
01087 if (!(child->anterior_iw == 2 || child->anterior_iw == 3)) {
01088 multi_anterior(child, post_mat_m, post_mat_f, i_peel);
01089 }
01090 }
01091 }
01092 }
|
|
||||||||||||||||||||
|
Definition at line 739 of file nufamily.cpp. References matvec::Individual::bet_sw, matvec::Individual::eps_sw, matvec::Population::F, father_indx, matvec::Individual::get_m_penetrance(), get_m_posterior(), matvec::Fpmm::getpr(), matvec::Individual::id(), matvec::Individual::index_sw, matvec::Individual::m_anterior, matvec::Individual::m_anterior_scale, mother_indx, multi_fullsibs_prob(), myfather, mymother, matvec::Individual::n_switches(), matvec::Individual::nspouse(), numoffs, matvec::Population::P_ndim, penetrance, pop, matvec::Population::popmember, matvec::Population::size(), matvec::Population::switch_table, matvec::Population::switch_table_gmt, matvec::Population::switch_table_prob, and wspace. Referenced by multi_ant_post(), and multi_terminal_peeling().
00740 {
00741
00742 unsigned i_q,i_sw,f_q,f_sw,m_q,m_sw,w_row,w_col,ndim,i_p,m_p,f_p;
00743 unsigned m_switches = mymother->n_switches();
00744 unsigned f_switches = myfather->n_switches();
00745 double val, scale=0, *ws, i_sum=0,prob_beta, prob_epsilon , t_val, Prob_PGN;
00746 int EQTable[4][4]={2,3,2,3,0,1,0,1,1,0,1,0,3,2,3,2};
00747 int BQTable[4][4]={2,2,3,3,0,0,1,1,1,1,0,0,3,3,2,2};
00748 int Findex_sw=myfather->index_sw;
00749 int Mindex_sw=mymother->index_sw;
00750 int Iindex_sw=I->index_sw;
00751 int eps=I->eps_sw;
00752 int bet=I->bet_sw;
00753 int Bqtl, beta_gamete, Eqtl,epsilon_gamete,Iswitch,Mswitch,Fswitch,i,j,tag_id,popsize;
00754 Individual *temp_mother, *temp_father, *temp_I;
00755 temp_father=myfather;
00756 temp_mother=mymother;
00757 temp_I=I;
00758 popsize=pop->size();
00759 ndim=pop->P_ndim;
00760
00761 if (i_peel) {
00762 //Iterative peeling so have to ensure that we are dealing with the
00763 //original offspring not the cut one.
00764 // check myfather
00765 tag_id=myfather->id();
00766 if (tag_id > popsize) {
00767 // this is a cut individual so need to get the original one
00768 tag_id -= popsize; // correct id;
00769 temp_father=pop->popmember[tag_id-1];
00770 }
00771 // check mymother
00772 tag_id=mymother->id();
00773 if (tag_id > popsize) {
00774 // this is a cut individual so need to get the original one
00775 tag_id -= popsize; // correct id;
00776 temp_mother=pop->popmember[tag_id-1];
00777 }
00778 // check I
00779 tag_id=I->id();
00780 if (tag_id > popsize) {
00781 // this is a cut individual so need to get the original one
00782 tag_id -= popsize; // correct id;
00783 temp_I=pop->popmember[tag_id-1];
00784 }
00785 }
00786 unsigned i_switches = temp_I->n_switches(); // ensure correct entry to switch table
00787
00788 // we are computing fullsibs_prob matrix with rows for mother and columns for father
00789 if ( numoffs > 1 )
00790 scale = multi_fullsibs_prob(I, post_mat_m, i_peel);
00791
00792 // std::cout << "scale after multi_fullsibs_prob " << scale << std::endl;
00793
00794 // don't send temp_I because it is not a child of this nu_family
00795 for (m_p=0; m_p < ndim; m_p++){
00796 for (m_sw=0;m_sw<m_switches; m_sw++){
00797 for (m_q=0;m_q<4;m_q++){
00798 post_mat_m[m_p][m_sw][m_q] = temp_mother->m_anterior[m_p][m_sw][m_q];
00799 }
00800 }
00801 }
00802 scale += temp_mother->m_anterior_scale;
00803
00804 if (temp_mother->nspouse() > 1) {
00805 scale += get_m_posterior(temp_mother,father_indx,post_mat_m);
00806 }
00807
00808 for (f_p=0; f_p < ndim; f_p++){
00809 for (f_sw=0;f_sw<f_switches; f_sw++){
00810 for (f_q=0;f_q<4;f_q++){
00811 post_mat_f[f_p][f_sw][f_q] = temp_father->m_anterior[f_p][f_sw][f_q];
00812 }
00813 }
00814 }
00815 scale += temp_father->m_anterior_scale;
00816
00817 if (temp_father->nspouse() > 1) {
00818 scale += get_m_posterior(temp_father,mother_indx,post_mat_f);
00819 }
00820
00821 scale += I->get_m_penetrance(penetrance);
00822
00823 if ( numoffs > 1 ) {
00824 for (i_p=0; i_p < ndim; i_p++){ // Individuals PGN
00825 for (i_q=0; i_q<4; i_q++) {
00826 for (i_sw=0;i_sw<i_switches; i_sw++) { // loop for genotypes of ind
00827 w_row = 0;
00828 Prob_PGN = 0.0;
00829 Iswitch=pop->switch_table[Iindex_sw][i_sw+1]; // switch for ind (i_sw)
00830 for (m_p=0; m_p < ndim; m_p++){ // loop for PGN of mother
00831 for (m_q=0;m_q<4;m_q++){
00832 Bqtl=BQTable[m_q][i_q]; // Find which type QTL allele from dam
00833 // if (Bqtl !=3) { // need to be compatible
00834 for (m_sw=0;m_sw<m_switches; m_sw++){ // loop for genotypes of mother
00835 Mswitch=pop->switch_table[Mindex_sw][m_sw+1]; // switch for mother (m_sw)
00836 beta_gamete=pop->switch_table_gmt[bet][(Mswitch^Iswitch)]; // get coded mat marker gamete
00837 prob_beta=pop->switch_table_prob[beta_gamete][Bqtl];// Get Pr of mat gamete with qtl
00838 w_col = 0;
00839 for (f_p=0; f_p < ndim; f_p++){ // loop for PGN of father
00840 val = 0.0;
00841 for (f_q=0;f_q<4;f_q++){
00842 Eqtl=EQTable[f_q][i_q]; // Find which type QTL allele from sire
00843 // if (Eqtl != 3) { // need to be compatible
00844 for (f_sw=0;f_sw<f_switches; f_sw++){ // loop for genotypes of father
00845 Fswitch=pop->switch_table[Findex_sw][f_sw+1]; // switch for father (f_sw)
00846 epsilon_gamete=pop->switch_table_gmt[eps][(Fswitch^Iswitch)]; // get coded paternal marker gamete
00847 prob_epsilon=pop->switch_table_prob[epsilon_gamete][Eqtl];// Get Pr of paternal gamete
00848 if ( (Eqtl != 3) && (Bqtl !=3)) {
00849 val += post_mat_m[m_p][m_sw][m_q]*wspace[w_row][w_col]*post_mat_f[f_p][f_sw][f_q] *prob_beta*prob_epsilon;
00850 }
00851 w_col++; // increment w_col only
00852 }
00853 }
00854 Prob_PGN += val*pop->F->getpr(f_p,m_p,i_p);
00855 }
00856 w_row++; // increment w_row only
00857 }
00858 }
00859 }
00860 t_val=Prob_PGN*penetrance[i_p][i_q];
00861 temp_I->m_anterior[i_p][i_sw][i_q] = t_val; // store value
00862 i_sum += t_val; // accumulate sum of values for scaling
00863 }
00864 }
00865 }
00866 }
00867 else {
00868 for (i_p=0; i_p < ndim; i_p++){ // Individuals PGN
00869 for (i_q=0; i_q<4; i_q++) {
00870 for (i_sw=0;i_sw<i_switches; i_sw++) { // loop for genotypes of ind
00871 Prob_PGN=0.0;
00872 Iswitch=pop->switch_table[Iindex_sw][i_sw+1]; // switch for ind (i_sw)
00873 for (m_p=0; m_p < ndim; m_p++){ // loop for PGN of mother
00874 for (m_q=0;m_q<4;m_q++){
00875 Bqtl=BQTable[m_q][i_q]; // Find which type QTL allele from dam
00876 if (Bqtl !=3) { // need to be compatible
00877 for (m_sw=0;m_sw<m_switches; m_sw++){ // loop for genotypes of mother
00878 Mswitch=pop->switch_table[Mindex_sw][m_sw+1]; // switch for mother (m_sw)
00879 beta_gamete=pop->switch_table_gmt[bet][(Mswitch^Iswitch)]; // get coded maternal marker gamete
00880 prob_beta=pop->switch_table_prob[beta_gamete][Bqtl]; // Get prob of maternal gamete with qtl
00881 for (f_p=0; f_p < ndim; f_p++){ // loop for PGN of father
00882 val = 0.0;
00883 for (f_q=0;f_q<4;f_q++){
00884 Eqtl=EQTable[f_q][i_q]; // Find which type QTL allele from sire
00885 if (Eqtl != 3) { // need to be compatible
00886 for (f_sw=0;f_sw<f_switches; f_sw++){ // loop for genotypes of father
00887 Fswitch=pop->switch_table[Findex_sw][f_sw+1]; // switch for father (f_sw)
00888 epsilon_gamete=pop->switch_table_gmt[eps][(Fswitch^Iswitch)]; // get coded paternal marker gamete
00889 prob_epsilon=pop->switch_table_prob[epsilon_gamete][Eqtl]; // Get prob of paternal gamete
00890 val += post_mat_m[m_p][m_sw][m_q]*post_mat_f[f_p][f_sw][f_q] *prob_beta*prob_epsilon;
00891 }
00892 }
00893 }
00894 Prob_PGN += val*pop->F->getpr(f_p,m_p,i_p);
00895 }
00896 }
00897 }
00898 }
00899 }
00900 t_val=Prob_PGN*penetrance[i_p][i_q];
00901 temp_I->m_anterior[i_p][i_sw][i_q] = t_val; // store value
00902 i_sum += t_val; // accumulate sum of values for scaling
00903 }
00904 }
00905 }
00906 }
00907 for (i_p=0; i_p < ndim; i_p++){ // Individuals PGN
00908 for (i_q=0; i_q<4; i_q++){
00909 for (i_sw=0;i_sw<i_switches; i_sw++) {
00910 temp_I->m_anterior[i_p][i_sw][i_q] /= i_sum;
00911 }
00912 }
00913 }
00914 scale += std::log(i_sum);
00915 temp_I->m_anterior_scale = scale;
00916 // std::cout << " id " << temp_I->id() << " Anterior scale " << scale << std::endl;
00917 }
|
|
||||||||||||||||
|
Definition at line 950 of file nufamily.cpp. References matvec::Matrix< double >::assign(), matvec::Individual::bet_sw, matvec::Individual::eps_sw, matvec::Population::F, matvec::Individual::get_m_penetrance(), get_m_posterior(), matvec::Fpmm::getpr(), matvec::Individual::id(), matvec::Individual::index_sw, matvec::Dblock::init(), myfather, mymother, myoffspring, matvec::Individual::n_switches(), matvec::Individual::nspouse(), numoffs, matvec::Population::P_ndim, penetrance, pop, matvec::Population::popmember, matvec::Population::size(), matvec::Population::switch_table, matvec::Population::switch_table_gmt, matvec::Population::switch_table_prob, and wspace. Referenced by multi_anterior(), and multi_posterior().
00950 {
00951 // fullsib probs with rows for mother's genotypes and columns for father's
00952
00953 wspace.assign(1.0);
00954 unsigned nsp;
00955 int j, eps, bet, Jswitch, j_switches, w_row, w_col, m_q, m_sw, Fswitch, Mswitch,
00956 f_q, f_sw, Bqtl, Eqtl, j_q, j_sw, Jindex_sw, beta_gamete, epsilon_gamete, m_p, f_p, j_p;
00957 double scale=0.0, ProbIP, ProbIQ, ProbIS, prob_beta, prob_epsilon;
00958 Individual *J;
00959 int EQTable[4][4]={2,3,2,3,0,1,0,1,1,0,1,0,3,2,3,2};
00960 int BQTable[4][4]={2,2,3,3,0,0,1,1,1,1,0,0,3,3,2,2};
00961 int tag_id;
00962 int popsize=pop->size();
00963 int m_switches = mymother->n_switches();
00964 int f_switches = myfather->n_switches();
00965 int Findex_sw=myfather->index_sw;
00966 int Mindex_sw=mymother->index_sw;
00967 int ndim = pop->P_ndim;
00968
00969 for (j=0; j<numoffs; j++) {
00970 J = myoffspring[j];
00971
00972 if (i_peel) {
00973 //Iterative peeling so have to ensure that we are dealing with the
00974 //original offspring not the cut one.
00975 tag_id=J->id();
00976 if (tag_id > popsize) {
00977 // this is a cut individual so need to get the original one
00978 tag_id -= popsize; // correct id;
00979 J=pop->popmember[tag_id-1];
00980 }
00981 }
00982 if (J != excludeI) {
00983 nsp = J->nspouse();
00984 if (nsp) {
00985 post_mat.init(1.0);
00986 scale += get_m_posterior(J,-1,post_mat);
00987 }
00988
00989 scale += J->get_m_penetrance(penetrance);
00990 eps=J->eps_sw;
00991 bet=J->bet_sw;
00992 Jindex_sw=J->index_sw;
00993 j_switches = J->n_switches();
00994 w_row=0; // wspace row index
00995 for (m_p=0; m_p < ndim; m_p++) { // for each mother PGN
00996 for (m_q=0; m_q < 4; m_q++) { // for each dam qtl
00997 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of mother
00998 Mswitch=pop->switch_table[Mindex_sw][m_sw+1]; // get current dam switch
00999 w_col=0; // wspace column index
01000 for (f_p=0; f_p < ndim; f_p++) { //For each sire PGN
01001 for (f_q=0; f_q < 4; f_q++) { //For each sire QTL
01002 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of father
01003 Fswitch=pop->switch_table[Findex_sw][f_sw+1]; // switch for father (f_sw)
01004 ProbIP=0.0;
01005 for (j_p=0; j_p < ndim; j_p++) { // for each child PGN
01006 ProbIQ=0.0;
01007 for (j_q=0; j_q < 4; j_q++) { // for each child QTL
01008 Bqtl=BQTable[m_q][j_q]; // Find which type QTL allele from dam
01009 Eqtl=EQTable[f_q][j_q]; // Find which type QTL allele from sire
01010 if ((Eqtl != 3) && (Bqtl !=3)) { // need to be compatible
01011 ProbIS=0.0;
01012 for (j_sw=0; j_sw < j_switches; j_sw++) { // for each child switch
01013 Jswitch=pop->switch_table[Jindex_sw][j_sw+1];
01014 beta_gamete=pop->switch_table_gmt[bet][(Mswitch^Jswitch)]; // get coded maternal marker gamete
01015 prob_beta=pop->switch_table_prob[beta_gamete][Bqtl]; // Get prob of maternal gamete
01016 epsilon_gamete=pop->switch_table_gmt[eps][(Fswitch^Jswitch)]; // get coded paternal marker gamete
01017 prob_epsilon=pop->switch_table_prob[epsilon_gamete][Eqtl]; // Get prob of paternal gamete
01018 if (nsp) {
01019 ProbIS += (prob_epsilon*prob_beta*post_mat[j_p][j_sw][j_q]);
01020 }
01021 else {
01022 ProbIS += (prob_epsilon*prob_beta);
01023 }
01024 }
01025 ProbIQ += (ProbIS*penetrance[j_p][j_q]);
01026 }
01027 }
01028 ProbIP += ProbIQ*pop->F->getpr(f_p,m_p,j_p);
01029 }
01030 wspace[w_row][w_col] *= ProbIP;
01031 w_col++;
01032 }
01033 }
01034 }
01035 w_row++;
01036 }
01037 }
01038 }
01039 }
01040 }
01041 // std::cout << "scale from multi_fullsib " << scale << std::endl;
01042 return scale;
01043 }
|
|
||||||||||||||||||||||||
|
Definition at line 1108 of file nufamily.cpp. References matvec::Individual::bet_sw, matvec::Individual::eps_sw, matvec::Individual::index_sw, matvec::Individual::n_switches(), pop, matvec::Population::switch_table, matvec::Population::switch_table_gmt, matvec::Population::switch_table_prob, and tr. Referenced by multi_sumint_offspring().
01108 {
01109
01110 int j_q,j_sw,j_switches,Jswitch,beta_gamete,epsilon_gamete,bet,eps,Eqtl,Bqtl;
01111 double prob_beta,prob_epsilon;
01112 int EQTable[4][4]={2,3,2,3,0,1,0,1,1,0,1,0,3,2,3,2};
01113 int BQTable[4][4]={2,2,3,3,0,0,1,1,1,1,0,0,3,3,2,2};
01114
01115 eps=J->eps_sw;
01116 bet=J->bet_sw;
01117 int Jindex_sw=J->index_sw;
01118 j_switches = J->n_switches();
01119 for (j_q=0; j_q < 4; j_q++) { // for each child QTL
01120 Bqtl=BQTable[m_q][j_q]; // Find which type QTL allele from dam
01121 Eqtl=EQTable[f_q][j_q]; // Find which type QTL allele from sire
01122 if ((Eqtl != 3) && (Bqtl !=3)) { // need to be compatible
01123 for (j_sw=0; j_sw < j_switches; j_sw++) { // for each child switch
01124 Jswitch=pop->switch_table[Jindex_sw][j_sw+1];
01125 beta_gamete=pop->switch_table_gmt[bet][(Mswitch^Jswitch)]; // get coded maternal marker gamete
01126 prob_beta=pop->switch_table_prob[beta_gamete][Bqtl]; // Get prob of maternal gamete
01127 epsilon_gamete=pop->switch_table_gmt[eps][(Fswitch^Jswitch)]; // get coded paternal marker gamete
01128 prob_epsilon=pop->switch_table_prob[epsilon_gamete][Eqtl]; // Get prob of paternal gamete
01129 tr[j_q][j_sw] = prob_beta*prob_epsilon;
01130 }
01131 }
01132 else {
01133 for (j_sw=0; j_sw < j_switches; j_sw++) {
01134 tr[j_q][j_sw] = 0.0;
01135 }
01136 }
01137 }
01138 }
|
|
|
Definition at line 2258 of file nufamily.cpp. References matvec::Individual::anterior_iw, father_indx, matvec::Individual::initial_multi_anterior(), matvec::Individual::initial_multi_posterior(), mother_indx, myfather, mymother, myoffspring, numoffs, and matvec::Individual::posterior_iw. Referenced by matvec::Population::multipoint_likelihood().
02258 {
02259 // Initialize ALL members of the Nuclear family.
02260 // Else values are wrong for profile and maximization steps!
02261 // Would like to do just founders but cut individuals do not have this flag set.
02262 // WARNING I have reset flags may not be correct
02263 int nchild;
02264 Individual *I;
02265
02266 // if (myfather->base()) {
02267
02268 myfather->initial_multi_anterior(pen);
02269 myfather->initial_multi_posterior(-1);
02270 myfather->posterior_iw[mother_indx] = 0;
02271 // }
02272
02273 // if (mymother->base()) {
02274
02275 mymother->initial_multi_anterior(pen);
02276 mymother->initial_multi_posterior(-1);
02277 mymother->posterior_iw[father_indx] = 0;
02278 // }
02279
02280 for (nchild=0; nchild < numoffs; nchild++) {
02281 I = myoffspring[nchild];
02282 I->initial_multi_anterior(pen);
02283 // I->initial_multi_posterior(-1); //BRS Child doesn't need posterior unless parent?
02284 I->anterior_iw = 0;
02285 // std::cout << I->id() << " initial ant scale " << I->m_anterior_scale << std::endl;
02286 }
02287 }
|
|
|
Definition at line 585 of file nufamily.cpp. References get_m_posterior(), matvec::Dblock::init(), kernal, matvec::Individual::m_anterior, matvec::Individual::m_anterior_scale, multi_posterior(), myfather, mymother, matvec::Individual::n_switches(), matvec::Population::P_ndim, and pop. Referenced by matvec::Population::multipoint_likelihood().
00585 {
00586 unsigned i,j,k;
00587 double scale;
00588
00589 // std::cout << "entering multi_llh " << std::endl;
00590 if (!kernal) throw exception(" NuFamily::log_likelihood(): must be kernal nufamily");
00591 double llh = myfather->m_anterior_scale;
00592 // std::cout << "LLh " << llh << " " << (myfather->m_anterior) << std::endl;
00593 // std::cout << "multi_llh m_anterior_scale " << llh << std::endl;
00594 multi_posterior(myfather,mymother,post_mat);
00595 post_mat.init(1.0);
00596 llh += get_m_posterior(myfather,-1,post_mat);
00597 // std::cout << "llh from get m posterior " << llh << std::endl;
00598 scale = 0.0;
00599 for (k=0; k<pop->P_ndim; k++) {
00600 for (i=0;i< myfather->n_switches(); i++) {
00601 for (j=0;j<4; j++) {
00602 scale += ((post_mat[k][i][j])*(myfather->m_anterior[k][i][j]));
00603 }
00604 }
00605 }
00606 // std::cout << "Scale " << scale << std::endl;
00607 llh += std::log(scale);
00608 // std::cout << "exiting multi_llh " << std::endl;
00609 return llh;
00610 }
|
|
|
Definition at line 2210 of file nufamily.cpp. References matvec::Individual::anterior_iw, matvec::Individual::base(), father_indx, matvec::Individual::id(), matvec::Individual::loop_connector, mother_indx, multi_m_anterior(), multi_m_posterior(), myfather, mymother, myoffspring, numoffs, pop, matvec::Population::popmember, matvec::Individual::posterior_iw, and matvec::Population::size(). Referenced by matvec::Population::multi_m_geno_dist_peeling(), and matvec::Population::multi_m_log_likelihood_peeling().
02210 {
02211
02212 // indx[0] is the index of the mother in the father's spouse list
02213 // indx[1] is the index of the father in the mother's spouse list
02214
02215 int spouse_index,i,i_peel=1, tag;
02216 Individual* child;
02217 int popsize=pop->size();
02218
02219 // calculate posterior for father through mother
02220
02221 spouse_index = mother_indx;
02222 if ((!(myfather->base()) && (myfather->loop_connector)) || (myfather->id() > popsize)) {
02223 if (!(myfather->posterior_iw[spouse_index] == 2 ||
02224 myfather->posterior_iw[spouse_index] == 3)) {
02225 multi_m_posterior(myfather ,mymother, i_peel);
02226 }
02227 }
02228
02229 // calculate posterior for mother through father
02230
02231 spouse_index = father_indx;
02232 if ((!(mymother->base()) && (mymother->loop_connector)) || (mymother->id() > popsize)) {
02233 if (!(mymother->posterior_iw[spouse_index] == 2 ||
02234 mymother->posterior_iw[spouse_index] == 3)) {
02235 multi_m_posterior(mymother ,myfather, i_peel);
02236 }
02237 }
02238
02239 // calculate anteriors for children
02240
02241 for (i=0;i<numoffs; i++){
02242 child = myoffspring[i];
02243 tag=child->id();
02244 // std::cout << "Extending child original " << tag << " ";
02245 if (tag > popsize) {
02246 child=pop->popmember[tag-popsize-1];
02247 tag -= popsize;
02248 }
02249 // std::cout << "Extending child " << tag << std::endl;
02250 if (child->loop_connector) {
02251 if (!(child->anterior_iw == 2 || child->anterior_iw == 3)) {
02252 multi_m_anterior(child, i_peel);
02253 }
02254 }
02255 }
02256 }
|
|
||||||||||||
|
Definition at line 1659 of file nufamily.cpp. References matvec::Individual::bet_sw, child_matrix, matvec::DataNode::double_val(), matvec::Individual::eps_sw, matvec::Population::F, father_indx, matvec::Individual::id(), matvec::Individual::index_sw, matvec::Vector< Vector< Sym2x2 > >::initialize(), m_weight, matvec::Population::mean_for_genotype, matvec::DataNode::missing, matvec::Individual::mix_anterior, matvec::Individual::mix_posterior, matvec::MIXED_TOL, mother_indx, multi_sumint_offspring(), myfather, myfather_matrix, mymother, mymother_matrix, myoffspring, matvec::Individual::n_switches(), matvec::Individual::nspouse(), numoffs, offspring(), pm, pop, matvec::Population::popmember, matvec::Individual::record(), matvec::Population::residual_var, matvec::Population::size(), matvec::Population::switch_table, matvec::Population::switch_table_gmt, matvec::Population::switch_table_prob, tn_qtl, and matvec::Fpmm::var. Referenced by multi_m_ant_post(), and multi_m_terminal_peeling().
01659 {
01660
01661 double a11,a12,a13,a14,a22,a23,a24,a33,a34,a44,k;
01662 double u11, u12, u22, sum1, sum2, v, nu, y, v_y, k_y, nu_y;
01663 double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0)); // pi = 2.0*std::asin(1.0)
01664 unsigned i_q,m_q,f_q,j,independent=0;
01665 int m_sw, f_sw, i_sw, offspring;
01666 int m_switches = mymother->n_switches();
01667 int f_switches = myfather->n_switches();
01668 int Mindex_sw = mymother->index_sw;
01669 int Findex_sw = myfather->index_sw;
01670 int EQTable[4][4]={2,3,2,3,0,1,0,1,1,0,1,0,3,2,3,2};
01671 int BQTable[4][4]={2,2,3,3,0,0,1,1,1,1,0,0,3,3,2,2};
01672 int eps=I->eps_sw;
01673 int bet=I->bet_sw;
01674 int Iindex_sw=I->index_sw;
01675 int i_switches = I->n_switches();
01676 int wrow=0, wcol=0, Iswitch, Bqtl, Mswitch, beta_gamete, Eqtl, Fswitch, epsilon_gamete;
01677 double prob_beta, prob_epsilon, temp;
01678 Individual *temp_mother, *temp_father, *temp_I;
01679 temp_father=myfather;
01680 temp_mother=mymother;
01681 temp_I=I;
01682 int popsize=pop->size();
01683 int tag_id;
01684 double P_var = pop->F->var; // polygenic varaince
01685 // indx[0] is the index of the mother in the father's spouse list
01686 // indx[1] is the index of the father in the mother's spouse list
01687
01688 unsigned excM = mother_indx;
01689 unsigned excF = father_indx;
01690 double ve=pop->residual_var[0][0];
01691
01692 if (i_peel) {
01693 //Iterative peeling so have to ensure that we are dealing with the
01694 //original individual not the cut one.
01695 // check myfather
01696 tag_id=myfather->id();
01697 if (tag_id > popsize) {
01698 // this is a cut individual so need to get the original one
01699 tag_id -= popsize; // correct id;
01700 temp_father=pop->popmember[tag_id];
01701 }
01702 // check mymother
01703 tag_id=mymother->id();
01704 if (tag_id > popsize) {
01705 // this is a cut individual so need to get the original one
01706 tag_id -= popsize; // correct id;
01707 temp_mother=pop->popmember[tag_id];
01708 }
01709 // check I
01710 tag_id=I->id();
01711 if (tag_id > popsize) {
01712 // this is a cut individual so need to get the original one
01713 tag_id -= popsize; // correct id;
01714 temp_I=pop->popmember[tag_id];
01715 }
01716 }
01717 // std::cout << "Anterior for " << pop->ind_name(temp_I->id()) << std::endl;
01718 // First, integrate and sum over each offspring
01719 // Approximate results are stored for each gm and f_q
01720 // in child_matrix
01721
01722 // Initialize child_matrix for this anterior
01723 for (m_q=0;m_q<tn_qtl;m_q++){
01724 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of I
01725 wcol = 0;
01726 for (f_q=0;f_q<tn_qtl; f_q++){
01727 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of J
01728 child_matrix[wrow][wcol].initialize();
01729 wcol++;
01730 }
01731 }
01732 wrow++;
01733 }
01734 }
01735
01736 // sum and int. each offspring except I
01737 for (offspring=0; offspring<numoffs; offspring++) { // sum and int. each offspring
01738 if ( temp_I != myoffspring[offspring]){
01739 multi_sumint_offspring(offspring);
01740 }
01741 }
01742
01743 // Collect results to integrate uf for each f_q
01744
01745 for (f_q=0;f_q<tn_qtl;f_q++){
01746 // contributions from penetrance function of father
01747 // skip if phenotype of father is missing
01748
01749 if (!(temp_father->record()[0].missing)) {
01750 y = temp_father->record()[0].double_val();
01751 v_y = 1.0/ve;
01752 k_y = std::log(inverse_sqrt2pi*std::sqrt(v_y));
01753 nu_y = pop->mean_for_genotype[f_q] - y;
01754 }
01755 else{
01756 nu_y = 0.0;
01757 v_y = 0.0;
01758 k_y = 0.0;
01759 }
01760
01761 for (f_sw=0; f_sw < f_switches; f_sw++) { // loop for genotypes of J
01762
01763 a11 = nu_y*nu_y*v_y;
01764 a13 = nu_y*v_y;
01765 a33 = v_y;
01766 k = k_y;
01767 // contributions from Father's anterior
01768
01769 v = temp_father->mix_anterior[f_q][f_sw].tsq;
01770 nu = temp_father->mix_anterior[f_q][f_sw].nu;
01771 k += temp_father->mix_anterior[f_q][f_sw].k;
01772 a11 += nu*nu*v;
01773 a13 += -nu*v;
01774 a33 += v;
01775
01776 // Contributions from posteriors of father except mother.
01777 // Skip over any posteriors not yet available
01778
01779 unsigned numspouse = temp_father->nspouse();
01780
01781 for (j=0; j<numspouse; j++){
01782 if ((temp_father->mix_posterior[j].done) && (j!=excM) ) {
01783 v = temp_father->mix_posterior[j].postvec[f_q][f_sw].tsq;
01784 nu = temp_father->mix_posterior[j].postvec[f_q][f_sw].nu;
01785 k += temp_father->mix_posterior[j].postvec[f_q][f_sw].k;
01786 a11 += nu*nu*v;
01787 a13 += -nu*v;
01788 a33 += v;
01789 }
01790 }
01791
01792 // Put results into temp_father matrix
01793
01794 myfather_matrix[f_q][f_sw].a11 = a11;
01795 myfather_matrix[f_q][f_sw].a13 = a13;
01796 myfather_matrix[f_q][f_sw].a33 = a33;
01797 myfather_matrix[f_q][f_sw].k = k;
01798 }
01799 }
01800
01801 // Collect results to integrate um for each m_q
01802
01803 for (m_q=0;m_q<tn_qtl;m_q++){
01804
01805 // contributions from penetrance funtion of mother
01806 // skip if phenotype of mother is missing
01807
01808 if (!(temp_mother->record()[0].missing)) {
01809 y = temp_mother->record()[0].double_val();
01810 v_y = 1.0/ve;
01811 k_y = std::log(inverse_sqrt2pi*std::sqrt(v_y));
01812 nu_y = pop->mean_for_genotype[m_q] - y;
01813 }
01814 else{
01815 nu_y = 0.0;
01816 v_y = 0.0;
01817 k_y = 0.0;
01818 }
01819 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for switches of mother
01820
01821 a11 = nu_y*nu_y*v_y;
01822 a14 = nu_y*v_y;
01823 a44 = v_y;
01824 k = k_y;
01825 // contributions from Mother's anterior
01826
01827 v = temp_mother->mix_anterior[m_q][m_sw].tsq;
01828 nu = temp_mother->mix_anterior[m_q][m_sw].nu;
01829 k += temp_mother->mix_anterior[m_q][m_sw].k;
01830 a11 += nu*nu*v;
01831 a14 += -nu*v;
01832 a44 += v;
01833
01834 // Contributions from posteriors of mother except father.
01835 // Skip over any posteriors not yet available
01836
01837 unsigned numspouse = temp_mother->nspouse();
01838 for (j=0; j<numspouse; j++){
01839 if (( temp_mother->mix_posterior[j].done) && (j!=excF) ) {
01840 v = temp_mother->mix_posterior[j].postvec[m_q][m_sw].tsq;
01841 nu = temp_mother->mix_posterior[j].postvec[m_q][m_sw].nu;
01842 k += temp_mother->mix_posterior[j].postvec[m_q][m_sw].k;
01843 a11 += nu*nu*v;
01844 a14 += -nu*v;
01845 a44 += v;
01846 }
01847 }
01848
01849 // Put results into mother_matrix
01850
01851 mymother_matrix[m_q][m_sw].a11 = a11;
01852 mymother_matrix[m_q][m_sw].a14 = a14;
01853 mymother_matrix[m_q][m_sw].a44 = a44;
01854 mymother_matrix[m_q][m_sw].k = k;
01855 }
01856 }
01857
01858 // For each i_q and j_q complete integration of um and uf
01859 // and store results in pm
01860 sum1=0.0;
01861 wrow=0;
01862 for (m_q=0; m_q<tn_qtl; m_q++) {
01863 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for switches of mother
01864 wcol=0;
01865 for (f_q=0; f_q<tn_qtl; f_q++){
01866 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for switches of father
01867
01868 a11=a12=a13=a14=a22=a23=a24=a33=a34=a44=k=0.0;
01869
01870 // Get contributions from mother_vector for genotype m_q
01871
01872 a11 = mymother_matrix[m_q][m_sw].a11;
01873 a14 = mymother_matrix[m_q][m_sw].a14;
01874 a44 = mymother_matrix[m_q][m_sw].a44;
01875 k = mymother_matrix[m_q][m_sw].k ;
01876
01877 // Get contributions from child_matrix for genotypes f_q and m_q
01878
01879 a11 += child_matrix[wrow][wcol].a11;
01880 a13 += child_matrix[wrow][wcol].a12;
01881 a14 += child_matrix[wrow][wcol].a12;
01882 a33 += child_matrix[wrow][wcol].a22;
01883 a34 += child_matrix[wrow][wcol].a22;
01884 a44 += child_matrix[wrow][wcol].a22;
01885 k += child_matrix[wrow][wcol].k;
01886
01887 // Contributions from transition function of I to integrate um
01888
01889 v = 1.0/(0.5*P_var);
01890 a22 += v;
01891 a23 += -0.5*v;
01892 a24 += -0.5*v;
01893 a33 += 0.25*v;
01894 a34 += 0.25*v;
01895 a44 += 0.25*v;
01896 k += std::log(inverse_sqrt2pi*std::sqrt(v));
01897
01898 // Integrate um
01899
01900 v = 1.0/a44;
01901 a11 = a11 - a14*a14*v;
01902 a12 = a12 - a14*a24*v;
01903 a13 = a13 - a14*a34*v;
01904 a22 = a22 - a24*a24*v;
01905 a23 = a23 - a24*a34*v;
01906 a33 = a33 - a34*a34*v;
01907 k += std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01908
01909 // Get contributions from father
01910
01911 a11 += myfather_matrix[f_q][f_sw].a11;
01912 a13 += myfather_matrix[f_q][f_sw].a13;
01913 a33 += myfather_matrix[f_q][f_sw].a33;
01914 k += myfather_matrix[f_q][f_sw].k ;
01915
01916 // Integrate uf
01917
01918 v = 1.0/a33;
01919 u11 = a11 - a13*a13*v;
01920 u12 = a12 - a13*a23*v;
01921 u22 = a22 - a23*a23*v;
01922 k += std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01923
01924 if(u22 < MIXED_TOL) {
01925 pm[wrow][wcol].k = k - 0.5*u11;
01926 independent = 1;
01927 }
01928 else {
01929
01930 // Calculate Normal mean, variance and k for f_q and m_q
01931
01932 v = 1.0 /u22;
01933 nu = -u12*v;
01934 pm[wrow][wcol].tsq = v;
01935 pm[wrow][wcol].nu = nu;
01936 pm[wrow][wcol].k = k - 0.5*(u11 - nu*nu*u22)
01937 + std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01938 sum1 += pm[wrow][wcol].k;
01939 }
01940 wcol++;
01941 }
01942 }
01943 wrow++;
01944 }
01945 }
01946 sum1 /= (tn_qtl*f_switches*tn_qtl*m_switches);
01947
01948 if (independent) {
01949 // The posterior and anterior likelihoods are independent
01950 // This happens due to missing data
01951
01952 for (i_q=0; i_q<tn_qtl; i_q++) {
01953 for (i_sw=0; i_sw < i_switches; i_sw++){ // loop for genotypes of I
01954 Iswitch=pop->switch_table[Iindex_sw][i_sw+1]; // switch for ind (i_sw)
01955 wrow=0;
01956 sum2=0.0;
01957 for(m_q=0;m_q<tn_qtl; m_q++){
01958 Bqtl=BQTable[m_q][i_q]; // Find which type QTL allele from dam
01959 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of J
01960 Mswitch=pop->switch_table[Mindex_sw][m_sw+1]; // switch for mother (m_sw)
01961 beta_gamete=pop->switch_table_gmt[bet][(Mswitch^Iswitch)]; // get coded maternal marker gamete
01962 prob_beta=pop->switch_table_prob[beta_gamete][Bqtl]; // Get prob of maternal gamete
01963 wcol=0;
01964 for (f_q=0; f_q<tn_qtl; f_q++) {
01965 Eqtl=EQTable[f_q][i_q]; // Find which type QTL allele from sire
01966 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of father
01967 Fswitch=pop->switch_table[Findex_sw][f_sw+1]; // switch for father (f_sw)
01968 epsilon_gamete=pop->switch_table_gmt[eps][(Fswitch^Iswitch)]; // get coded paternal marker gamete
01969 prob_epsilon=pop->switch_table_prob[epsilon_gamete][Eqtl]; // Get prob of paternal gamete
01970 if ( (Eqtl != 3) && (Bqtl !=3)) {
01971 sum2 += prob_beta*prob_epsilon*std::exp(pm[wrow][wcol].k - sum1);
01972 }
01973 wcol++;
01974 }
01975 }
01976 wrow++;
01977 }
01978 }
01979
01980 temp_I->mix_anterior[i_q][i_sw].nu = 0.0;
01981 temp_I->mix_anterior[i_q][i_sw].tsq = 0.0;
01982 temp_I->mix_anterior[i_q][i_sw].k = std::log(sum2) + sum1;
01983 // std::cout << i_q << " " << i_sw << " k " << temp_I->mix_anterior[i_q][i_sw].k << std::endl;
01984 }
01985 }
01986 return;
01987 }
01988
01989 // For each i_q, approximate the mixture that results from summing over
01990 // f_q and m_q by an univariate normal with the same mean and variance
01991 // as the mixture
01992
01993 for (i_q=0; i_q<tn_qtl; i_q++) {
01994 for (i_sw=0; i_sw < i_switches; i_sw++){ // loop for genotypes of I
01995 Iswitch=pop->switch_table[Iindex_sw][i_sw+1]; // switch for ind (i_sw)
01996 wrow=0;
01997
01998 // Weights of mixture are calculated here.
01999 sum2=0.0;
02000 for(m_q=0;m_q<tn_qtl; m_q++){
02001 Bqtl=BQTable[m_q][i_q]; // Find which type QTL allele from dam
02002 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of J
02003 Mswitch=pop->switch_table[Mindex_sw][m_sw+1]; // switch for mother (m_sw)
02004 beta_gamete=pop->switch_table_gmt[bet][(Mswitch^Iswitch)]; // get coded maternal marker gamete
02005 prob_beta=pop->switch_table_prob[beta_gamete][Bqtl]; // Get prob of maternal gamete
02006 wcol=0;
02007 for (f_q=0; f_q<tn_qtl; f_q++) {
02008 Eqtl=EQTable[f_q][i_q]; // Find which type QTL allele from sire
02009 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of father
02010 Fswitch=pop->switch_table[Findex_sw][f_sw+1]; // switch for father (f_sw)
02011 epsilon_gamete=pop->switch_table_gmt[eps][(Fswitch^Iswitch)]; // get coded paternal marker gamete
02012 prob_epsilon=pop->switch_table_prob[epsilon_gamete][Eqtl]; // Get prob of paternal gamete
02013 if ( (Eqtl != 3) && (Bqtl !=3)) {
02014 temp = prob_beta*prob_epsilon*std::exp(pm[wrow][wcol].k - sum1);
02015 sum2 += temp;
02016 m_weight[wrow][wcol]=temp;
02017 }
02018 else {
02019 m_weight[wrow][wcol]=0.0; // because the means are only due the penetrance function
02020 }
02021 wcol++;
02022 }
02023 }
02024 wrow++;
02025 }
02026 }
02027
02028 nu = v = 0.0;
02029 wrow=0;
02030 for (m_q=0; m_q<tn_qtl; m_q++) {
02031 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of mother
02032 wcol=0;
02033 for (f_q=0; f_q<tn_qtl; f_q++){
02034 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of father
02035 m_weight[wrow][wcol] /= sum2;
02036 nu += pm[wrow][wcol].nu*m_weight[wrow][wcol];
02037 wcol++;
02038 }
02039 }
02040 wrow++;
02041 }
02042 }
02043
02044 // Now the variance for mixture is calculated
02045 wrow=0;
02046 for (m_q=0; m_q<tn_qtl; m_q++){
02047 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of mother
02048 wcol=0;
02049 for (f_q=0; f_q<tn_qtl; f_q++){
02050 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of father
02051 v += pm[wrow][wcol].tsq*m_weight[wrow][wcol]+m_weight[wrow][wcol]*
02052 (pm[wrow][wcol].nu-nu)*(pm[wrow][wcol].nu-nu);
02053 wcol++;
02054 }
02055 }
02056 wrow++;
02057 }
02058 }
02059
02060 // k is the constant for a univariate normal (includes 1/sqrt(2*pi*sigma^2))
02061
02062 k = std::log(sum2) + sum1 - std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
02063
02064 temp_I->mix_anterior[i_q][i_sw].nu = nu;
02065 temp_I->mix_anterior[i_q][i_sw].tsq = 1.0/v;
02066 temp_I->mix_anterior[i_q][i_sw].k = k;
02067 // std::cout << i_q << " " << i_sw << " sum2 " << log(sum2) << " sum1 " << sum1 << " k " << k << std::endl;
02068 }
02069 }
02070 }
|
|
|
Definition at line 2289 of file nufamily.cpp. References matvec::Individual::anterior_iw, father_indx, matvec::Individual::initial_multi_m_anterior(), matvec::Individual::initial_multi_m_posterior(), mother_indx, myfather, mymother, myoffspring, numoffs, and matvec::Individual::posterior_iw. Referenced by matvec::Population::multi_m_log_likelihood_peeling().
02289 {
02290 // Initialize ALL members of the Nuclear family.
02291 // Else values are wrong for profile and maximization steps!
02292 // Would like to do just founders but cut individuals do not have this flag set.
02293
02294 int nchild;
02295 Individual *I;
02296
02297 // if (myfather->base()) {
02298 myfather->initial_multi_m_posterior(t_qtl,-1);
02299 myfather->initial_multi_m_anterior(t_qtl);
02300 myfather->posterior_iw[mother_indx] = 0;
02301 // }
02302
02303 // if (mymother->base()) {
02304 mymother->initial_multi_m_posterior(t_qtl,-1);
02305 mymother->initial_multi_m_anterior(t_qtl);
02306 mymother->posterior_iw[father_indx] = 0;
02307 // }
02308
02309 for (nchild=0; nchild < numoffs; nchild++) {
02310 I = myoffspring[nchild];
02311 I->initial_multi_m_posterior(t_qtl,-1);
02312 I->initial_multi_m_anterior(t_qtl);
02313 I->anterior_iw = 0;
02314 }
02315 }
|
|
|
Definition at line 2094 of file nufamily.cpp. References matvec::DataNode::double_val(), kernal, m_weight, matvec::Population::mean_for_genotype, matvec::DataNode::missing, matvec::Individual::mix_anterior, matvec::Individual::mix_posterior, multi_m_posterior(), myfather, mymother, matvec::Individual::n_switches(), matvec::Individual::nspouse(), pop, matvec::Individual::record(), matvec::Population::residual_var, and tn_qtl. Referenced by matvec::Population::multi_m_log_likelihood_peeling().
02094 {
02095
02096 unsigned j,gi;
02097 int f_q, f_sw, f_switches;
02098 double v, nu, y, v_y, nu_y, k_y;
02099 double a11,a12,a22,k,sum1=0.0,sum2=0.0;
02100 double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0)); // pi = 2.0*std::asin(1.0)
02101 f_switches=myfather->n_switches();
02102 if (!kernal) throw exception(" NuFamily::log_likelihood(): must be kernal nufamily");
02103 multi_m_posterior(myfather,mymother);
02104
02105 // Collect results to integrate uf for each f_q
02106
02107 for (f_q=0; f_q < tn_qtl; f_q++){
02108 // contributions from penetrance function of father
02109 // skip if phenotype of father is missing
02110
02111 if (!(myfather->record()[0].missing)) {
02112 y = myfather->record()[0].double_val();
02113 v_y = 1.0/(pop->residual_var[0][0]);
02114 k_y = std::log(inverse_sqrt2pi*std::sqrt(v_y));
02115 nu_y = pop->mean_for_genotype[f_q] - y;
02116 }
02117 else{
02118 nu_y = 0.0;
02119 v_y = 0.0;
02120 k_y = 0.0;
02121 }
02122
02123 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of J
02124
02125 a11 = nu_y*nu_y*v_y;
02126 a12 = nu_y*v_y;
02127 a22 = v_y;
02128 k = k_y;
02129 // contributions from Father's anterior
02130
02131 v = myfather->mix_anterior[f_q][f_sw].tsq;
02132 nu = myfather->mix_anterior[f_q][f_sw].nu;
02133 k += myfather->mix_anterior[f_q][f_sw].k;
02134 a11 += nu*nu*v;
02135 a12 += -nu*v;
02136 a22 += v;
02137
02138 // Contributions from posteriors of father except mother.
02139 // Skip over any posteriors not yet available
02140 unsigned numspouse = myfather->nspouse();
02141 for (j=0; j<numspouse; j++){
02142 if ((myfather->mix_posterior[j].done)) {
02143 v = myfather->mix_posterior[j].postvec[f_q][f_sw].tsq;
02144 nu = myfather->mix_posterior[j].postvec[f_q][f_sw].nu;
02145 k += myfather->mix_posterior[j].postvec[f_q][f_sw].k;
02146 a11 += nu*nu*v;
02147 a12 += -nu*v;
02148 a22 += v;
02149 }
02150 }
02151
02152 // compute "log likelihood for gi"
02153
02154 v = 1.0/a22;
02155 m_weight[f_q][f_sw] = -0.5*(a11 - a12*a12*v)
02156 + k + std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
02157 sum1 += m_weight[f_q][f_sw];
02158 }
02159 }
02160
02161 // compute sum of the "likelihoods"
02162
02163 sum1 /= (tn_qtl*f_switches);
02164 for (f_q=0;f_q<tn_qtl;f_q++){
02165 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of J
02166 sum2 += (std::exp(m_weight[f_q][f_sw] - sum1));
02167 }
02168 }
02169
02170 return (std::log(sum2) + sum1);
02171 }
|
|
||||||||||||||||
|
Definition at line 1344 of file nufamily.cpp. References child_matrix, matvec::DataNode::double_val(), father_indx, matvec::Individual::id(), matvec::Individual::index_sw, matvec::Vector< Vector< Sym2x2 > >::initialize(), m_weight, matvec::Population::mean_for_genotype, matvec::DataNode::missing, matvec::Individual::mix_anterior, matvec::Individual::mix_posterior, matvec::MIXED_TOL, mother_indx, multi_sumint_offspring(), myfather, matvec::Individual::n_switches(), matvec::Individual::nspouse(), numoffs, offspring(), pm, pop, matvec::Population::popmember, matvec::Individual::record(), matvec::Population::residual_var, matvec::Population::size(), spouse_matrix, and tn_qtl. Referenced by multi_m_ant_post(), multi_m_log_likelihood(), and multi_m_terminal_peeling().
01344 {
01345 ///////////////////////////////////////////////////
01346 // I and J must be parents of this nuclei family //
01347 ///////////////////////////////////////////////////
01348 // adding the switch stuff
01349
01350
01351 // A separate posterior has to be calculated for each discrete
01352 // genotype of I.
01353 unsigned j, excI, jj, independent=0;
01354 double a11, a12, a13, a22, a23, a33,k,ve;
01355 double u11, u12, u22, sum1, sum2, v, nu,y, v_y, k_y, nu_y, temp;
01356 double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0)); // pi = 2.0*std::asin(1.0)
01357 int m_q, f_q, j_q, i_q, offspring;
01358 int popsize=pop->size();
01359 // indx[0] is the index of the mother in the father's spouse list
01360 // indx[1] is the index of the father in the mother's spouse list
01361
01362 if (I==myfather) {jj = mother_indx; excI = father_indx; }
01363 else {jj = father_indx; excI = mother_indx; }
01364
01365 // switch information for parents
01366
01367 int i_switches = I->n_switches();
01368 int j_switches = J->n_switches();
01369 int Iindex_sw=I->index_sw;
01370 int Jindex_sw=J->index_sw;
01371 int wrow, wcol, drow, dcol, i_sw, j_sw, tag_id;
01372 ve=pop->residual_var[0][0];
01373
01374 if (i_peel) {
01375 //Iterative peeling so have to ensure that we are dealing with the
01376 //original individual not the cut one.
01377 // check I
01378 tag_id=I->id();
01379 if (tag_id > popsize) {
01380 // this is a cut individual so need to get the original one
01381 tag_id -= popsize; // correct id;
01382 I=pop->popmember[tag_id];
01383 }
01384 // check J
01385 tag_id=J->id();
01386 if (tag_id > popsize) {
01387 // this is a cut individual so need to get the original one
01388 tag_id -= popsize; // correct id;
01389 J=pop->popmember[tag_id];
01390 }
01391 }
01392 // std::cout << "Post for " << pop->ind_name(I->id()) << " thru " << pop->ind_name(J->id()) << std::endl;
01393
01394 // First, integrate and sum over each offspring
01395 // Approximate results are stored for each i_q and j_q
01396 // in child_matrix
01397
01398 // Initialize child_matrix for this posterior
01399 wrow = 0;
01400 for (i_q=0;i_q<tn_qtl;i_q++){
01401 for (i_sw=0; i_sw < i_switches; i_sw++){ // loop for genotypes of I
01402 wcol = 0;
01403 for (j_q=0;j_q<tn_qtl; j_q++){
01404 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01405 // First make sure get right access to child_matrix
01406 // MOTHER is the ROW
01407 // FATHER is the COLUMN
01408 if (I == myfather) {
01409 drow=wcol;
01410 dcol=wrow;
01411 }
01412 else {
01413 drow=wrow;
01414 dcol=wcol;
01415 }
01416 child_matrix[drow][dcol].initialize();
01417 wcol++;
01418 }
01419 }
01420 wrow++;
01421 }
01422 }
01423
01424 for (offspring=0; offspring<numoffs; offspring++) { // sum and int. each offspring
01425 multi_sumint_offspring(offspring);
01426 }
01427
01428 // Collect results to integrate uj for each (j_q and sj) in spouse
01429
01430 for (j_q=0;j_q<tn_qtl;j_q++){
01431
01432 a11=a12=a13=a22=a23=a33=k=0.0;
01433
01434 // contributions from penetrance funtion of J
01435 // skip if phenotype of J is missing
01436
01437 if (!(J->record()[0].missing)) {
01438 y = J->record()[0].double_val();
01439 v_y = 1.0/ve;
01440 k_y = std::log(inverse_sqrt2pi*std::sqrt(v_y));
01441 nu_y = pop->mean_for_genotype[j_q] - y;
01442
01443 }
01444 else{
01445 nu_y = 0.0;
01446 v_y = 0.0;
01447 k_y = 0.0;
01448 }
01449
01450 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01451
01452 a11 = nu_y*nu_y*v_y;
01453 a13 = nu_y*v_y;
01454 a33 = v_y;
01455 k = k_y;
01456
01457 // contributions from J's anterior
01458
01459 v = J->mix_anterior[j_q][j_sw].tsq;
01460 nu = J->mix_anterior[j_q][j_sw].nu;
01461 k += J->mix_anterior[j_q][j_sw].k;
01462 a11 += nu*nu*v;
01463 a13 += -nu*v;
01464 a33 += v;
01465
01466
01467 // Contributions from posteriors of J except I.
01468 // Skip over any posteriors not yet available
01469
01470 unsigned numspouse = J->nspouse();
01471 for (j=0; j<numspouse; j++){
01472 if ((J->mix_posterior[j].done) && (j!=excI) ) {
01473 v = J->mix_posterior[j].postvec[j_q][j_sw].tsq;
01474 nu = J->mix_posterior[j].postvec[j_q][j_sw].nu;
01475 k += J->mix_posterior[j].postvec[j_q][j_sw].k;
01476 a11 += nu*nu*v;
01477 a13 += -nu*v;
01478 a33 += v;
01479 }
01480 }
01481
01482 // Put results into spouse_matrix
01483
01484 spouse_matrix[j_q][j_sw].a11 = a11;
01485 spouse_matrix[j_q][j_sw].a13 = a13;
01486 spouse_matrix[j_q][j_sw].a33 = a33;
01487 spouse_matrix[j_q][j_sw].k = k;
01488 // std::cout << "Spouse matrix j_q " << j_q << " j_sw " << j_sw << " a11 " << a11 << " a13 " << a13 << " a33 " << a33 << " K " << k << std::endl;
01489 }
01490 }
01491
01492 // For each (i_q,i_sw) and (j_q,j_sw) complete integration of uj and store results
01493 // in pm
01494 wrow = 0;
01495 for (i_q=0;i_q<tn_qtl;i_q++){
01496 for (i_sw=0; i_sw < i_switches; i_sw++){ // loop for genotypes of I
01497 wcol = 0;
01498 for (j_q=0;j_q<tn_qtl; j_q++){
01499 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01500
01501 a11=a12=a13=a22=a23=a33=k=0.0;
01502
01503 // Get contributions from spouse_matrix for genotype j_q
01504
01505 a11 = spouse_matrix[j_q][j_sw].a11;
01506 a13 = spouse_matrix[j_q][j_sw].a13;
01507 a33 = spouse_matrix[j_q][j_sw].a33;
01508 k = spouse_matrix[j_q][j_sw].k;
01509
01510 // Get contributions from child_matrix for genotypes i_q and j_q
01511 // First make sure get right access to child_matrix
01512 // MOTHER is the ROW
01513 // FATHER is the COLUMN
01514 if (I == myfather) {
01515 drow=wcol;
01516 dcol=wrow;
01517 }
01518 else {
01519 drow=wrow;
01520 dcol=wcol;
01521 }
01522 a11 += child_matrix[drow][dcol].a11;
01523 a12 += child_matrix[drow][dcol].a12;
01524 a13 += child_matrix[drow][dcol].a12;
01525 a22 += child_matrix[drow][dcol].a22;
01526 a23 += child_matrix[drow][dcol].a22;
01527 a33 += child_matrix[drow][dcol].a22;
01528 k += child_matrix[drow][dcol].k;
01529
01530 // complete integration of uj
01531
01532 v = 1.0/a33;
01533 u11 = a11 - a13*a13*v;
01534 u12 = a12 - a13*a23*v;
01535 u22 = a22 - a23*a23*v;
01536 k += std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01537
01538 if (u22 < MIXED_TOL) {
01539 pm[wrow][wcol].k = k - 0.5*u11;
01540 independent = 1;
01541 }
01542 else {
01543
01544 // Calculate Normal mean, variance, and k for i_q and j_q
01545
01546 v = 1.0/u22;
01547 nu = -u12*v;
01548 pm[wrow][wcol].tsq = v;
01549 pm[wrow][wcol].nu = nu;
01550 pm[wrow][wcol].k = k - 0.5*(u11 - nu*nu*u22) + std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01551 }
01552 wcol++;
01553 }
01554 }
01555 wrow++;
01556 }
01557 }
01558 if (independent) {
01559
01560 // The posterior and anterior likelihoods are independent
01561 // This happens due to missing data
01562 wrow=0;
01563 for (i_q=0; i_q<tn_qtl; i_q++){
01564 for (i_sw=0; i_sw < i_switches; i_sw++){ // loop for genotypes of I
01565 wcol=0;
01566 sum1=0.0;
01567 for(j_q=0;j_q<tn_qtl; j_q++){
01568 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01569 sum1 += pm[wrow][wcol].k;
01570 wcol++;
01571 }
01572 }
01573 sum1 /= (tn_qtl*j_switches);
01574 wcol=0;
01575 sum2=0.0;
01576 for(j_q=0;j_q<tn_qtl; j_q++){
01577 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01578 sum2 += std::exp(pm[wrow][wcol].k - sum1);
01579 wcol++;
01580 }
01581 }
01582
01583 I->mix_posterior[jj].postvec[i_q][i_sw].nu = 0.0;
01584 I->mix_posterior[jj].postvec[i_q][i_sw].tsq = 0.0;
01585 I->mix_posterior[jj].postvec[i_q][i_sw].k = std::log(sum2) + sum1;
01586 wrow++;
01587 }
01588 }
01589 I->mix_posterior[jj].done=1;
01590 return;
01591 }
01592
01593
01594 // Now for each i_q sum over j_q. The mixture normal that results from summing over
01595 // j_q is approximated by an univariate normal with the same mean and variance as
01596 // the mixture
01597
01598
01599 wrow=0;
01600 for (i_q=0; i_q<tn_qtl; i_q++){
01601 for (i_sw=0; i_sw < i_switches; i_sw++){ // loop for genotypes of I
01602 // Calculate weights of mixture after scaling by sum1
01603 wcol=0;
01604 sum1=0.0;
01605 for(j_q=0;j_q<tn_qtl; j_q++){
01606 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01607 sum1 += pm[wrow][wcol].k;
01608 wcol++;
01609 }
01610 }
01611 sum1 /= (tn_qtl*j_switches);
01612 sum2=0.0;
01613 wcol=0;
01614 for(j_q=0;j_q<tn_qtl; j_q++){
01615 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01616 temp= std::exp(pm[wrow][wcol].k - sum1);
01617 sum2 +=temp;
01618 m_weight[j_q][j_sw]=temp;
01619 wcol++;
01620 }
01621 }
01622 for (j_q=0; j_q<tn_qtl;j_q++){
01623 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01624 m_weight[j_q][j_sw] /= sum2;
01625 }
01626 }
01627
01628 // Now the mean and variance for mixture are calculated
01629
01630 nu = v = k = 0;
01631 wcol=0;
01632 for (j_q=0; j_q<tn_qtl;j_q++) {
01633 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01634 nu += pm[wrow][wcol].nu*m_weight[j_q][j_sw];
01635 wcol++;
01636 }
01637 }
01638 wcol=0;
01639 for (j_q=0; j_q<tn_qtl;j_q++){
01640 for (j_sw=0; j_sw < j_switches; j_sw++){ // loop for genotypes of J
01641 v += pm[wrow][wcol].tsq*m_weight[j_q][j_sw]
01642 + m_weight[j_q][j_sw]*(pm[wrow][wcol].nu - nu)*(pm[wrow][wcol].nu - nu);
01643 wcol++;
01644 }
01645 }
01646 // k is the constant for a univariate normal (includes 1/sqrt(2*pi*sigma^2))
01647
01648 k = std::log(sum2) + sum1 - std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01649 I->mix_posterior[jj].postvec[i_q][i_sw].nu = nu;
01650 I->mix_posterior[jj].postvec[i_q][i_sw].tsq = 1.0/v;
01651 I->mix_posterior[jj].postvec[i_q][i_sw].k = k;
01652 wrow++;
01653 // std::cout << i_q << " " << i_sw << " sum2 " << log(sum2) << " sum1 " << sum1 << " k " << k << std::endl;
01654 }
01655 }
01656 I->mix_posterior[jj].done=1;
01657 }
|
|
|
Definition at line 2073 of file nufamily.cpp. References connectors, father_indx, mother_indx, multi_m_anterior(), multi_m_posterior(), myfather, mymother, and matvec::Individual::posterior_iw. Referenced by matvec::Population::multi_m_log_likelihood_peeling().
02073 {
02074 ////////////////////////////////////////////////////////////
02075 // REQUIREMENTS: initilization for anterior and posterior
02076 ////////////////////////////////////////////////////////////
02077 std::list<Individual *>::iterator pos;
02078 pos = connectors.begin();
02079 if (pos == connectors.end()) throw exception(" NuFamily::mutlti_m_terminal_peeling(): no connector");
02080 if (*pos == myfather) {
02081 multi_m_posterior(*pos, mymother);
02082 myfather->posterior_iw[mother_indx] = iw;
02083 }
02084 else if (*pos == mymother) {
02085 multi_m_posterior(*pos,myfather);
02086 mymother->posterior_iw[father_indx] = iw;
02087 }
02088 else {
02089 multi_m_anterior(*pos);
02090 (*pos)->anterior_iw = iw;
02091 }
02092 }
|
|
||||||||||||||||||||
|
Definition at line 633 of file nufamily.cpp. References father_indx, get_m_posterior(), matvec::Individual::id(), matvec::Individual::m_anterior, matvec::Individual::m_anterior_scale, matvec::Individual::m_posterior, matvec::Individual::m_posterior_scale, mother_indx, multi_fullsibs_prob(), myfather, matvec::Individual::n_switches(), matvec::Individual::nspouse(), matvec::Population::P_ndim, pop, matvec::Population::popmember, matvec::Population::size(), and wspace. Referenced by multi_ant_post(), multi_llh(), and multi_terminal_peeling().
00635 {
00636 /////////////////////////////////////////////////////
00637 // I and J must be parents of this nuclei family
00638 // posterior should have enouph space
00639 /////////////////////////////////////////////////////
00640
00641 unsigned jj,ui,uj,excI,i,j,k,i_sw,i_p,i_q,j_sw,j_p,j_q,i_w,j_w,tag_id,popsize,ndim, drow,dcol;
00642 unsigned n_switches = I->n_switches();
00643 unsigned j_switches = J->n_switches();
00644 double val, scale=0, *ws, *post_j, i_sum=0, temp;
00645
00646 popsize=pop->size();
00647 ndim = pop->P_ndim;
00648
00649 if (I==myfather) {
00650 jj = mother_indx;
00651 excI = father_indx;
00652 scale = (multi_fullsibs_prob(0,post_mat,i_peel));
00653 // std::cout << "scale in mp " << setprecision(14) << scale << std::endl;
00654 }
00655 else {
00656 jj = father_indx;
00657 excI = mother_indx;
00658 scale = multi_fullsibs_prob(0,post_mat,i_peel);
00659 }
00660 // std::cout << wspace << std::endl;
00661 if (i_peel) {
00662 //Iterative peeling so have to ensure that we are dealing with the
00663 //original offspring not the cut one.
00664 // check I
00665 tag_id=I->id();
00666 if (tag_id > popsize) {
00667 // this is a cut individual so need to get the original one
00668 tag_id -= popsize; // correct id;
00669 I=pop->popmember[tag_id-1];
00670 }
00671 // check J
00672 tag_id=J->id();
00673 if (tag_id > popsize) {
00674 // this is a cut individual so need to get the original one
00675 tag_id -= popsize; // correct id;
00676 J=pop->popmember[tag_id-1];
00677 }
00678 }
00679 for (k=0;k<ndim;k++){
00680 for (i=0;i<j_switches; i++){
00681 for (j=0;j<4;j++){
00682 post_mat[k][i][j] = J->m_anterior[k][i][j];
00683 }
00684 }
00685 }
00686 scale += J->m_anterior_scale;
00687 // std::cout << "ant scale " << J->m_anterior_scale << " post mat " << post_mat << " Workspace " << wspace << std::endl;
00688 // std::cout << "scale after anterior added " << scale << std::endl;
00689 if (J->nspouse() > 1) {
00690 temp= get_m_posterior(J,excI,post_mat);
00691 scale += temp;
00692 }
00693
00694 // std::cout << "scale after get m posterior " << scale << std::endl;
00695 i_w = 0;
00696 for (i_p=0; i_p<ndim; i_p++){
00697 for (i_q=0; i_q<4; i_q++){
00698 for (i_sw=0;i_sw<n_switches; i_sw++) {
00699 j_w = 0;
00700 val = 0.0;
00701 for (j_p=0; j_p<ndim; j_p++){
00702 for (j_q=0; j_q<4; j_q++){
00703 for (j_sw=0;j_sw<j_switches; j_sw++) {
00704 if (I == myfather) { // mother is the row, father is the column
00705 drow=j_w;
00706 dcol=i_w;
00707 }
00708 else {
00709 drow=i_w;
00710 dcol=j_w;
00711 }
00712
00713 val += (post_mat[j_p][j_sw][j_q]*wspace[drow][dcol]);
00714 // std::cout << j_p << " " << j_q << " " << j_sw << " " << post_mat[j_p][j_sw][j_q] << " " << wspace[i_w][j_w] << std::endl;
00715 j_w++;
00716 }
00717 }
00718 }
00719 i_w++;
00720 I->m_posterior[jj][i_p][i_sw][i_q] = val;
00721 i_sum += val;
00722 // std::cout << " jj " << jj << " i_sw " << i_sw << " i_q " << i_q << " m_post " << I->m_posterior[jj][i_sw][i_q] << std::endl;
00723 }
00724 }
00725 }
00726
00727 for (i_p=0; i_p<ndim; i_p++){
00728 for (i_q=0; i_q<4; i_q++){
00729 for (i_sw=0;i_sw<n_switches; i_sw++) {
00730 I->m_posterior[jj][i_p][i_sw][i_q] /= i_sum;
00731 }
00732 }
00733 }
00734 scale += std::log(i_sum);
00735 I->m_posterior_scale[jj] = scale;
00736 // std::cout << " Posterior scale " << scale << std::endl;
00737 }
|
|
||||||||||||
|
Definition at line 1140 of file nufamily.cpp. References child_matrix, matvec::DataNode::double_val(), matvec::Population::F, matvec::Individual::id(), matvec::Individual::index_sw, m_weight, matvec::Population::mean_for_genotype, matvec::DataNode::missing, matvec::Individual::mix_posterior, matvec::MIXED_TOL, multi_get_tr(), myfather, mymother, myoffspring, matvec::Individual::n_switches(), matvec::Individual::nspouse(), pop, matvec::Population::popmember, matvec::Individual::record(), matvec::Population::residual_var, matvec::Population::size(), matvec::Population::switch_table, tn_qtl, tr, matvec::Fpmm::var, and wksp_for_gen. Referenced by multi_m_anterior(), and multi_m_posterior().
01140 {
01141 //////////////////////////////////////////////////////////////////////////
01142 // Summation over discrete genotypes and integration over continuous //
01143 // genotypes for offspring i. Integration is done first. Then, for each //
01144 // discrete genotype of mother and father, the mixture that results //
01145 // from summing over the discrete genotype of i is approximated by a //
01146 // univariate normal with the same mean and variance as the mixture. //
01147 // The contributions from this approximation to the posterior or //
01148 // anterior are accumulated in "unormal_aprox_for_gen" //
01149 //////////////////////////////////////////////////////////////////////////
01150
01151 // We are adding the summataion over switches to the previous code
01152
01153 unsigned j,independent;
01154 Individual *child = myoffspring[i];
01155 double inverse_sqrt2pi = 1.0/std::sqrt(4.0*std::asin(1.0)); // pi = 2.0*std::asin(1.0)
01156 double u11,u12,u22,v,nu, ve;
01157 double a11,a12,a13,a22,a23,a33,k;
01158 double P_var = pop->F->var; // polygenic varaince
01159 int i_q=0,m_q=0,f_q=0;
01160 int Iindex_sw=child->index_sw;
01161 int i_switches = child->n_switches();
01162 double sum1=0.0, sum2=0.0;
01163 double a33_t, k_t, k_q, a11_q, a13_q, a33_q;
01164 int i_sw, tag_id;
01165 int popsize=pop->size();
01166 ve=pop->residual_var[0][0];
01167 if (i_peel) {
01168 //Iterative peeling so have to ensure that we are dealing with the
01169 //original offspring not the cut one.
01170 tag_id=child->id();
01171 if (tag_id > popsize) {
01172 // this is a cut individual so need to get the original one
01173 tag_id -= popsize; // correct id;
01174 child=pop->popmember[tag_id];
01175 }
01176 }
01177
01178 ////// contributions from transition function of offspring
01179
01180 v = 1.0/(0.5*P_var);
01181 a22 = 0.25*v;
01182 a23 = -0.5 *v;
01183 a33_t = v;
01184 k_t = std::log(inverse_sqrt2pi*std::sqrt(v));
01185
01186 for (i_q=0;i_q<tn_qtl;i_q++) {
01187 ///// contributions from penetrance function of offspring
01188 ///// skip if phenotype of offspring is missing
01189
01190 if (!(child->record()[0].missing)) {
01191 v = 1.0/ve;
01192 nu = pop->mean_for_genotype[i_q] - child->record()[0].double_val();
01193 k_q = std::log(inverse_sqrt2pi*std::sqrt(v));
01194 a11_q = nu*nu*v;
01195 a13_q = nu*v;
01196 a33_q = v;
01197 }
01198 else {
01199 k_q = 0.0;
01200 a11_q = 0.0;
01201 a13_q = 0.0;
01202 a33_q = 0.0;
01203 }
01204
01205 for (i_sw=0; i_sw < i_switches; i_sw++) {
01206 a11=a12=a13=a33=k=0.0;
01207
01208 ///// contributions from posteriors of offspring.
01209 ///// Skip any posteriors that are not yet calculated.
01210 ///// Skipping will happen only in iterative peeling.
01211 ///// In termainal and recursive peeling, posteriors
01212 ///// are always available when required.
01213
01214 unsigned numspouse = child->nspouse();
01215
01216 for (j=0; j<numspouse; j++){
01217 if (child->mix_posterior[j].done) {
01218 v = child->mix_posterior[j].postvec[i_q][i_sw].tsq;
01219 nu = child->mix_posterior[j].postvec[i_q][i_sw].nu;
01220 k += child->mix_posterior[j].postvec[i_q][i_sw].k;
01221 a11 += nu*nu*v;
01222 a13 += -nu*v;
01223 a33 += v;
01224 }
01225 }
01226
01227 // Accumulate a's where needed
01228 a11 += (a11_q);
01229 a13 += (a13_q);
01230 a33 += (a33_q+a33_t);
01231 k += (k_q+k_t);
01232 ///// Now calculate U matrix and K
01233
01234 v = 1.0/(a33);
01235 u11 = a11 - a13*a13*v;
01236 u12 = a12 - a13*a23*v;
01237 u22 = a22 - a23*a23*v;
01238 k += std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01239 independent = 0;
01240 if (u22*u22 < MIXED_TOL) {
01241 wksp_for_gen[i_q][i_sw].k = k - 0.5*u11;
01242 sum1 += wksp_for_gen[i_q][i_sw].k;
01243 independent =1;
01244 }
01245 else {
01246 //// Calculate Normal mean, variance, and k for genotype g
01247 v = 1.0/u22;
01248 nu = -u12*v;
01249 wksp_for_gen[i_q][i_sw].tsq = v;
01250 wksp_for_gen[i_q][i_sw].nu = nu;
01251 wksp_for_gen[i_q][i_sw].k = k - 0.5*(u11 - nu*nu*u22) + std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01252 sum1 += wksp_for_gen[i_q][i_sw].k;
01253 }
01254 }
01255 }
01256
01257 //// For each genotype of mother (gm) and of father (f_q):
01258 //// Calculate mean and variance of mixture. Also calculate k for mixture.
01259 //// First calculate weights. sum1 is for scaling; scaling does not
01260 //// affect weights
01261
01262
01263
01264 sum1 /= (tn_qtl*i_switches);
01265
01266 // now wksp_for_gen has exp(wksp_for_gen.k - sum1) !!!
01267
01268 for (i_q=0; i_q<tn_qtl;i_q++) {
01269 for (i_sw=0; i_sw < i_switches; i_sw++) {
01270 wksp_for_gen[i_q][i_sw].k = std::exp(wksp_for_gen[i_q][i_sw].k - sum1);
01271 }
01272 }
01273
01274 int m_switches = mymother->n_switches();
01275 int f_switches = myfather->n_switches();
01276 int Findex_sw=myfather->index_sw;
01277 int Mindex_sw=mymother->index_sw;
01278 int m_sw, Mswitch, f_sw, Fswitch;
01279
01280 // Loops over the switches for parents go here
01281
01282 int w_row=0, w_col;
01283 for (m_q=0; m_q<tn_qtl;m_q++) {
01284 for (m_sw=0; m_sw < m_switches; m_sw++){ // loop for genotypes of mother
01285 w_col = 0;
01286 Mswitch=pop->switch_table[Mindex_sw][m_sw+1]; // get current dam switch
01287 for (f_q=0; f_q<tn_qtl;f_q++) {
01288 for (f_sw=0; f_sw < f_switches; f_sw++){ // loop for genotypes of father
01289 Fswitch=pop->switch_table[Findex_sw][f_sw+1]; // switch for father (f_sw)
01290 multi_get_tr(child,m_q,f_q,Mswitch,Fswitch); // results are stored in tr
01291 // std::cout << tr;
01292 sum2=0.0;
01293 for (i_q=0; i_q<tn_qtl;i_q++) {
01294 for (i_sw=0; i_sw < i_switches; i_sw++) {
01295 m_weight[i_q][i_sw]= tr[i_q][i_sw] * wksp_for_gen[i_q][i_sw].k;
01296 // std::cout << i_q << " i_sw " << i_sw << " m_ weight = " << m_weight[i_q][i_sw] << " tr " << tr[i_q][i_sw] << " wksp_for_gen " << wksp_for_gen[i_q][i_sw].k << std::endl;
01297 sum2 += m_weight[i_q][i_sw];
01298 }
01299 }
01300 // std::cout << " sum1= " << sum1 << " sum2= " << sum2 << std::endl;
01301 //// Now the mean and variance for mixture are calculated
01302 if (!independent) {
01303 nu=0.0;
01304 v=0.0;
01305 for (i_q=0; i_q<tn_qtl;i_q++) {
01306 for (i_sw=0; i_sw < i_switches; i_sw++) {
01307 m_weight[i_q][i_sw] /= sum2;
01308 nu += wksp_for_gen[i_q][i_sw].nu *m_weight[i_q][i_sw];
01309 }
01310 }
01311 for (i_q=0; i_q<tn_qtl;i_q++) {
01312 for (i_sw=0; i_sw < i_switches; i_sw++) {
01313 v += wksp_for_gen[i_q][i_sw].tsq*m_weight[i_q][i_sw]
01314 + m_weight[i_q][i_sw]*(wksp_for_gen[i_q][i_sw].nu - nu)
01315 *(wksp_for_gen[i_q][i_sw].nu - nu);
01316 }
01317 }
01318 //// k is the constant for a univariate normal (includes 1/sqrt(2*pi*sigma^2))
01319
01320 k = std::log(sum2) + sum1 - std::log(std::sqrt(4.0*std::asin(1.0)*v)); // pi = 2.0*std::asin(1.0)
01321
01322 //// contributions to anterior or posterior are accumulated in child_matrix
01323
01324 child_matrix[w_row][w_col].a11 += nu*nu/v;
01325 child_matrix[w_row][w_col].a12 += -nu/v;
01326 child_matrix[w_row][w_col].a22 += 1.0/v;
01327 child_matrix[w_row][w_col].k += k;
01328 }
01329 else {
01330 child_matrix[w_row][w_col].k += std::log(sum2) + sum1;
01331 }
01332 w_col++;
01333 }
01334 }
01335 w_row++;
01336 }
01337 }
01338 }
|
|
||||||||||||||||
|
Definition at line 612 of file nufamily.cpp. References connectors, father_indx, mother_indx, multi_anterior(), multi_posterior(), myfather, mymother, and matvec::Individual::posterior_iw. Referenced by matvec::Population::multipoint_likelihood().
00612 {
00613
00614 // std::cout << "entering multi_terminal_peeling " << std::endl;
00615 std::list<Individual *>::iterator pos;
00616 pos = connectors.begin();
00617 if (pos == connectors.end()) throw exception(" NuFamily::multi_terminal_peeling(): no connector");
00618 if (*pos == myfather) {
00619 multi_posterior(myfather,mymother,post_mat_f);
00620 myfather->posterior_iw[mother_indx] = iw;
00621 }
00622 else if (*pos == mymother) {
00623 multi_posterior(mymother,myfather,post_mat_m);
00624 mymother->posterior_iw[father_indx] = iw;
00625 }
00626 else {
00627 multi_anterior(*pos, post_mat_m, post_mat_f);
00628 (*pos)->anterior_iw = iw;
00629 }
00630 // std::cout << "exiting multi_terminal_peeling " << std::endl;
00631 }
|
|
||||||||||||
|
Definition at line 2173 of file nufamily.cpp. References matvec::Matrix< double >::resize(). Referenced by matvec::Population::maxant_maxpost(), and matvec::Population::multi_m_log_likelihood_peeling().
02173 {
02174 int i,j,k;
02175
02176 NuFamily::child_matrix.resize(t_qtl*max_switches);
02177 NuFamily::pm.resize(t_qtl*max_switches);
02178
02179 for (i=0; i< (t_qtl*max_switches); i++){
02180 NuFamily::child_matrix[i].resize(t_qtl*max_switches);
02181 NuFamily::pm[i].resize(t_qtl*max_switches);
02182 }
02183
02184 NuFamily::wksp_for_gen.resize(t_qtl);
02185 NuFamily::spouse_matrix.resize(t_qtl);
02186 NuFamily::myfather_matrix.resize(t_qtl);
02187 NuFamily::mymother_matrix.resize(t_qtl);
02188
02189 for (i=0; i < t_qtl; i++ ) {
02190 NuFamily::wksp_for_gen[i].resize(max_switches);
02191 NuFamily::spouse_matrix[i].resize(max_switches);
02192 NuFamily::myfather_matrix[i].resize(max_switches);
02193 NuFamily::mymother_matrix[i].resize(max_switches);
02194 }
02195
02196 NuFamily::m_weight.resize(t_qtl*max_switches,t_qtl*max_switches);
02197 NuFamily::tr.resize(t_qtl,max_switches,0.0);
02198 }
|
|
|
Definition at line 252 of file nufamily.cpp. References connectors. Referenced by matvec::Population::detect_loop().
00253 {
00254 std::list<Individual *>::iterator pos;
00255 pos=connectors.begin(); //BRS flag
00256 while (pos != connectors.end()) {
00257 if ((*pos)->group_id <= 1) { pos = connectors.erase(pos);}
00258 else { pos++; } //BRS
00259 }
00260 return connectors.size();
00261 }
|
|
|
Definition at line 119 of file nufamily.h. References numoffs.
00119 {return numoffs;}
|
|
|
Definition at line 99 of file nufamily.h. References numoffs. Referenced by matvec::Population::build_nufamily(), and offspring().
00099 {numoffs = n;}
|
|
|
Definition at line 122 of file nufamily.h.
00122 {return myoffspring;}
|
|
||||||||||||
|
Definition at line 157 of file nufamily.cpp. References myoffspring, noffs(), and numoffs. Referenced by matvec::Population::build_nufamily(), multi_m_anterior(), and multi_m_posterior().
00158 {
00159 if (myoffspring) {
00160 delete [] myoffspring;
00161 myoffspring=0;
00162 }
00163 numoffs = noffs;
00164 if(numoffs>0){
00165 myoffspring = new Individual* [numoffs];
00166 }
00167 else {
00168 myoffspring = 0;
00169 }
00170 for (unsigned i=0; i<numoffs; i++) myoffspring[i] = offs[i];
00171 }
|
|
||||||||||||||||||||
|
Definition at line 449 of file nufamily.cpp. References matvec::Individual::anterior, matvec::Vector< double >::begin(), father_indx, fullsibs_prob(), get_posterior(), mother_indx, mymother, matvec::Individual::nspouse(), matvec::Individual::posterior, and tn_genotype. Referenced by log_likelihood(), and terminal_peeling().
00451 {
00452 /////////////////////////////////////////////////////
00453 // I and J must be parents of this nuclei family
00454 // posterior should have enouph space
00455 /////////////////////////////////////////////////////
00456 if (!trans) throw exception("NuFamily::posterior(arg1,arg2,arg3,arg4): null arg3");
00457 unsigned jj,ui,uj,excI;
00458 double val, scale, *ws, *post_j;
00459 Vector<double> post_vec(tn_genotype);
00460 Vector<double> apj(tn_genotype);
00461
00462 if (I==mymother) {jj = father_indx; excI = mother_indx; }
00463 // if (I==myfather) {jj = mates_indx[0]; excI = mates_indx[1]; } //old version
00464 else {jj = mother_indx; excI = father_indx; }
00465
00466 scale = fullsibs_prob(0,trans,WSpace);
00467 for (uj=0; uj<tn_genotype; uj++) apj[uj] = J->anterior[uj];
00468 scale += J->anterior[tn_genotype];
00469 if (J->nspouse() > 1) {
00470 scale += get_posterior(J,excI,post_vec);
00471 for (uj=0; uj<tn_genotype; uj++) apj[uj] *= post_vec[uj];
00472 }
00473 post_j = I->posterior[jj].begin();
00474 for (ui=0; ui<tn_genotype; ui++) {
00475 ws = WSpace[ui];
00476 for (val=0.0, uj=0; uj<tn_genotype; uj++) val += *ws++ * apj[uj];
00477 post_j[ui] = val;
00478 }
00479 for (val=0.0, uj=0; uj<tn_genotype; uj++) val += post_j[uj];
00480 for (uj=0; uj<tn_genotype; uj++) post_j[uj] /= val;
00481 scale += std::log(val);
00482 post_j[tn_genotype] = scale;
00483 }
|
|
||||||||||||
|
Definition at line 230 of file nufamily.cpp. References matvec::Vector< T >::begin(), myfather, mymother, myoffspring, numoffs, matvec::Individual::pretend_missing(), and matvec::Vector< T >::size().
00231 {
00232 unsigned tng = (unsigned) genotype_freq.size();
00233 const double* gfreq = genotype_freq.begin();
00234 myfather->pretend_missing(on,gfreq,tng );
00235 mymother->pretend_missing(on,gfreq,tng);
00236 for (unsigned i=0; i<numoffs; i++) {
00237 myoffspring[i]->pretend_missing(on,gfreq,tng);
00238 }
00239 }
|
|
|
Definition at line 2200 of file nufamily.cpp. References myfather, mymother, myoffspring, numoffs, matvec::Individual::pretend_multi_m_missing(), and tn_qtl.
02200 {
02201 // unsigned tng = (unsigned) genotype_freq.size();
02202 // const double* gfreq = genotype_freq.ve;
02203 myfather->pretend_multi_m_missing(on,tn_qtl);
02204 mymother->pretend_multi_m_missing(on,tn_qtl);
02205 for (unsigned i=0; i<numoffs; i++) {
02206 myoffspring[i]->pretend_multi_m_missing(on,tn_qtl);
02207 }
02208 }
|
|
|
Definition at line 1094 of file nufamily.cpp. References myfather, mymother, myoffspring, numoffs, penetrance, and matvec::Individual::pretend_multi_missing().
01094 {
01095 // unsigned tng = (unsigned) genotype_freq.size();
01096 // const double* gfreq = genotype_freq.ve;
01097 myfather->pretend_multi_missing(on,penetrance);
01098 mymother->pretend_multi_missing(on,penetrance);
01099 for (unsigned i=0; i<numoffs; i++) {
01100 myoffspring[i]->pretend_multi_missing(on,penetrance);
01101 }
01102 }
|
|
|
Definition at line 119 of file nufamily.cpp. References matvec::Individual::id(), myfather, mymother, myoffspring, numoffs, pop, and matvec::Population::size(). Referenced by ~NuFamily().
00120 {
00121 if (pop) {
00122 unsigned popsize = pop->size();
00123 if (mymother->id() > popsize) {
00124 if(mymother){
00125 delete mymother;
00126 mymother=0;
00127 }
00128 }
00129 if (myfather->id() > popsize) {
00130 if(myfather){
00131 delete myfather;
00132 myfather=0;
00133 }
00134 }
00135 if (myoffspring) {
00136 for (unsigned i=0; i<numoffs; i++) {
00137 if (myoffspring[i]->id() > popsize) {
00138 if(myoffspring[i]){
00139 delete myoffspring[i];
00140 myoffspring[i]=0;
00141 }
00142 }
00143 }
00144 if(myoffspring){
00145 delete [] myoffspring;
00146 myoffspring=0;
00147 }
00148 }
00149 }
00150 else {
00151 if (myoffspring) {
00152 delete [] myoffspring; myoffspring = 0;
00153 }
00154 }
00155 }
|
|
||||||||||||||||
|
Definition at line 518 of file nufamily.cpp. References anterior(), connectors, father_indx, mother_indx, myfather, mymother, posterior(), and matvec::Individual::posterior_iw.
00519 {
00520 ///////////////////////////////////////////////////////////
00521 // REQUIREMENTS: initilization for anterior and posterior
00522 ///////////////////////////////////////////////////////////
00523 if (!trans) throw exception("NuFamily::terminal_peeling(trans,...): null trans");
00524 std::list<Individual *>::iterator pos;
00525 pos = connectors.begin();
00526 if (pos == connectors.end()) throw exception(" NuFamily::terminal_peeling(): no connector");
00527 if (*pos == mymother) {
00528 posterior(*pos,myfather,trans,WSpace);
00529 mymother->posterior_iw[father_indx] = iw;
00530 }
00531 else if (*pos == myfather) {
00532 posterior(*pos,mymother,trans,WSpace);
00533 myfather->posterior_iw[mother_indx] = iw;
00534 }
00535 else {
00536 anterior(*pos,trans,WSpace);
00537 (*pos)->anterior_iw = iw;
00538 }
00539 }
|
|
||||||||||||
|
Definition at line 498 of file nufamily.cpp. References anterior(), connectors, myfather, mymother, and posterior().
00499 {
00500 ////////////////////////////////////////////////////////////
00501 // REQUIREMENTS: initilization for anterior and posterior
00502 ////////////////////////////////////////////////////////////
00503 if (!trans) throw exception("NuFamily::terminal_peeling(trans,...): null trans");
00504 std::list<Individual *>::iterator pos;
00505 pos = connectors.begin();
00506 if (pos == connectors.end()) throw exception(" NuFamily::terminal_peeling(): no connector");
00507 if (*pos == mymother) {
00508 posterior(*pos,myfather,trans,WSpace);
00509 }
00510 else if (*pos == myfather) {
00511 posterior(*pos,mymother,trans,WSpace);
00512 }
00513 else {
00514 anterior(*pos,trans,WSpace);
00515 }
00516 }
|
|
|
Definition at line 77 of file nufamily.h. |
|
|
Definition at line 39 of file nufamily.cpp. Referenced by multi_m_anterior(), multi_m_posterior(), and multi_sumint_offspring(). |
|
|
Definition at line 81 of file nufamily.h. Referenced by matvec::Population::break_loop(), build_connectors(), cutting(), matvec::Population::detect_loop(), matvec::Population::maxant_maxpost(), matvec::Population::maxant_maxpost_old(), multi_m_terminal_peeling(), multi_terminal_peeling(), nconnector(), and terminal_peeling(). |
|
|
Definition at line 84 of file nufamily.h. Referenced by anterior(), matvec::Population::build_nufamily(), cutting(), multi_ant_post(), multi_anterior(), multi_initialize(), multi_m_ant_post(), multi_m_anterior(), multi_m_initialize(), multi_m_posterior(), multi_m_terminal_peeling(), multi_posterior(), multi_terminal_peeling(), posterior(), and terminal_peeling(). |
|
|
Definition at line 89 of file nufamily.h. Referenced by matvec::Population::break_loop(), matvec::Population::maxant_maxpost(), matvec::Population::maxant_maxpost_old(), NuFamily(), and matvec::Population::peeling_sequence(). |
|
|
Definition at line 79 of file nufamily.h. Referenced by matvec::Population::detect_loop(), log_likelihood(), multi_llh(), multi_m_log_likelihood(), matvec::Population::multi_m_log_likelihood_peeling(), matvec::Population::multipoint_likelihood(), and NuFamily(). |
|
|
Definition at line 46 of file nufamily.cpp. Referenced by multi_m_anterior(), multi_m_log_likelihood(), multi_m_posterior(), and multi_sumint_offspring(). |
|
|
Definition at line 84 of file nufamily.h. Referenced by anterior(), matvec::Population::build_nufamily(), cutting(), multi_ant_post(), multi_anterior(), multi_initialize(), multi_m_ant_post(), multi_m_anterior(), multi_m_initialize(), multi_m_posterior(), multi_m_terminal_peeling(), multi_posterior(), multi_terminal_peeling(), posterior(), and terminal_peeling(). |
|
|
Definition at line 80 of file nufamily.h. Referenced by anterior(), build_connectors(), cutting(), eliminateGenotypes(), father_gid(), father_id(), iterative_peeling(), log_likelihood(), multi_ant_post(), multi_anterior(), multi_fullsibs_prob(), multi_initialize(), multi_llh(), multi_m_ant_post(), multi_m_anterior(), multi_m_initialize(), multi_m_log_likelihood(), multi_m_posterior(), multi_m_terminal_peeling(), multi_posterior(), multi_sumint_offspring(), multi_terminal_peeling(), NuFamily(), pretend_missing(), pretend_multi_m_missing(), pretend_multi_missing(), release(), and terminal_peeling(). |
|
|
Definition at line 43 of file nufamily.cpp. Referenced by multi_m_anterior(). |
|
|
Definition at line 80 of file nufamily.h. Referenced by anterior(), build_connectors(), cutting(), eliminateGenotypes(), log_likelihood(), mother_gid(), mother_id(), multi_ant_post(), multi_anterior(), multi_fullsibs_prob(), multi_initialize(), multi_llh(), multi_m_ant_post(), multi_m_anterior(), multi_m_initialize(), multi_m_log_likelihood(), multi_m_terminal_peeling(), multi_sumint_offspring(), multi_terminal_peeling(), NuFamily(), posterior(), pretend_missing(), pretend_multi_m_missing(), pretend_multi_missing(), release(), and terminal_peeling(). |
|
|
Definition at line 44 of file nufamily.cpp. Referenced by multi_m_anterior(). |
|
|
Definition at line 80 of file nufamily.h. Referenced by build_connectors(), cutting(), eliminateGenotypes(), fullsibs_prob(), iterative_peeling(), multi_ant_post(), multi_fullsibs_prob(), multi_initialize(), multi_m_ant_post(), multi_m_anterior(), multi_m_initialize(), multi_sumint_offspring(), NuFamily(), offspring(), pretend_missing(), pretend_multi_m_missing(), pretend_multi_missing(), and release(). |
|
|
Definition at line 79 of file nufamily.h. |
|
|
Definition at line 79 of file nufamily.h. Referenced by anterior(), build_connectors(), cutting(), eliminateGenotypes(), fullsibs_prob(), iterative_peeling(), multi_ant_post(), multi_anterior(), multi_fullsibs_prob(), multi_initialize(), multi_m_ant_post(), multi_m_anterior(), multi_m_initialize(), multi_m_posterior(), noffs(), NuFamily(), offspring(), pretend_missing(), pretend_multi_m_missing(), pretend_multi_missing(), and release(). |
|
|
Definition at line 89 of file nufamily.h. Referenced by matvec::Population::break_loop(), matvec::Population::maxant_maxpost(), matvec::Population::maxant_maxpost_old(), NuFamily(), and matvec::Population::peeling_sequence(). |
|
|
Definition at line 49 of file nufamily.cpp. Referenced by multi_anterior(), multi_fullsibs_prob(), and pretend_multi_missing(). |
|
|
Definition at line 40 of file nufamily.cpp. Referenced by multi_m_anterior(), and multi_m_posterior(). |
|
|
Definition at line 90 of file nufamily.h. Referenced by matvec::Population::build_nufamily(), cutting(), display(), get_m_posterior(), iterative_peeling(), multi_ant_post(), multi_anterior(), multi_fullsibs_prob(), multi_get_tr(), multi_llh(), multi_m_ant_post(), multi_m_anterior(), multi_m_log_likelihood(), multi_m_posterior(), multi_posterior(), multi_sumint_offspring(), NuFamily(), and release(). |
|
|
Definition at line 88 of file nufamily.h. Referenced by matvec::Population::break_loop(), matvec::compare_seq_id(), matvec::Population::detect_loop(), NuFamily(), and matvec::Population::peeling_sequence(). |
|
|
Definition at line 42 of file nufamily.cpp. Referenced by multi_m_posterior(). |
|
|
Definition at line 79 of file nufamily.h. Referenced by matvec::Population::detect_loop(), and NuFamily(). |
|
|
Definition at line 88 of file nufamily.h. Referenced by anterior(), cutting(), fullsibs_prob(), get_posterior(), log_likelihood(), NuFamily(), matvec::Population::peeling_sequence(), and posterior(). |
|
|
Definition at line 88 of file nufamily.h. Referenced by cutting(), multi_m_anterior(), multi_m_log_likelihood(), multi_m_posterior(), multi_sumint_offspring(), NuFamily(), matvec::Population::peeling_sequence(), and pretend_multi_m_missing(). |
|
|
Definition at line 47 of file nufamily.cpp. Referenced by multi_get_tr(), and multi_sumint_offspring(). |
|
|
Definition at line 45 of file nufamily.cpp. |
|
|
Definition at line 41 of file nufamily.cpp. Referenced by multi_sumint_offspring(). |
|
|
Definition at line 87 of file nufamily.h. Referenced by matvec::Population::break_loop(), matvec::Population::maxant_maxpost(), and matvec::Population::maxant_maxpost_old(). |
|
|
Definition at line 48 of file nufamily.cpp. Referenced by multi_anterior(), multi_fullsibs_prob(), and multi_posterior(). |
1.2.16