utilities.c
/************************************************************/
/* Utility functions */
/* Author: Wayne Zhang */
/* actuary_zhang@hotmail.com */
/************************************************************/
/**
* @file utilities.c
* @brief Utility functions for simple arithmetic calculation,
* matrix computation and numerical derivatives.
* @author Wayne Zhang
*/
#include "common.h"
#include "utilities.h"
/** constant used in numerical differentiation */
#define DIFF_EPS 0.001
/************************************************************/
/* Simple arithmetic utility */
/************************************************************/
/**
* Compute sample variance
*
* @param x samples in long vector
* @param n number of samples
*
* @return sample variance
*/
double var(double *x, int n){
double ans = 0.0, m = mean(x, n) ;
for (int i = 0; i < n; i++)
ans += (x[i] - m) * (x[i] - m);
ans /= n - 1.0 ;
return ans;
}
/**
* Compute sample covariance matrix
*
* @param n number of samples
* @param p number of variables (columns), p >2
* @param x samples in long vector
* @param ans vector to store computed covariance matrix
*
*/
void cov(int n, int p, double *x, double *ans){
double *one = Calloc(n * n, double),
*x2 = Calloc(n * p, double),
*x3 = Calloc(n * p, double);
double alpha = -1.0 / n, beta = 1.0, beta2 = 0.0;
// subtract mean
for (int i = 0; i < n * n; i++) one[i] = 1.0 ;
Memcpy(x2, x, n * p) ;
Memcpy(x3, x, n * p);
F77_CALL(dgemm)("N", "N", &n, &p, &n, &alpha, one,
&n, x2, &n, &beta, x3, &n);
Memcpy(x2, x3, n * p) ;
AZERO(ans, p * p) ;
// compute covariance
F77_CALL(dgemm)("T", "N", &p, &p, &n, &beta, x2,
&n, x3, &n, &beta2, ans, &p);
for (int i = 0; i < p * p; i++) ans[i] /= (n - 1) * 1.0 ;
Free(one) ;
Free(x2) ;
Free(x3) ;
}
/************************************************************/
/* Matrix computations */
/************************************************************/
/**
* Multiply a matrix and a vector
*
* @param trans transpose of matrix?
* @param m row count of matrix
* @param n column count of matrix
* @param A input matrix
* @param x input vector
* @param out output vector
*
*/
void mult_mv(char *trans, int m, int n, double *A,
double *x, double *out){
double one = 1.0, zero = 0.0 ;
int incx = 1;
F77_CALL(dgemv)(trans, &m, &n, &one, A, &m, x, &incx,
&zero, out, &incx) ;
}
/**
* compute t(x) * x
*
* @param m row dimension of the matrix
* @param n column dimension of the matrix
* @param x the input matrix
* @param out output results
*
*/
void mult_xtx(int m, int n, double *x, double *out){
double alpha = 1.0, beta = 0.0, *x2 = Calloc(m * n, double);
Memcpy(x2, x, m * n) ;
F77_CALL(dgemm)("T", "N", &n, &n, &m, &alpha, x2, &m,
x, &m, &beta, out, &n) ;
Free(x2) ;
}
/**
* compute the lower cholesky factor
*
* @param d dimension of the matrix
* @param v input matrix
* @param iv output cholesky factor
*
*/
void chol(int d, double *v, double *iv){
int info = 0;
// cholesky factor of v
Memcpy(iv, v, d * d) ;
F77_CALL(dpotrf)("L", &d, iv, &d, &info) ;
if (info) error(_("Error %d in Cholesky decomposition."), info) ;
}
/**
* invert a positive symmetric matrix
*
* @param d dimension of the matrix
* @param v input matrix
* @param iv output inverse of the matrix
*
*/
void solve_po(int d, double *v, double *iv){
int info = 0;
// cholesky factor of v
chol(d, v, iv) ;
// compute inverse
F77_CALL(dpotri)("L", &d, iv, &d, &info) ;
if (info) error(_("Error %d in inverting matrix."), info) ;
// fill upper triangle
for (int i = 0; i < d - 1; i++){
for (int j = i + 1; j < d; j++)
iv[j * d + i] = iv[i * d + j] ;
}
}
/************************************************************/
/* Compute numerical derivatives */
/************************************************************/
/**
* Compute numerical gradient
*
* @param n length of parmaters
* @param x values at which to evaluate the gradient
* @param myfunc user specified function
* @param data struct used in myfunc
* @param ans vector to store the gradient
*
*/
void grad(int n, double *x, double (*myfunc)(double *x, void *data),
void *data, double *ans){
double y1, y2 ;
for (int i = 0; i < n; i++){
x[i] += DIFF_EPS ;
y1 = myfunc(x, data) ;
x[i] -= 2 * DIFF_EPS ;
y2 = myfunc(x, data) ;
ans[i] = (y1 - y2) / DIFF_EPS * 0.5 ;
x[i] += DIFF_EPS ;
}
}
/**
* Compute numerical hessian matrix
*
* @param n length of parmaters
* @param x values at which to evaluate the hessian
* @param myfunc user specified function
* @param data struct used in myfunc
* @param ans n*n vector to store the hessian matrix
*
*/
void hess(int n, double *x, double (*myfunc)(double *x, void *data),
void *data, double *ans){
double *y1 = Calloc(n, double),
*y2 = Calloc(n, double) ;
for (int i = 0; i < n; i++){
x[i] += DIFF_EPS ;
grad(n, x, myfunc, data, y1) ;
x[i] -= 2 * DIFF_EPS ;
grad(n, x, myfunc, data, y2) ;
for (int j = 0; j < n; j++)
ans[j + i * n] = (y1[j] - y2[j]) / DIFF_EPS * 0.5 ;
x[i] += DIFF_EPS ;
}
Free(y1) ; Free(y2) ;
}
/************************************************************/
/* Simulate a Wishart variable */
/************************************************************/
/**
* Simulate the Cholesky factor of a standardized Wishart variate with
* dimension p and nu degrees of freedom.
*
* @param nu degrees of freedom
* @param p dimension of the Wishart distribution
* @param upper if 0 the result is lower triangular, otherwise upper
triangular
* @param ans array of size p * p to hold the result
*
* @return ans
*/
static double *std_rWishart_factor(double nu, int p, int upper, double ans[])
{
int pp1 = p + 1;
if (nu < (double) p || p <= 0)
error(_("inconsistent degrees of freedom and dimension"));
AZERO(ans, p * p);
for (int j = 0; j < p; j++) { /* jth column */
ans[j * pp1] = sqrt(rchisq(nu - (double) j));
for (int i = 0; i < j; i++) {
int uind = i + j * p, /* upper triangle index */
lind = j + i * p; /* lower triangle index */
ans[(upper ? uind : lind)] = norm_rand();
ans[(upper ? lind : uind)] = 0;
}
}
return ans;
}
/**
* Simulate a sample of random matrix from a Wishart distribution
*
* @param d row (=column) dimension of the matrix
* @param nu Degrees of freedom
* @param scal Positive-definite scale matrix
* @param out simulated matrix (d*d)
*
*/
void rwishart(int d, double nu, double *scal, double *out)
{
int info, psqr;
double *scCp, *tmp, one = 1, zero = 0;
psqr = d*d;
tmp = Calloc(psqr, double);
scCp = Calloc(psqr, double);
Memcpy(scCp, scal, psqr);
AZERO(tmp, psqr);
F77_CALL(dpotrf)("U", &d, scCp, &d, &info);
if (info)
error(_("scale matrix is not positive-definite"));
GetRNGstate();
std_rWishart_factor(nu, d, 1, tmp);
F77_CALL(dtrmm)("R", "U", "N", "N", &d, &d,
&one, scCp, &d, tmp, &d);
F77_CALL(dsyrk)("U", "T", &d, &d, &one, tmp, &d,
&zero, out, &d);
for (int i = 1; i < d; i++){
for (int k = 0; k < i; k++)
out[i + k * d] = out[k + i * d];
}
PutRNGstate();
Free(tmp) ;
Free(scCp) ;
}
