Skip to content

Commit

Permalink
updated updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@lp2
lp2 committed Nov 27, 2010
1 parent 90efd49 commit 29db6c5
Show file tree
Hide file tree
Showing 5 changed files with 81 additions and 38 deletions.
20 changes: 12 additions & 8 deletions bayesnet.cpp
View file
Expand Up @@ -18,8 +18,8 @@ using alglib::psi;


double cNode::entropy() {return 0;}
double cNode::calcBound(cBayesNet *net) {return 0;}
void cNode::update(cBayesNet *net) {}
double cNode::calcBound(cBayesNet &net) {return 0;}
void cNode::update(cBayesNet &net) {}


cDirichletNode::cDirichletNode(int dim, float prior_val){
Expand All @@ -42,7 +42,7 @@ double cDirichletNode::entropy(){
}


double cDirichletNode::calcBound(cBayesNet *net){
double cDirichletNode::calcBound(cBayesNet &net){
double temp = 0;
double bound = ((a0 - VectorXf::Ones(a0.rows()))*lnE).sum() - lngamma(a0.sum(), temp);
for (int i=0; i < a0.rows(); i++) { bound += lngamma(a0(i), temp); }
Expand All @@ -51,7 +51,7 @@ double cDirichletNode::calcBound(cBayesNet *net){


/** Not used, so not implemented */
void cDirichletNode::update(cBayesNet *net){}
void cDirichletNode::update(cBayesNet &net){}



Expand All @@ -67,7 +67,7 @@ cGammaNode::cGammaNode(int dim, float prior_val_a, float prior_val_b, VectorXf *
b = pb*VectorXf::Ones(dim);
E1 = VectorXf::Zero(dim);
lnE = VectorXf::Zero(dim);
update((cBayesNet*)NULL);
updateParams();
if (E1_val != NULL){
E1 = *E1_val;
}
Expand All @@ -90,7 +90,7 @@ double cGammaNode::entropy(){



double cGammaNode::calcBound(cBayesNet *net){
double cGammaNode::calcBound(cBayesNet &net){

double bound = 0;
double temp = 0;
Expand All @@ -104,14 +104,18 @@ double cGammaNode::calcBound(cBayesNet *net){
}


/** Standard update equation */
void cGammaNode::update(cBayesNet *net){
void cGammaNode::updateParams(){
for (int i=0; i < lnE.rows(); i++){
E1(i) = a(i)/b(i);
lnE(i) = psi(a(i)) - log(b(i));
}
}

/** Standard update equation */
void cGammaNode::update(cBayesNet &net){
this->updateParams();
}



/*
Expand Down
13 changes: 7 additions & 6 deletions bayesnet.h
View file
Expand Up @@ -32,8 +32,8 @@ class cNode {
double bound;

virtual double entropy();
virtual double calcBound(cBayesNet *net);
virtual void update(cBayesNet *net);
virtual double calcBound(cBayesNet &net);
virtual void update(cBayesNet &net);
};


Expand All @@ -49,8 +49,8 @@ class cDirichletNode : public cNode {
public:
cDirichletNode(int dim, float prior);
double entropy();
double calcBound(cBayesNet *net);
void update(cBayesNet *net);
double calcBound(cBayesNet &net);
void update(cBayesNet &net);
};


Expand All @@ -69,8 +69,9 @@ class cGammaNode : public cNode {
cGammaNode();
cGammaNode(int dim, float prior_val_a, float prior_val_b, VectorXf *E1_val);
double entropy();
double calcBound(cBayesNet *net);
void update(cBayesNet *net);
double calcBound(cBayesNet &net);
void update(cBayesNet &net);
void updateParams();
};


Expand Down
31 changes: 28 additions & 3 deletions main.cpp
View file
Expand Up @@ -38,6 +38,7 @@ void playing()
}



MatrixXf randn(int n, int m)
/* create a randn matrix */
{
Expand All @@ -49,6 +50,30 @@ MatrixXf randn(int n, int m)
return rv;
}

/** Trying out update equations etc */
void play_matrix(){
int N = 3;
int P = 4;
int K = 2;

MatrixXf alpha = MatrixXf::Zero(K,K);//randn(K,1);
alpha.diagonal() = randn(K,1);
MatrixXf XE2s = randn(K,K);
MatrixXf X = randn(N,K);
MatrixXf eps = randn(P,1);
MatrixXf E1 = randn(P,K);
MatrixXf E2S = randn(K,K);
MatrixXf pheno = randn(N,P);

for(int i=0; i<P; ++i){
MatrixXf prec = alpha + XE2s*eps(i);
MatrixXf cov = prec.inverse();
}
}




MatrixXf simulate_expression(int N, int D, int K,double sigma=0.1)
/*
Simulate an expresison matrix with N entries, D dimensions and K factors
Expand All @@ -69,11 +94,11 @@ MatrixXf simulate_expression(int N, int D, int K,double sigma=0.1)


int main (int argc, char * const argv[]) {

play_matrix();
//1. simulate small dataset
MatrixXf Y=simulate_expression(10,100,5);
// MatrixXf Y=simulate_expression(10,100,5);
//2. create object
cVBFA vb(Y);
// cVBFA vb(Y);


}
21 changes: 15 additions & 6 deletions vbfa.cpp
View file
Expand Up @@ -25,7 +25,7 @@ cXNode::cXNode(MatrixXf E1){
}


double cWNode::calcBound(cBayesNet *net){
double cWNode::calcBound(cBayesNet &net){
return 0.;
}

Expand All @@ -35,12 +35,21 @@ double cWNode::entropy(){
}


void cWNode::update(cBayesNet *net){
void cWNode::update(cBayesNet &net){
cVBFA n = (cVBFA&)net;

// // for each phenotype, calculate covariance and mean of weight
for(int i = 0; i < n.Np; i++){
// MatrixXf diagAE1 = MatrixXf::Zero(n.
// cov = (diagAE1 + n.X.E2S*n.Eps.E1(i)).inverse(); // linalg.inv(diag(Alpha.E1) + Eps[d]*M)
// this->E1(d) = cov*(net->Eps.E1(i)&net->S.E1&net->pheno.E1(i)); // self.E1[d,:] = S.dot(dcov[:,:],Eps[d]*S.dot(_S.E1.T,net.dataNode.E1[ :,d]))
// this->E2S += (cov + this->E1(d)*this->E1(d)); // E2 = dcov + outer(self.E1[d], self.E1[d])
}
}



double cXNode::calcBound(cBayesNet *net){
double cXNode::calcBound(cBayesNet &net){
return 0.;
}

Expand All @@ -50,15 +59,15 @@ double cXNode::entropy(){
}


void cXNode::update(cBayesNet *net){
void cXNode::update(cBayesNet &net){
}


void cAlphaNode::update(cBayesNet *net){
void cAlphaNode::update(cBayesNet &net){
}


void cEpsNode::update(cBayesNet *net){
void cEpsNode::update(cBayesNet &net){
}


Expand Down
34 changes: 19 additions & 15 deletions vbfa.h
View file
Expand Up @@ -12,41 +12,45 @@

#include "bayesnet.h"





/** Weight node (vector normal with shared precision) */
class cWNode : public cNode {
protected: MatrixXf E2;
protected: MatrixXf E2S;
public:
cWNode(); // default
cWNode(MatrixXf E1); // from mean
void update(cBayesNet *net);
double calcBound(cBayesNet *net);
void update(cBayesNet &net);
double calcBound(cBayesNet &net);
double entropy();
};


/** Factor node (vector normal) */
class cXNode : public cNode {
protected: MatrixXf E2;
protected: MatrixXf E2S;
public:
cXNode(); // default
cXNode(MatrixXf E1); // from mean
void update(cBayesNet *net);
double calcBound(cBayesNet *net);
void update(cBayesNet &net);
double calcBound(cBayesNet &net);
double entropy();
};


/** Weight precision node (vector Gamma) */
class cAlphaNode : public cGammaNode {
public:
void update(cBayesNet *net);
void update(cBayesNet &net);
};


/** Noise precision node (vector Gamma) */
class cEpsNode : public cGammaNode {
public:
void update(cBayesNet *net);
void update(cBayesNet &net);
};


Expand All @@ -63,8 +67,11 @@ protected: MatrixXf E1,E2;

/** Main class for variational Bayesian factor analysis */
class cVBFA : public cBayesNet {

protected:


public:
/** Dimensions */
int Nj; //individuals
int Np; //phenotypes
Expand All @@ -77,14 +84,10 @@ class cVBFA : public cBayesNet {
cAlphaNode Alpha;
cPhenoNode pheno;




public:
/** Inference parameters */
int Niterations;
int initialisation; // type

//default constructor
//cVBFA();

Expand All @@ -93,10 +96,11 @@ class cVBFA : public cBayesNet {
//constructor that take variance into account
cVBFA(MatrixXf pheno_mean,MatrixXf pheno_var);


void update();


};


#endif

0 comments on commit 29db6c5

Please sign in to comment.