/src/C/hierarchical.c

  1/*
  2	Copyright (c) 2011, Yahoo! Inc.  All rights reserved.
  3	Copyrights licensed under the New BSD License. See the accompanying LICENSE file for terms.
  4
  5    Author: Bee-Chung Chen
  6*/
  7
  8
  9#include <Rinternals.h>
 10#include "util.h"
 11#include <stdio.h>
 12
 13void hierarchical_smoothing_1D_2Levels(
 14	// Output
 15	double *a_mean, double *a_var, // nItems x 1
 16	double *b_mean, double *b_var, double *b_cov, // nItems x nCategories
 17	// Input
 18	const int *itemIndex, const int *categoryIndex, const double *obs, const double *var_obs, // nObs x 1
 19	const double *q, const double *var_b, // nCategories x 1
 20	const double *var_a, // 1x1
 21	const int *numObs, const int *numItems, const int *numCategories,
 22	// Control
 23	const int *verbose, const int *debug
 24){
 25	int nObs = (*numObs), nItems = (*numItems), nCategories = (*numCategories);
 26
 27	// Compute sufficient statistics
 28	// 		b_mean[i,k] = sum_j { o_{ijk}/var_obs_{ijk} } = C_ik
 29	//       b_var[i,k] = sum_j {       1/var_obs_{ijk} } = F_ik
 30	for(int n=0; n < nItems*nCategories; n++){ b_mean[n] = 0; b_var[n] = 0; }
 31	for(int j=0; j<nObs; j++){
 32		double o = obs[j];
 33		int i = itemIndex[j] - 1;
 34		int k = categoryIndex[j] - 1;
 35		CHK_C_INDEX(i,nItems); CHK_C_INDEX(k,nCategories);
 36		// double var_obs_plus_var_b = var_obs[j] + var_b[k];
 37		b_mean[C_MAT(i,k,nItems)] += o / var_obs[j];
 38		b_var[C_MAT(i,k,nItems)]  += 1 / var_obs[j];
 39	}
 40
 41	// Compute E[a[i] | obs] and Var[a[i] | obs]
 42	for(int i=0; i<nItems; i++){
 43		double sum_for_mean = 0; // sum_k E[a_i | obs_ik]/Var[a_i | obs_ik]
 44		double sum_for_var  = 0; // sum_k (1/Var[a_i | obs_ik] - 1/var_a)
 45		for(int k=0; k<nCategories; k++){
 46			double F = b_var[C_MAT(i,k,nItems)];
 47			if(F == 0) continue;
 48			double C = b_mean[C_MAT(i,k,nItems)];
 49			double A = (1/var_b[k]) + F;
 50			// E[a_i | obs_ik]/Var[a_i | obs_ik] = (C * q[k]) / (A * var_b[k])
 51			// (1/Var[a_i | obs_ik] - 1/var_a)   = (q[k]^2 / var_b[k]) * (1 - 1/(A * var_b[k]))
 52			sum_for_mean += (C * q[k]) / (A * var_b[k]);
 53			sum_for_var  += (q[k]*q[k] / var_b[k]) * (1 - 1/(A * var_b[k]));
 54		}
 55		if((*verbose) >= 100) printf("   i=%d:  sum.for.a.mean=%f, sum.for.a.var=%f\n", i+1, sum_for_mean, sum_for_var);
 56
 57		// Compute E[a[i] | obs] and Var[a[i] | obs]
 58		a_var[i]  = 1/( (1/var_a[0]) + sum_for_var );
 59		a_mean[i] = a_var[i] * sum_for_mean;
 60	}
 61
 62	// Compute E[b[i,k] | obs], Var[b[i,k] | obs], Cov[b[i,k], a[i] | obs]
 63	for(int i=0; i<nItems; i++){ for(int k=0; k<nCategories; k++){
 64		int ik = C_MAT(i,k,nItems);
 65		double A = 1 / ( (1/var_b[k]) + b_var[ik] );
 66		double D = (A * q[k]) / var_b[k];
 67		b_mean[ik] = D * a_mean[i] + A * b_mean[ik];
 68		b_var[ik]  = A + D*D*a_var[i];
 69		b_cov[ik]  = D * a_var[i];
 70	}}
 71}
 72SEXP hierarchical_smoothing_1D_2Levels_Call(
 73  // Output
 74  SEXP a_mean, SEXP a_var, // nItems x 1
 75  SEXP b_mean, SEXP b_var, SEXP b_cov, // nItems x nCategories
 76  // Input
 77  SEXP itemIndex, SEXP categoryIndex, SEXP obs, SEXP var_obs, // nObs x 1
 78  SEXP q, SEXP var_b, // nCategories x 1
 79  SEXP var_a, // 1x1
 80  SEXP numObs, SEXP numItems, SEXP numCategories,
 81  // Control
 82  SEXP verbose, SEXP debug
 83){
 84  hierarchical_smoothing_1D_2Levels(
 85    // Output
 86    MY_REAL(a_mean), MY_REAL(a_var), // nItems x 1
 87    MY_REAL(b_mean), MY_REAL(b_var), MY_REAL(b_cov), // nItems x nCategories
 88    // Input
 89    MY_INTEGER(itemIndex), MY_INTEGER(categoryIndex), MY_REAL(obs), MY_REAL(var_obs), // nObs x 1
 90    MY_REAL(q), MY_REAL(var_b), // nCategories x 1
 91    MY_REAL(var_a), // 1x1
 92    MY_INTEGER(numObs), MY_INTEGER(numItems), MY_INTEGER(numCategories),
 93    // Control
 94    MY_INTEGER(verbose), MY_INTEGER(debug)
 95  );
 96  return R_NilValue;
 97}
 98
 99
100void hierarchical_smoothing_1D_2Levels_incorrect(
101	// Output
102	double *a_mean, double *a_var, // nItems x 1
103	double *b_mean, double *b_var, double *b_cov, // nItems x nCategories
104	// Input
105	const int *itemIndex, const int *categoryIndex, const double *obs, const double *var_obs, // nObs x 1
106	const double *q, const double *var_b, // nCategories x 1
107	const double *var_a, // 1x1
108	const int *numObs, const int *numItems, const int *numCategories,
109	// Control
110	const int *verbose, const int *debug
111){
112	int nObs = (*numObs), nItems = (*numItems), nCategories = (*numCategories);
113
114	// Compute sufficient statistics
115	// 		  a_mean[i] = sum{ (q_k * o_{ijk})/(var_obs_{ijk} + var_b_{k}) }
116	//         a_var[i] = sum{ (q_k * q_k)/(var_obs_{ijk} + var_b_{k}) }
117	// 		b_mean[i,k] = sum{ o_{ijk}/var_obs_{ijk} }
118	//       b_var[i,k] = sum{       1/var_obs_{ijk} }
119	for(int i=0; i<nItems; i++){ a_mean[i] = 0; a_var[i] = 0; }
120	for(int n=0; n < nItems*nCategories; n++){ b_mean[n] = 0; b_var[n] = 0; }
121	for(int j=0; j<nObs; j++){
122		double o = obs[j];
123		int i = itemIndex[j] - 1;
124		int k = categoryIndex[j] - 1;
125		CHK_C_INDEX(i,nItems); CHK_C_INDEX(k,nCategories);
126		double var_obs_plus_var_b = var_obs[j] + var_b[k];
127		a_mean[i] += (q[k]*o) / var_obs_plus_var_b;
128		a_var[i]  += (q[k]*q[k]) / var_obs_plus_var_b;
129		b_mean[C_MAT(i,k,nItems)] += o / var_obs[j];
130		b_var[C_MAT(i,k,nItems)]  += 1 / var_obs[j];
131	}
132
133	// Compute E[a[i] | obs] and Var[a[i] | obs]
134	for(int i=0; i<nItems; i++){
135		if((*verbose) >= 100) printf("   i=%d:  sum.for.a.mean=%f, sum.for.a.var=%f\n", i+1, a_mean[i], a_var[i]);
136		a_var[i]  = 1/( (1/var_a[0]) + a_var[i] );
137		a_mean[i] = a_var[i] * a_mean[i];
138	}
139
140	// Compute E[b[i,k] | obs], Var[b[i,k] | obs], Cov[b[i,k], a[i] | obs]
141	for(int i=0; i<nItems; i++){ for(int k=0; k<nCategories; k++){
142		int ik = C_MAT(i,k,nItems);
143		double A = 1 / ( (1/var_b[k]) + b_var[ik] );
144		double D = (A * q[k]) / var_b[k];
145		b_mean[ik] = D * a_mean[i] + A * b_mean[ik];
146		b_var[ik]  = A + D*D*a_var[i];
147		b_cov[ik]  = D * a_var[i];
148	}}
149}