normal_conjugate.cweb

@q $Id$ @>
@q Copyright 2007, Ondra Kamenik @>

@ Start of {\tt normal\_conjugate.cpp} file.

@c

#include "normal_conjugate.h"
#include "kord_exception.h"

@<|NormalConj| diffuse prior constructor@>;
@<|NormalConj| data update constructor@>;
@<|NormalConj| copy constructor@>;
@<|NormalConj::update| one observation code@>;
@<|NormalConj::update| multiple observations code@>;
@<|NormalConj::update| with |NormalConj| code@>;
@<|NormalConj::getVariance| code@>;

@ 
@<|NormalConj| diffuse prior constructor@>=
NormalConj::NormalConj(int d)
	: mu(d), kappa(0), nu(-1), lambda(d,d)
{
	mu.zeros();
	lambda.zeros();
}

@ 
@<|NormalConj| data update constructor@>=
NormalConj::NormalConj(const ConstTwoDMatrix& ydata)
	: mu(ydata.numRows()), kappa(ydata.numCols()), nu(ydata.numCols()-1),
	  lambda(ydata.numRows(), ydata.numRows())
{
	mu.zeros();
	for (int i = 0; i < ydata.numCols(); i++)
		mu.add(1.0/ydata.numCols(), ConstVector(ydata, i));

	lambda.zeros();
	for (int i = 0; i < ydata.numCols(); i++) {
		Vector diff(ConstVector(ydata, i));
		diff.add(-1, mu);
		lambda.addOuter(diff);
	}
}

@ 
@<|NormalConj| copy constructor@>=
NormalConj::NormalConj(const NormalConj& nc)
	: mu(nc.mu), kappa(nc.kappa), nu(nc.nu), lambda(nc.lambda)
{
}

@ The method performs the following:
$$\eqalign{
  \mu_1 = &\; {\kappa_0\over \kappa_0+1}\mu_0 + {1\over \kappa_0+1}y\cr
  \kappa_1 = &\; \kappa_0 + 1\cr
  \nu_1 = &\; \nu_0 + 1\cr
  \Lambda_1 = &\; \Lambda_0 + {\kappa_0\over\kappa_0+1}(y-\mu_0)(y-\mu_0)^T,
}$$

@<|NormalConj::update| one observation code@>=
void NormalConj::update(const ConstVector& y)
{
	KORD_RAISE_IF(y.length() != mu.length(),
				  "Wrong length of a vector in NormalConj::update");

	mu.mult(kappa/(1.0+kappa));
	mu.add(1.0/(1.0+kappa), y);

	Vector diff(y);
	diff.add(-1, mu);
	lambda.addOuter(diff, kappa/(1.0+kappa));

	kappa++;
	nu++;
}

@ The method evaluates the formula in the header file.

@<|NormalConj::update| multiple observations code@>=
void NormalConj::update(const ConstTwoDMatrix& ydata)
{
	NormalConj nc(ydata);
	update(nc);
}


@ 
@<|NormalConj::update| with |NormalConj| code@>=
void NormalConj::update(const NormalConj& nc)
{
	double wold = ((double)kappa)/(kappa+nc.kappa);
	double wnew = 1-wold;

	mu.mult(wold);
	mu.add(wnew, nc.mu);

	Vector diff(nc.mu);
	diff.add(-1, mu);
	lambda.add(1.0, nc.lambda);
	lambda.addOuter(diff);

	kappa = kappa + nc.kappa;
	nu = nu + nc.kappa;
}


@ This returns ${1\over \nu-d-1}\Lambda$, which is the mean of the
variance in the posterior distribution. If the number of degrees of
freedom is less than $d$, then NaNs are returned.

@<|NormalConj::getVariance| code@>=
void NormalConj::getVariance(TwoDMatrix& v) const
{
	if (nu > getDim()+1) {
		v = (const TwoDMatrix&)lambda;
		v.mult(1.0/(nu-getDim()-1));
	} else
		v.nans();
}


@ End of {\tt normal\_conjugate.cpp} file.