RandomVar.h

Go to the documentation of this file.

// -*- C++ -*-


// PLearn (A C++ Machine Learning Library)
// Copyright (C) 1998 Pascal Vincent
// Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
//

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// 
//  1. Redistributions of source code must retain the above copyright
//     notice, this list of conditions and the following disclaimer.
// 
//  2. Redistributions in binary form must reproduce the above copyright
//     notice, this list of conditions and the following disclaimer in the
//     documentation and/or other materials provided with the distribution.
// 
//  3. The name of the authors may not be used to endorse or promote
//     products derived from this software without specific prior written
//     permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
// NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
// 
// This file is part of the PLearn library. For more information on the PLearn
// library, go to the PLearn Web site at www.plearn.org


 

/* *******************************************************      
   * $Id: RandomVar.h,v 1.6 2004/07/21 16:30:54 chrish42 Exp $
   * AUTHORS: Pascal Vincent & Yoshua Bengio
   * This file is part of the PLearn library.
   ******************************************************* */

/* RandomVar.h

   Random Variables package:

   RandomVar class, helper classes such as RandomVariable,
   RVInstance, RVArray, RVInstanceArray, ConditionalExpression,
   and many subclasses of RandomVariable,
   plus helper global functions and operators.

 */

/* TUTORIAL ON THE RANDOMVAR PACKAGE

   Whereas a Var represents a variable in the mathematical sense, 
   a RandomVar represents a random variable in the mathematical sense. 
   A random variable is generally defined in terms of other random variables
   through deterministic transformations, or in terms of other random
   variables which are the parameters of a its distribution.
   A RandomVar represents a node in a graphical
   model, and its distribution is defined in terms
   of the values of its parents in the model.
   For example, its parents may be the parameter
   of its distribution or it may be the variables
   which when combined deterministically give rise
   to its value. The set of classes provided here
   allow to build such a network of random variables,
   and to make limited inferences, probability computations,
   gradient computations, and learning.

   Examples of use of RandomVars:

   Var u(1),lv(1); // constants, but changing their value will change
   Var w; // the distributions associated to the RandomVar's
   ...
   // things like *, +, log, tanh, etc... must do the proper thing for RandomVars
   RandomVar X = gamma(0.5)*log(normal(u,lv));
   RandomVar Y = tanh(product(w,X) + u);
   RandomVar LW(2); // unnormalized log-weights of the mixture
   RandomVar Z = mixture(Y&X,LW); 

   // see the comment on operator&, operator[] and mixtures below
   ...
   // (conditionned on the RHS of the |)
   Vec x,y; // put some value in x and y
   Vec z = sample(Z|(X==x && Y==y));
   // This is achieved by redefining operator==, operator&& and operator|
   // to represent the data structure which the function sample
   // takes as argument. In particular, note that == creates the
   // RVInstance data structure (which contains a RV and a Var instance),
   // while && or & makes an array of these structures (RVInstanceArray), and
   // the "|",  builds a ConditionalExpression, that contains a RVInstanceArray
   // for the left hand side and a RVInstanceArray for the right hand side. Note
   // the use of parentheses because of the default precedence of operators.
   ...
   // Using these operators, you can also express the computation of the probability of
   // a value for the random variable, or its conditional probability. In fact the
   // statement below defines a functional relationship (in the usual "Variable"
   // sense between the variable y and the variable p or log_p). Note that
   // no actual numerical value has yet been computed. 
   Var y(1);
   Var p = P(Y==y);
   Var log_p = logP(Y==y);
   // or
   Var y,x;
   Var log_p = logP((Y==y)|(X==x)); // note again ()'s for precedence
   // e.g. of use
   Vec actual_value_for_x_and_y = ...;
   Func f = log_p(x&y);
   real prob = f(actual_value_for_x_and_y)[0];
   cout << "log(P(Y|X))=" << prob << endl;
   ...
   // in the case of discrete distributions, the whole distribution can be returned
   // by P and logP, which are also defined on RandomVars.
   Vec p = P(Y);
   Vec log_p = logP(Y);
   ...
   // note that if Y has RandomVar parents, the above var can only be computed
   // if these parents are given a particular value or if
   // they are integrated over (with the function marginalize below).
   // This call will therefore automatically try to marginalize Y
   // (by integrating over the parents which are not observed).
   // To make a RandomVar observed, simply use the conditioning
   // notation V|(X==x && Y==y && Z==z), e.g. to condition V on X=x,Y=y,Z=z.
   ...
   // Similarly to P, logP, and sample, other functions of RVs are defined:
   // construct a Var that is functionally dependent on the Var x
   // and represents the expectation of Y given that X==x.
   Var e=E(Y|(X==x));
   // Similarly for covariance matrix:
   Var v=V(Y|(X==x && Z==z));
   // and for the cumulative distribution function
   // (which here depends on the Vars y and x)
   Var c=P((Y<y)|(X==x));
   // Note that derivatives through all these functional relationships
   // can automatically be computed. For example, to compute the
   // gradients of the log-probability wrt some parameters W & B & LogVariance
   RandomVar W,B,X,LogVariance; // all are "non-random", X is the input
   Var w,b,x,lv; // values that the above will take
   // e.g. to give values to these Vars
   w[0]=1; b[0]=0;x[0]=3;lv[0]=0;
   RandomVar Y=normal(W*X+B,LV); // the model (i.e. a regression)
   // establish the functional relationship of interest,
   // which goes from (y & x & w & b & lv) to logp:
   Var logp = logP((Y==y)|(X==x && W==w && B==b && LV==lv));
   // to actually compute logp, do
   logp->fprop_from_all_sources(); // a source is "constant" Var
   // to compute gradients, use the propagationPath function to find
   // the path of Vars from, say w&b, to logp:
   VarArray prop_path = propagationPath(w&b,logp);
   prop_path.bprop(); // compute dlogP/dparams in w->gradient and b->gradient
   ...
   // By default a RandomVar represents a "non-random" random variable 
   // (of class NonRandomVariable or a FunctionalRandomVariable which depends
   // only on NonRandomVariable). This is not the same as a "constant"
   // variable. It only means that its value is deterministic, but
   // its value may be a (Var) function of other Vars:
   RandomVar X;
   X->value = 1+exp(y+w*z); // y, w and z are Vars here
   ...
   // Parameters of the distributions that define the random variables can be
   // learned, most generally by optimizing a criterion to optimize, e.g.
   Mat observed_xy_matrix; // each row = concatenation of an X and a Y obs.
   Var cost =  -logP((Y==y)|(X==x && Params==params));
   // below, establish the functional relationship between x & y & params
   // and the "totalcost" which is the sum of "cost" when x & y are
   // "sampled" from the given empirical distribution.
   Var totalcost = meanOf(cost,x&y,VMat(observed_xy_matrix),
                          observed_xy_matrix.length(),
                          params);
   // construct an optimizer that optimizes totalcost by varying params
   ConjugateGradientOptimizer opt(params,totalcost, 1e-5, 1e-3, 100);
   // do the actual optimization
   real train_NLL = opt.optimize();
   // now we can test on a particular point setting values for
   // x and y and params and doing an fprop with
   propagationPath(x&y&params,cost).fprop();
   ...
   // Sometimes, the parameters can be estimated more efficiently
   // using an internal mechanism for estimation (usually the analytical 
   // solution of maximum likelihood or the EM algorithm):
   real avgNegLogLik=EM(Y|X,W&B&LV,VMat(observed_xy_matrix),
                         observed_xy_matrix.length(),4,0.001);
   // where the first argument specifies which conditional distribution
   // is of interest (ignoring the parameters), the second argument
   // gives the parameters to estimate, and the third one specifies
   // a training set. Note that the order of the variables in the 
   // observed_xy_matrix must be (1) inputs: all the variables on the RHS 
   // of the conditioning |, (2) outputs: all the variables on the LHS of the |.
   ...
   // Arrays of RVs can be formed with the operator &:
   RVArray a = X & Y & V;
   // and they can be automatically cast into JointRandomVariables
   // (whose value is the concatenation of the values of its parents)
   RandomVar Z = X & Y & V;
   ...
   // Marginals can sometimes be obtained (when X is discrete or 
   // the integral is feasible analytically and the code knows how to do it...).
   // For example, suppose X is one of the parents of Y. Then
   RandomVar mY = marginalize(Y,X);
   // is a random variable such that P(mY=y) = int_x P(Y=y|X=x)P(X=x) dx
   // this is obtained by summing over the values of X if it is discrete,
   // by doing the integral if we know how to do it, or otherwise, by the
   // Laplace approximation, or by some numerical integration method
   // such as Monte-Carlo.
   ...
   // The operator() is defined on RandomVar as follows:
   // If i is an integer, X(i) extracts a RandomVar that is scalar 
   // and corresponds to the i-th element of the vector random variable X 
   // (similarly, if the underlying Var is a Var, X(i,j), 
   // extracts the random element (i,j)). These two operators
   // are also defined for the case in which the index is a Var
   // (treated like the integer), and the case in which it is a RandomVar.
   // The last case, X(I), actually represents a mixture of the elements
   // of the vector X, with weights given by the parameters of I
   // (which must be discrete).
   ...
   // The operator[] is defined on RVArrays
   // and it allows to extract the i-th random variable in the array.
   // With I a RandomVar and A an RVArray, A[I], is very interesting because 
   // it represents the graphical model of a mixture in which I is the index,
   // and it is not (yet) integrated over. 
   RVArray A(3);
   A[0]=X; A[1]=Y; A[2]=Z;
   RandomVar XYZ(A); // joint distribution
   // or equivalently
   RandomVar XYZ = X & Y & Z;
   ...
   // A MultinomialRandomVariable is a subclass of RandomVariable
   // that represents discrete-valued variables in the range 0 ... N-1.
   RandomVar LW(3); // unnormalized log-probabilities of I
   RandomVar I = multinomial(LW); // N=3 possible values here
   // The parameters of a discrete random variable are the "log-probabilities" 
   // (more precisely the discrete probabilities are obtained with a softmax 
   // of the parameters, LW here). The discrete random variable will be 
   // conditional if LW is not a NonRandomVariable but rather depends
   // on some other RVs.
   ...
   // Let us consider a random variable that is obtained by selecting
   // one of several random variables (on the same space). We call
   // such a random variable an RVArrayRandomElementRandomVariable and it is 
   // obtained with the operator[] acting on a RVArray, 
   // with a discrete randomvariable as argument:
   RandomVar V = A[I]; // will take either the distribution of X, Y or Z according to value of I
   // therefore P(V=v|I=i)=P(A[i]=v)
   ...
   // A mixture is the marginalization of an IndexedRandomVariable with
   // respect to the random index:
   RandomVar LW(3); // unnormalized log-weights of the mixture
   RandomVar M = mixture(a,LW);
   // which is exactly the same thing as
   RandomVar I = multinomial(LW); 
   RandomVar M = marginalize(A[I],I);
   ...
   // Example of conditional mixture of n d-dimensional diagonal 
   // Gaussians with neural network expectations:
   // (1) define the neural network
   RandomVar X(n_inputs); // X is the network input RV
   Var x(n_inputs); // x will be its value
   int n_inputs=4, n_hidden=5,n_outputs=d*n;
   Var layer1W(n_hidden,n_inputs), layer2W(n_outputs,n_hidden);
   Var layer1bias(n_hidden), layer2bias(n_outputs);
   RandomVar NetOutput = layer2bias+layer2W*tanh(layer1bias+layer1W*X);
   // (2) define the gaussian mixture
   // (2.1) define the gaussians
   RVArray normals(n);
   RVArray mu(n),logsigma(n);
   for (int i=0;i<n;i++) {
     mu[i]=NetOutput.subVec(i*d,d); // extract subvector as i-th mean vector
     normals[i]=Normal(mu[i],logsigma[i]);
   }
   // (2.2) build the mixture itself
   Var lw(n);
   lw->value.fill(1.0/n);
   RandomVar Y = mixture(normals,lw); //the "target" output random variable
   Var y(d); // its value
   VarArray tunable_parameters = 
      lw & layer1W & layer2W & layer1bias & layer2bias;
   // each row of Mat is the concatenation of an x (input) and a y 
   Mat observed_xy_matrix; 
   // logP returns the path that goes from all "source" (constant) variables
   // into the computation of the given conditional probability
   Var cost =  -logP((Y==y)|(X==x)); 
   // note that we don't need to condition on the "non-random" parameters
   // such as the log-weights of the mixture (lw), but they will
   // occur as tunable parameters.
   // Below, the order of x and y in the observed_xy_matrix must
   // match their order in the second argument of meanOf.
   Var totalcost = meanOf(cost,x&y,VMat(observed_xy_matrix),
                          observed_xy_matrix.length(),tunable_parameters);
   ...
   // Example in which some parameters W & B have to be fitted
   // to some data, while the hyper-parameters gamma that control
   // the distribution of W & B should be fitted to maximize
   // the likelihood of the data.
   int npoints = 10; // there are 10 (x,y) pairs in each observation
   Var muW, logvarW, muB, logvarB; // parameters of the prior
   RandomVar W = normal(muW,logvarW); // prior on W
   RandomVar B = normal(muB,logvarB); // prior on B
   Var log_var(1); // log-variance of Y
   Var ones(npoints);
   ones->value.fill(1.0);
   Var log_vars = ones*log_var; // make vector of npoints copies of log_var
   Var x(npoints); // input
   Var muXint, logvXint; // parameters of Xinterval
   RandomVar Xinterval = normal(muXint,logvXint); // prior on Xinterval
   RandomVar Y = normal(tanh(W*x+B),log_vars); 
   Var cost = -logP(Y==y && Xinterval==vconcat(min(x) & max(x));
   // note that the above requires marginalizing over W & B
   VarArray gamma = muXint & logvXint & log_var & muW & logvarW & muB & logvarB;
   Var totalcost = meanOf(cost,x&y,VMat(observed_xy_matrix),
                          observed_xy_matrix.length(), gamma);
   ConjugateGradientOptimizer opt(gamma,totalcost, 1e-5, 1e-3, 100);
   real train_NLL = opt.optimize();
   // to obtain a fit of Theta for a particular value of x and y, optimize
   Var w,b;
   Var fitcost = 
     -logP((Y==y && Xinterval==vconcat(min(x) & max(x)))|(W==w && B==b));
   ConjugateGradientOptimizer fitopt(w & b,fitcost, 1e-5, 1e-3, 100);
   real fit_NLL = fitopt.optimize();
   // and the fitted parameters can be read in w->value and b->value.
*/


#ifndef RANDOMVAR_INC
#define RANDOMVAR_INC

#include <plearn/opt/Optimizer.h>

#include "SampleVariable.h"

#include <plearn/vmat/VMat.h>

namespace PLearn {
using namespace std;


class RandomVariable;
class RVArray;
class RVInstance;
class RVInstanceArray;
class ConditionalExpression;

class RandomVar: public PP<RandomVariable>
{
 public:
  RandomVar();
  RandomVar(int length, int width=1);
  RandomVar(RandomVariable* v);
  RandomVar(const RandomVar& other);

  RandomVar(const Vec& vec);
  RandomVar(const Mat& mat);
  RandomVar(const Var& var);

  RandomVar(const RVArray& vars); 
 
  RandomVar operator[](RandomVar index); 
#if 0
  RandomVar operator[](int i); 
#endif

  void operator=(const RVArray& vars);

  void operator=(real f); 
  void operator=(const Vec& v); 
  void operator=(const Mat& m);
  void operator=(const Var& v); 

  RVInstance operator==(const Var& v) const;

  bool operator==(const RandomVar& rv) const { return rv.ptr == this->ptr; }
  bool operator!=(const RandomVar& rv) const { return rv.ptr != this->ptr; }

  RVArray operator&(const RandomVar& v) const;

  ConditionalExpression operator|(RVArray rhs) const;

  ConditionalExpression operator|(RVInstanceArray rhs) const;

#if 0
  RandomVar operator[](RandomVar index);
  RandomVar operator[](int i);
  RandomVar operator()(RandomVar i, RandomVar j);
  RandomVar operator()(int i, int j);
#endif

};

typedef RandomVar MatRandomVar;


class RVArray: public Array<RandomVar>
{
 public:
  RVArray();
  RVArray(int n, int n_extra_allocated=0);
  RVArray(const Array<RandomVar>& va);
  RVArray(const RandomVar& v, int n_extra_allocated=0);
  RVArray(const RandomVar& v1, const RandomVar& v2, int n_extra_allocated=0);
  RVArray(const RandomVar& v1, const RandomVar& v2, const RandomVar& v3, 
          int n_extra_allocated=0);

  int length() const;

  VarArray values() const;



  RandomVar operator[](RandomVar index);

  RandomVar& operator[](int i)
    { return Array<RandomVar>::operator[](i); }

  const RandomVar& operator[](int i) const
    { return Array<RandomVar>::operator[](i); }

  static int compareRVnumbers(const RandomVar* v1, const RandomVar* v2);

  void sort();
};


class RVInstance 
{
 public:
  RandomVar V;
  Var v;

  RVInstance(const RandomVar& VV, const Var& vv);
  RVInstance();

  RVInstanceArray operator&&(RVInstance rvi);

  ConditionalExpression operator|(RVInstanceArray a);

  void swap_v_and_Vvalue();

};

class RVInstanceArray: public Array<RVInstance>
{
 public:
  RVInstanceArray();
  RVInstanceArray(int n, int n_extra_allocated=0);
  RVInstanceArray(const Array<RVInstance>& a);
  RVInstanceArray(const RVInstance& v, int n_extra_allocated=0);
  RVInstanceArray(const RVInstance& v1, const RVInstance& v2,
                  int n_extra_allocated=0);
  RVInstanceArray(const RVInstance& v1, const RVInstance& v2, 
                  const RVInstance& v3, int n_extra_allocated=0);

  int length() const; 

  RVInstanceArray operator&&(RVInstance rhs);

  ConditionalExpression operator|(RVInstanceArray rhs);

  RVArray random_variables() const;

  VarArray values() const;
  VarArray instances() const;

  void swap_v_and_Vvalue()
    { for (int i=0;i<size();i++) (*this)[i].swap_v_and_Vvalue(); }

  static int compareRVnumbers(const RVInstance* rvi1, const RVInstance* rvi2);

  void sort();

};

class ConditionalExpression
{
 public:
  RVInstance LHS;
  RVInstanceArray RHS;

  ConditionalExpression(RVInstance lhs, RVInstanceArray rhs);
  ConditionalExpression(RVInstance lhs);
  ConditionalExpression(RandomVar lhs);
  ConditionalExpression(RVInstanceArray lhs);
};

class RandomVariable: public PPointable
{
  friend class RandomVar;
  friend class RVInstanceArray;
  friend class RVArray;

  static int rv_counter; 


protected:

  const int rv_number; 

 public:
  const RVArray parents; 

  Var value; 
  
 protected:
  bool marked;

  bool EMmark; 
  bool pmark; 

  bool* learn_the_parameters;

 public:
  RandomVariable(int thelength, int thewidth=1);
  RandomVariable(const Vec& the_value);
  RandomVariable(const Mat& the_value);
  RandomVariable(const Var& the_value);
  RandomVariable(const RVArray& parents, int thelength);
  RandomVariable(const RVArray& parents, int thelength, int thewidth);
  
  virtual char* classname() = 0;

  virtual int length() { return value->length(); }
  virtual int width() { return value->width(); }
  int nelems() { return value->nelems(); }
  bool isScalar() { return length()==1 && width()==1; }
  bool isVec() { return width()==1 || length()==1; }
  bool isColumnVec() { return width()==1; }
  bool isRowVec() { return length()==1; }

  virtual bool isNonRandom() = 0;
  inline bool isConstant() { return isNonRandom() && value->isConstant(); }

  virtual bool isDiscrete() = 0;

  RandomVar subVec(int start, int length);




  virtual void setValueFromParentsValue() = 0;

  void markRHSandSetKnownValues(const RVInstanceArray& RHS)
    {
      for (int i=0;i<RHS.size();i++)
        RHS[i].V->mark(RHS[i].v);
      setKnownValues();
    }
      
  virtual void EMBprop(const Vec obs, real posterior) = 0;

  virtual void EMUpdate();

  virtual bool canStopEM();

  virtual void EMTrainingInitialize(const RVArray& parameters_to_learn);

  virtual void EMEpochInitialize();


  virtual void mark(Var v) { marked = true; value = v; }
  virtual void mark() { marked = true; }
  virtual void unmark() { marked = false; }
  virtual void clearEMmarks();

  virtual void unmarkAncestors();

  virtual bool isMarked() { return marked; }

  virtual void setKnownValues();

  virtual Var logP(const Var& obs, const RVInstanceArray& RHS,
      RVInstanceArray* parameters_to_learn=0) = 0;
  virtual Var P(const Var& obs, const RVInstanceArray& RHS);

  virtual Var ElogP(const Var& obs, RVArray& parameters_to_learn, 
                    const RVInstanceArray& RHS);




  virtual real EM(const RVArray& parameters_to_learn,
                   VarArray& prop_path, 
                   VarArray& observedVars, 
                   VMat distr, int n_samples, 
                   int max_n_iterations, 
                   real relative_improvement_threshold,
                   bool accept_worsening_likelihood=false);

  virtual real epoch(VarArray& prop_path, 
                      VarArray& observed_vars, const VMat& distr, 
                      int n_samples, 
                      bool do_EM_learning=true);

  virtual ~RandomVariable();

};



RandomVar operator*(RandomVar a, RandomVar b);

RandomVar operator+(RandomVar a, RandomVar b);

RandomVar operator-(RandomVar a, RandomVar b);

RandomVar operator/(RandomVar a, RandomVar b);

RandomVar exp(RandomVar x);

RandomVar log(RandomVar x);

RandomVar extend(RandomVar v, real extension_value = 1.0, int n_extend = 1);

RandomVar hconcat(const RVArray& a);




//real EMbyEMBprop(ConditionalExpression conditional_expression, 
real EM(ConditionalExpression conditional_expression, 
         RVArray parameters_to_learn,
         VMat distr, int n_samples, int max_n_iterations=1, 
         real relative_improvement_threshold=0.001,
         bool accept_worsening_likelihood=false,
         bool compute_final_train_NLL=true);

real oEM(ConditionalExpression conditional_expression,
         RVArray parameters_to_learn,
         VMat distr, int n_samples, int max_n_iterations, 
         real relative_improvement_threshold=0.001,
         bool compute_final_train_NLL=true);

real oEM(ConditionalExpression conditional_expression,
         RVArray parameters_to_learn,
         VMat distr, int n_samples, 
         Optimizer& MStepOptimizer,
         int max_n_iterations,
         real relative_improvement_threshold=0.001,
         bool compute_final_train_NLL=true);

Var logP(ConditionalExpression conditional_expression, 
         bool clearMarksUponReturn=true,
         RVInstanceArray* parameters_to_learn=0);

Var P(ConditionalExpression conditional_expression,
         bool clearMarksUponReturn=true);

Var ElogP(ConditionalExpression conditional_expression, 
    RVInstanceArray& parameters_to_learn,
    bool clearMarksUponReturn=true);

RandomVar marginalize(const RandomVar& RV, const RandomVar& hiddenRV);
Vec sample(ConditionalExpression conditional_expression);

Var Sample(ConditionalExpression conditional_expression);

void sample(ConditionalExpression conditional_expression,Mat& samples);


RandomVar normal(real mean=0, real standard_dev=1, int d=1,
                 real minimum_standard_deviation=1e-6);

RandomVar normal(RandomVar mean, RandomVar log_variance,
                 real minimum_standard_deviation=1e-6);

RandomVar mixture(RVArray components, RandomVar log_weights);

RandomVar multinomial(RandomVar log_probabilities);

class StochasticRandomVariable: public RandomVariable
{
 public:
  StochasticRandomVariable(int length=1);
  StochasticRandomVariable(const RVArray& params,int length);
  StochasticRandomVariable(const RVArray& params,int length, int width);


  virtual bool isNonRandom() { return false; }

  virtual bool isDiscrete() { return false; }

  virtual void setKnownValues();

  
};

class FunctionalRandomVariable: public RandomVariable {
 public:
  FunctionalRandomVariable(int length);
  FunctionalRandomVariable(int length, int width);
  FunctionalRandomVariable(const Vec& the_value);
  FunctionalRandomVariable(const Mat& the_value);
  FunctionalRandomVariable(const Var& the_value);
  FunctionalRandomVariable(const RVArray& parents,int length);
  FunctionalRandomVariable(const RVArray& parents,int length, int width);


  virtual Var logP(const Var& obs, const RVInstanceArray& RHS,
      RVInstanceArray* parameters_to_learn);

  bool isNonRandom();

  virtual bool isDiscrete();


  
  virtual bool invertible(const Var& obs, 
                             RVInstanceArray& unobserved_parents,
                             Var** JacobianCorrection);

  virtual void setValueFromParentsValue() = 0;

  
};

class NonRandomVariable: public FunctionalRandomVariable
{
public:
  NonRandomVariable(int thelength);
  NonRandomVariable(int thelength, int thewidth);
  NonRandomVariable(const Var& v);

  virtual char* classname() { return "NonRandomVariable"; }

  void setValueFromParentsValue() { }
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection) 
    { return true; }
  void EMBprop(const Vec obs, real post) { }
};

class JointRandomVariable: public FunctionalRandomVariable
{
 public:
  JointRandomVariable(const RVArray& variables);

  virtual char* classname() { return "JointRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
};

class RandomElementOfRandomVariable: public FunctionalRandomVariable
{
public:
  RandomElementOfRandomVariable(const RandomVar& v, const RandomVar& index);

  virtual char* classname() { return "RandomElementOfRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);

  inline const RandomVar& v() { return parents[0]; }
  inline const RandomVar& index() { return parents[1]; }

};


class RVArrayRandomElementRandomVariable: public FunctionalRandomVariable
{
 public:
  RVArrayRandomElementRandomVariable(const RVArray& table, const RandomVar& index);

  virtual char* classname() { return "RVArrayRandomElementRandomVariable"; }

  void setValueFromParentsValue();
  virtual Var logP(const Var& obs, const RVInstanceArray& RHS,
      RVInstanceArray* parameters_to_learn=0);
  void EMBprop(const Vec obs, real post);

  inline const RandomVar& index() { return parents[parents.size()-1]; }

};

class NegRandomVariable: public FunctionalRandomVariable
{
 public:
  NegRandomVariable(RandomVariable* input);
  
  virtual char* classname() { return "NegRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
};

class ExpRandomVariable: public FunctionalRandomVariable
{
 public:
  ExpRandomVariable(RandomVar& input);
  
  virtual char* classname() { return "ExpRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
};


class LogRandomVariable: public FunctionalRandomVariable
{
 public:
  LogRandomVariable(RandomVar& input);
  
  virtual char* classname() { return "LogRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
};

class DiagonalNormalRandomVariable: public StochasticRandomVariable
{
 protected:
  real minimum_variance;
  real normfactor; 
  bool shared_variance; 

 public:
  DiagonalNormalRandomVariable(const RandomVar& mean, 
                               const RandomVar& log_variance,
                               real minimum_standard_deviation = 1e-10);

  virtual char* classname() { return "DiagonalNormalRandomVariable"; }

  Var logP(const Var& obs, const RVInstanceArray& RHS,
      RVInstanceArray* parameters_to_learn);
  void setValueFromParentsValue();
  void EMUpdate();
  void EMBprop(const Vec obs, real posterior);
  void EMEpochInitialize();

  inline const RandomVar& mean() { return parents[0]; }
  inline const RandomVar& log_variance() { return parents[1]; }
  inline bool& learn_the_mean() { return learn_the_parameters[0]; }
  inline bool& learn_the_variance() { return learn_the_parameters[1]; }

 protected:
  Vec mu_num; 
  Vec sigma_num; 
  real denom; 
};

class MixtureRandomVariable: public StochasticRandomVariable
{
 protected:
  RVArray components; 

 public:
  MixtureRandomVariable(const RVArray& components,
                        const RandomVar& log_weights);

  virtual char* classname() { return "MixtureRandomVariable"; }

  inline const RandomVar& log_weights() { return parents[0]; }
  inline bool& learn_the_weights() { return learn_the_parameters[0]; }

  virtual Var logP(const Var& obs, const RVInstanceArray& RHS,
      RVInstanceArray* parameters_to_learn);
  virtual Var ElogP(const Var& obs, RVInstanceArray& parameters_to_learn,
      const RVInstanceArray& RHS);

  virtual void setValueFromParentsValue();
  virtual void EMUpdate();
  virtual void EMBprop(const Vec obs, real posterior);
  virtual void EMEpochInitialize();
  virtual void EMTrainingInitialize(const RVArray& parameters_to_learn);
  virtual bool isDiscrete();
  virtual bool canStopEM();
  virtual void setKnownValues();
  virtual void unmarkAncestors();
  virtual void clearEMmarks();

 protected:
  Vec posteriors; 
  Vec sum_posteriors; 

  VarArray componentsLogP; 
  Var lw; 
  Var logp; 
};

class PlusRandomVariable: public FunctionalRandomVariable
{
 public:
  PlusRandomVariable(RandomVar input1, RandomVar input2);

  virtual char* classname() { return "PlusRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
  void EMTrainingInitialize(const RVArray& parameters_to_learn);
  void EMEpochInitialize();
  void EMUpdate();

  const RandomVar& X0() { return parents[0]; }
  const RandomVar& X1() { return parents[1]; }

  bool learn_X0() { return learn_the_parameters[0]; }
  bool learn_X1() { return learn_the_parameters[1]; }
  bool learn_something;
  RandomVar parent_to_learn; 
  RandomVar other_parent; 
  Vec numerator;
  Vec difference;
  real denom;
};

class MinusRandomVariable: public FunctionalRandomVariable
{
 public:
  MinusRandomVariable(RandomVar input1, RandomVar input2);

  virtual char* classname() { return "MinusRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
  void EMTrainingInitialize(const RVArray& parameters_to_learn);
  void EMEpochInitialize();
  void EMUpdate();

  const RandomVar& X0() { return parents[0]; }
  const RandomVar& X1() { return parents[1]; }

  bool learn_X0() { return learn_the_parameters[0]; }
  bool learn_X1() { return learn_the_parameters[1]; }
  bool learn_something;
  RandomVar parent_to_learn; 
  RandomVar other_parent; 
  Vec numerator;
  Vec difference;
  real denom;
};


class ElementWiseDivisionRandomVariable: public FunctionalRandomVariable
{
 public:
  ElementWiseDivisionRandomVariable(RandomVar input1, RandomVar input2);

  virtual char* classname() { return "ElementWiseDivisionRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
  void EMTrainingInitialize(const RVArray& parameters_to_learn);
  void EMEpochInitialize();
  void EMUpdate();

  const RandomVar& X0() { return parents[0]; }
  const RandomVar& X1() { return parents[1]; }

};


class ProductRandomVariable: public FunctionalRandomVariable
{
 public:
  int m,n,l; 

  ProductRandomVariable(MatRandomVar input1, MatRandomVar input2);

  virtual char* classname() { return "ProductRandomVariable"; }

  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real post);
  void EMTrainingInitialize(const RVArray& parameters_to_learn);
  void EMEpochInitialize();
  void EMUpdate();

  const RandomVar& X0() { return parents[0]; } 
  const RandomVar& X1() { return parents[1]; } 
  bool scalars; 

  bool learn_X0() { return learn_the_parameters[0]; }
  bool learn_X1() { return learn_the_parameters[1]; }
  bool learn_something;
  Mat X0numerator; 
  Mat X1numerator; 
  Mat denom; 
  Mat tmp1; 
  Mat tmp2; 
  Mat tmp3; 
  Vec vtmp3; 
  Vec tmp4;
};

class SubVecRandomVariable: public FunctionalRandomVariable
{
 protected:
  int start;
 public:
  SubVecRandomVariable(const RandomVar& parent,int start, int length);
  virtual char* classname() { return "SubvecRandomVariable"; }
  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real posterior);
};

class MultinomialRandomVariable: public StochasticRandomVariable
{
 public:
  MultinomialRandomVariable(const RandomVar& log_probabilities);

  inline const RandomVar& log_probabilities() { return parents[0]; }
  inline bool learn_the_probabilities() { return learn_the_parameters[0]; }

  virtual char* classname() { return "MultinomialRandomVariable"; }

  Var logP(const Var& obs, const RVInstanceArray& RHS, 
      RVInstanceArray* parameters_to_learn);
  void setValueFromParentsValue(); 
  void EMUpdate();
  void EMBprop(const Vec obs, real posterior);
  void EMEpochInitialize();
  bool isDiscrete();

 protected:
  Vec sum_posteriors; 
};


class ExtendedRandomVariable: public FunctionalRandomVariable
{
 protected:
  int n_extend;
  real fill_value;
 public:
  ExtendedRandomVariable(const RandomVar& parent, real fill_value=1.0,int n_extend=1);
  virtual char* classname() { return "ExtendedRandomVariable"; }
  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real posterior);
};

class ConcatColumnsRandomVariable: public FunctionalRandomVariable
{
 public:
  ConcatColumnsRandomVariable(const RVArray& vars);
  virtual char* classname() { return "ConcatColumnsRandomVariable"; }
  void setValueFromParentsValue();
  bool invertible(const Var& obs, RVInstanceArray& unobserved_parents,
                     Var** JacobianCorrection);
  void EMBprop(const Vec obs, real posterior);
};



class RandomVarVMatrix: public VMatrix
{
 protected:
  RandomVar rv;
  Var instance;
  VarArray prop_path;

 public:
  RandomVarVMatrix(ConditionalExpression conditional_expression);
  virtual int nVars() { return instance->length(); }
  virtual Vec sample()
    {
      prop_path.fprop();
      return instance->value;
    }
};


} // end of namespace PLearn

#endif

---