#### /src/C/hierarchical.c

http://github.com/beechung/Latent-Factor-Models
C | 149 lines | 103 code | 12 blank | 34 comment | 16 complexity | 87c878f1948eac516c44e4fc5cfa80f0 MD5 | raw file
```  1/*
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}
```