Version 1.1 of the MATLAB toolbox can be downloaded from here. Please see “examples_bayesreg.m” for examples on how to use the toolbox, or type “help bayesreg” within MATLAB.
An R version of this toolbox is now available on CRAN. To install the R package, type “install.packages(“bayesreg”)” within R.
To cite this toolbox:
Updates:
(23/01): Updated the bayesreg code to version 1.2. This version implements Zellner’s gprior for linear and logistic regression. The gprior only works with full rank matrices. The examples in “examples_bayesreg.m” have been updated to include a gprior example.
(17/01): Updated the bayesreg code to version 1.1. The display code is now in a separate “summary()” function so that it can be used many times on the same data and sampling results. We have updated “examples_bayesreg.m” to include examples of the new “summary()” function.
]]>I hope this document will be useful to researchers interested in learning more about the the MML principle.
Update (24/07/2014): I have written MATLAB code to compute MML estimates and the message length. I have also done minor updates to the document.
]]>TEAM  ELO  WINS  LOSSES  DRAWS  TOTAL  

1  Hawthorn  1475  869  951  10  1830 
2  Geelong  1469  1224  1042  21  2287 
3  Collingwood  1468  1460  910  26  2396 
4  Sydney  1428  1061  1199  23  2283 
5  West Coast  1360  339  269  5  613 
6  Brisbane Bears  1319  72  148  2  222 
7  Adelaide  1316  269  241  1  511 
8  St Kilda  1311  877  1338  25  2240 
9  North Melbourne  1293  803  1005  17  1825 
10  Fremantle  1292  168  236  0  404 
11  Carlton  1253  1385  939  33  2357 
12  Essendon  1194  1315  973  34  2322 
13  Richmond  1170  1058  1035  22  2115 
14  Brisbane Lions  1147  192  175  6  373 
15  Western Bulldogs  1099  803  976  22  1801 
16  Melbourne  1032  1034  1209  21  2264 
17  Port Adelaide  1032  183  181  5  369 
18  GW Sydney  1016  2  20  0  22 
19  Gold Coast  992  6  38  0  44 
20  Fitzroy  871  869  1034  25  1928 
21  University  781  27  97  2  126 
Note:
It would be of interest to use this data for prediction of the AFL 2013 season, perhaps using a methodology similar to the ELO++ work.
]]>Suppose you want to do a QR decomposition of a matrix X. There are three classes that can do this for you, one of which is ColPivHouseholderQR, and you may write something like
// perform a QR decomposition of the [n x p] X matrix; where n > p
Eigen::ColPivHouseholderQR< Eigen::MatrixXd > qr = X.colPivHouseholderQr();
Seems fairly straightforward, right? But, how do we get Q and R? If we look at the online documentation for ColPivHouseholderQR, we find a public method called matrixQR() which seems to return a reference to a matrix related to Q and R in some fashion. The documentation doesn’t state anything about how to get Q and R from the returned matrix. It turns out, that the matrix returned by matrixQR() contains both Q and R, with R occupying the ‘top’ block and Q the ‘bottom’ block. Thus, to get, say, R we need to use something like
// perform a QR decomposition of the [n x p] X matrix; where n > p
Eigen::ColPivHouseholderQR< Eigen::MatrixXd > qr = X.colPivHouseholderQr();
Eigen::Matrix R;
R = qr.matrixQR().topRightCorner(p,p).triangularView(); // R matrix
Here, topRightCorner() extracts the matrix of interest and triangularView
It didn’t take too long to figure out how QR works in Eigen. After discovering the solution I found a post on the Eigen forum that describes a similar solution to the aforementioned. I urge anyone who uses Eigen to contribute code and documentation to the project, if possible. It really is a fantastic opensource library and equals and surpasses many other commercial linear algebra packages in both performance and quality.
]]>where , is an i.i.d. zeromean multivariate Gaussian random variable with a symmetric, positive definite variancecovariance matrix ; that is, for all . It is thus assumed that each data vector follows a variate Gaussian distribution with moments
The aim of single factor analysis is to estimate the unknown factor scores , the mean , the variancecovariance matrix and the factor loadings .
I’ve written this note to show how one may derive the Fisher information for the single factor model. This turns out to be a fairly involved task requiring a lot of linear algebra and matrix differential calculus identities. I’ve uploaded a document describing all the necessary mathematical steps to the Publications section. Anyone interested in matrix differential calculus and linear algebra should find it useful.
]]>z=1.96;
2*(1normcdf(z))
However, the above only works for relatively small z, in the range z < 10. In order to get pvalues for z as large as 30, one should instead run
z=1.96;
2*normcdf(abs(z))
For those using R, the above can be accomplished by running:
z < 1.96
2*pnorm(z,lower.tail=F)
]]>#include
#include
#include
boost::mt19937 rng; // random number generator
boost::random::uniform_01 dist;
/*
* Generate n samples from a knomial distribution.
* Input:
* tab  vector of length k to store the counts;
*  tab[j] \in [0, n] for all j=1,...,k
*  sum_{j=0}^k tab[j] = n
* p  the probability vector of length k
*  does not assume p(0)+p(1)+...+p(k)=1;
* k  knomial (k>1)
* n  sample size
*/
void sample_knomial(unsigned int *tab, double *p, unsigned int k, unsigned int n)
{
assert(k>1);
// set all kcells to zero
memset(tab,0,sizeof(unsigned int)*k);
// normalise the probabilities
double norm_factor=0.0;
for(unsigned int j=0; j 0);
norm_factor += p[j];
}
// precompute cumulative sums [p(0),p(0)+p(1),...,p(0)+p(1)+...+p(k1)]
// (could sort p prior to this)
double *cumsum = new double[k];
cumsum[0]=p[0]/norm_factor;
for(unsigned int j=1; jcumsum[t]; t++);
// Increment count for cell j
tab[t]++;
}
// free memory
delete[] cumsum;
}
]]>Below is a simple C function to compute the HL estimator for a given array of observations.
// Comparison function needed for quicksort
int compare (const void * a, const void * b)
{
const double *da = (const double *) a;
const double *db = (const double *) b;
return (*da > *db)  (*da < *db);
}
// Function to compute the HodgesLehmann estimator
double hodgeslehmann(double *x, int n)
{
// store and compute Walsh averages
unsigned int nw=n*(n+1)/2;
double *walsh_avgs=(double*) malloc(sizeof(double)*nw);
unsigned int k=0;
for(unsigned int i=0; i
I have also written a MATLAB driver function that allows the aforementioned C code to be executed within MATLAB, see the Software section to download this implementation. An efficient (pure) MATLAB implementation that avoids the use of nested for loops is a little bit tricky. Here is a simple (memmory intensive!) implementation that relies on the ndgrid(.) function. This code works reasonably well as long as the sample size is not too big.
function m=hodglehm_slow(x)
%% Length of vector
x=x(:); % make sure x is a column vector
n=length(x);
%% compute cartesian product
[a,b]=ndgrid(x,x);
%% Sum of all pairs
A=triu((a+b)./2);
%% Extract upper diagonal
Ix=triu(true(n));
A=A(Ix);
%% Compute median of averages
m=median(A(:));
%% done
return;
We can get the standard error of the HL estimator in MATLAB by using the bootstrap:
% Generate 50 observations from a Studentt distribution with 2 dof
x=trnd(2,50,1);
% Compute standard errors for the mean and HL estimators
std(bootstrp(1000, @(a) mean(a), x))
std(bootstrp(1000, @(a) hodglehm_slow(a), x))
]]>
http://www.emakalic.org/blog/?feed=rss2&p=78
0

Hyperspherical Coordinates
http://www.emakalic.org/blog/?p=75
Tue, 14 Jun 2011 05:00:58 +0000
http://www.emakalic.org/blog/?p=75
Suppose we have a vector residing in an dimensional Euclidean space. How do we express this vector using (hyper)spherical coordinates in MATLAB? MATLAB implements a couple of functions, namely cart2pol() and cart2sph() , that handle the conversion in two and three dimensions, respectively. However, there is no general method for vectors of arbitrary dimension.
Fortunately, Wikipedia lists the formulas required to convert Euclidean coordinates to hyperspherical coordinates and these are relatively straightforward to implement in MATLAB. Below is a MATLAB function that takes an dimensional Euclidean vector as input and returns the corresponding radial coordinate as well as the angular coordinates ; note that, , while all the other angles are in the set .
%% The input vector x is a ndim. column vector.
function [r,theta]=cart2nsph(x)
%% setup
n = length(x);
theta = zeros(n1,1);
x2 = x.^2;
cumx = cumsum(x2(end:1:2));
%% Compute the radial coordinate
r = sqrt(sum(x2));
%% Compute the angles
theta(1:end1,1) = x(1:end2,1) ./ sqrt( cumx(end:1:2) ) ;
theta(end,1) = (sqrt(sum(x2(end:1:end1,1))) + x(end1,1)) ./ x(end,1);
theta = acot(theta);
theta(theta<0) = theta(theta<0)+pi; % Convert acot principal value to [0,pi]
theta(end,1) = 2*theta(end,1); % Last angle is between [0,2pi]
%% done
return;
Line 11 computes the radial coordinate which is simply the square root of the sum of the squares of the original, Euclidean coordinates. Lines 1418 are required to obtain the angular coordinates. Some care is necessary in order to get all the angles in the correct range, since the principal value of the function acot() in MATLAB is in and not in the (usual) range of . We correct for this in Line 17 by adding to all the "negative" angles.
The function to compute dimensional Euclidean coordinates from hyperspherical coordinates is somewhat simpler. One possible implementation is given below:
%% Radial coordinate, r>0, must be positive.
%% The angles theta_1,...,theta_{n2} must be in [0,pi], while
%% theta_{n1} should be in [0,2pi].
function x = nsph2cart(r,theta)
c=[cos(theta); 1];
s=[1 ;cumprod(sin(theta))];
x=r*(c.*s);
%% done
return;
Except in some special cases the transform is unique. Both functions were tested on MATLAB version 7.12.0 (R2011a) and dont have any error checking. I hope that some of you can find these functions useful in your work.
]]>