//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Bayesian Inference // // Filename: fa.c (was bayesapp.c) // // edited by DJCM Thu 1/4/04 to work with mainfa.c and implement simple // one-latent-variable factor analysis // // the one loading vector is called V. // The prob of data given V is a correlated normal along the lines of // covariance matrix ( I + VV^T ) // The data themselves are summarised by ~ k^2 sufficient statistics, where // k is the number of dimensions. // // These are only example programs, involving inefficient calculations that // repeatedly build an object from scratch with "UserBuild" instructions. // //============================================================================= // History: JS 2 Jan 2002, 4 Jan 2003, 10 Feb 2003, 20 Aug 2003 // Copyright (c) 2002,2003 Maximum Entropy Data Consultants Ltd. //----------------------------------------------------------------------------- #include // stdio added by DJCM #include #include "bayesys3.h" // userstr.h added by DJCM #include "userstr.h" static int UserBuild(double*, CommonStr*, ObjectStr*, int, int); //============================================================================= // SIMPLE CODE //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserEmpty // // Purpose: Set Lhood = logLikelihood(coords) for empty sample // with 0 atoms in Object, and initialise its other information. // // Lhood := log L(empty) // // I have already put the number 0 in Object->Natoms, // so you can use a simple "UserBuild" instruction. //----------------------------------------------------------------------------- // int UserEmpty( // O >=0, or -ve error code double* Lhood, // O loglikelihood CommonStr* Common, // I O general information ObjectStr* Object) // I O sample object, output new Lhood { UserBuild(Lhood, Common, Object, 0, 1); return 1; } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserTry1 // // Purpose: d(logLikelihood(coords)) after supposedly adding one new atom // to Object, without any side effects on it. // // If Valency = 0, // dLtry := d logL(...,x) // else mimic the above with cool // dLtry := (1/cool) log INTEGRAL dPrior(z) dlogL(...,x,z) // // I have already put the new atom x after the last location // of the Atoms list, at Object->Cubes[ Object->Natoms ], so // if Valency=0 you can use simple "UserBuild" instructions. //----------------------------------------------------------------------------- // int UserTry1( // O +ve = OK, 0 = DO NOT USE, -ve = error double* dLtry, // O trial d(logLikelihood) value CommonStr* Common, // I general information ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood) { int Natoms = Object->Natoms; double Lhood = Object->Lhood; // existing Lhood double Ltry; int OK; *dLtry = 0.0; OK = UserBuild(&Ltry, Common, Object, Natoms+1, 0); // trial if( OK > 0 ) *dLtry = Ltry - Lhood; // increment return OK; // OK? } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserTry2 // // Purpose: d(logLikelihood(coords)) after supposedly adding one, // then two new atoms to Object, without any side effects on it. // // If Valency = 0, // dLtry1 := d logL(...,x1) // // dLtry2 := d logL(...,x1,x2) // // else mimic the above with cool // dLtry1 := (1/cool) log INTEGRAL dPrior(z1) dlogL(...,x1,z1) // cool // dLtry2 := (1/cool) log INTEGRAL dPrior(z1,z2) dlogL(...,x1,z1,x2,z2) // // I have already put the new atoms x1,x2 after the last location // of the Atoms list, at Object->Cubes[ Object->Natoms ] and // at Object->Cubes[ Object->Natoms + 1 ] // so if Valency=0 you can use simple "UserBuild" instructions. //----------------------------------------------------------------------------- // int UserTry2( // O +ve = OK, 0 = DO NOT USE, -ve = error double* dLtry1, // O trial d(logLikelihood) value for 1st atom double* dLtry2, // O trial d(logLikelihood) value for both atoms CommonStr* Common, // I general information ObjectStr* Object) // I sample object (DO NOT UPDATE Lhood) { int Natoms = Object->Natoms; double Lhood = Object->Lhood; // existing Lhood double Ltry1, Ltry2; int OK; *dLtry1 = *dLtry2 = 0.0; OK = UserBuild(&Ltry1, Common, Object, Natoms+1, 0); //trial for 1 more if( OK > 0 ) { *dLtry1 = Ltry1 - Lhood; // increment for 1 OK = UserBuild(&Ltry2, Common, Object, Natoms+2, 0); //trial for 2 more if( OK > 0 ) *dLtry2 = Ltry2 - Lhood; // increment for 2 } return OK; // return OK only if both trials OK } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserInsert1 // // Purpose: Insert 1 new atom into Object, keeping it up-to-date, and // set d(loglikelihood(coords)). // // If Valency = 0, // dL := d logL(...,x) // // else // cool // sample z from Prior(z) L(...,x,z) // // and set dL := d logL(...,x,z) at the sampled z. // // I have already put the new atom x at the last location of // the updated Atoms list, at Object->Cubes[ Object->Natoms - 1 ], //----------------------------------------------------------------------------- // int UserInsert1( // O >=0, or -ve error code double* dL, // O d(loglikelihood) CommonStr* Common, // I general information ObjectStr* Object) // I O sample object { int Natoms = Object->Natoms; // new number double Lold = Object->Lhood; // not yet updated double Lnew; UserBuild(&Lnew, Common, Object, Natoms, 1); // new updated state *dL = Lnew - Lold; // I will update Object->Lhood return 1; } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserInsert2 // // Purpose: Insert 2 new atoms into Object, keeping it up-to-date, and // set d(loglikelihood(fluxes)). // // If Valency = 0, // dL := d logL(...,x1,x2) // // else // cool // sample z1,z2 from Prior(z1,z2) L(...,x1,z1,x2,z2) // // and set dL := d logL(...,x1,z1,x2,z2) at the sampled z1,z2. // // I have already put the new atoms x1,x2 at the last location // of the Atoms list, at Object->Cubes[ Object->Natoms - 2 ] and // at Object->Cubes[ Object->Natoms - 1 ] // so if Valency=0 you can use a simple "UserBuild" instruction. //----------------------------------------------------------------------------- // int UserInsert2( // O >=0, or -ve error code double* dL, // O d(loglikelihood) CommonStr* Common, // I general information ObjectStr* Object) // I O sample object { int Natoms = Object->Natoms; // new number double Lold = Object->Lhood; // not yet updated double Lnew; UserBuild(&Lnew, Common, Object, Natoms, 1); // new updated state *dL = Lnew - Lold; // I will update Object->Lhood return 1; } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserDelete1 // // Purpose: Delete 1 old atom from Object, keeping it up-to-date, and // set d(loglikelihood(fluxes)). // // dL := d logL(...) // // I have already put the old atom after the last location of // the updated Atoms list, at Object->Cubes[ Object->Natoms ], // so if Valency=0 you can use a simple "UserBuild" instruction. //----------------------------------------------------------------------------- // int UserDelete1( // O >=0, or -ve error code double* dL, // O d(loglikelihood) CommonStr* Common, // I general information ObjectStr* Object) // I O sample object { int Natoms = Object->Natoms; // new number double Lold = Object->Lhood; // not yet updated double Lnew; UserBuild(&Lnew, Common, Object, Natoms, 1); // new updated state *dL = Lnew - Lold; // I will update Object->Lhood return 1; } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Function: UserBuild // // Purpose: Build Lhood, and any other data info such as Mock from scratch. // Can be called with any number of atoms // // Thu 1/4/04 DJCM uses one atom with Ndim dimension to get normal variables. // // Prior(V) = Normal( 0 , 1.0 ) // //----------------------------------------------------------------------------- // static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error double* Lhood, // O loglikelihood CommonStr* Common, // I General information ObjectStr* Object, // I(O) Cubes in, perhaps Mock out int Natoms, // I # atoms int output) // I +ve if Mock to be written out, -ve if not { const double Sqrt2Pi = 2.50662827463100050240; double** Cubes = Object->Cubes; // I Cubes in [0,1) [Natoms][Ndim] int Ndata = Common->Ndata; // I # data double* Data = Common->Data; // I Data [Ndata] Contain covariance matrix sufficient stats double* Acc = Common->Acc; // I Accuracies [Ndata] Not used. #define KDIMENSION 2 // same as Ndim double V[KDIMENSION]; double My_f ; // double A[KDIMENSION][KDIMENSION]; double Lnorm = 0.0; // Likelihood normalisation double C = 0.0; // chisquared double flux; double* Cube; int k, k2 , kkk ; UserCommonStr* UserCommon = (UserCommonStr*)Common->UserCommon; double TRUE_N_data = UserCommon->TRUE_N_data ; double VVAR = UserCommon->VVAR ; double DVAR = UserCommon->DVAR ; if (!(Natoms==1)) { // fprintf( stderr , "WARNING, Natoms %d not 1!\n" , Natoms ) ; *Lhood = 0.0 ; return 1 ; } // Recover the loading vector, V. There should only be one atom. for( k = 0; k < Natoms; k++ ) { Cube = Cubes[k]; BayeShape( V , Cube , KDIMENSION , 5 ); } // Compute the factor appearing in all the formulae for( My_f = 1.0 , k = 0; k < KDIMENSION; k++ ) My_f += VVAR * V[k]/DVAR * V[k] ; // Accumulate logLikelihood - use explicit formula based on sufficient statistics for( kkk=0 , k = 0; k < KDIMENSION; k++ ) for( k2 = 0 ; k2 < KDIMENSION ; k2++ ) { C += (k==k2?Data[kkk]/DVAR:0.0) - VVAR * V[k]/DVAR * V[k2]/DVAR * Data[kkk] / My_f ; kkk ++ ; } Lnorm += TRUE_N_data * log( 1.0 / Sqrt2Pi); // (fussy to include this) Lnorm -= TRUE_N_data * 0.5 * log( My_f + KDIMENSION*log(DVAR) ) ; *Lhood = Lnorm - C / 2.0; if( output ) for( kkk=0 , k = 0; k < KDIMENSION; k++ ) for( k2 = 0; k2 < KDIMENSION ; k2++ ) { Object->Mock[kkk] = TRUE_N_data * ( VVAR * V[k] * V[k2] + ( k2==k ? DVAR : 0.0 ) ) ; kkk ++ ; } return 1; } //============================================================================= // Dummy procedure to link to BayeSys3 when not running MassInf //============================================================================= int UserFoot( double* Cube, CommonStr* Common, int* ibits, double* zbits, int* nbits) { return (Cube && Common && ibits && zbits && nbits) ? 1 : 1; }