swh:1:snp:ffdd0a7d2c8ea15ad41d45b3b178f668bd942287
Tip revision: ce0bcae35573fa8b1dfb9f24425358852d41882a authored by Derek Young on 10 March 2017, 07:50:20 UTC
version 1.1.0
version 1.1.0
Tip revision: ce0bcae
KDElocscale.c
#include <R.h>
#include <Rmath.h>
/* Implement kernel density-estimation step for location-scale
mixture model (17) in Benaglia et al (2008) where each
component and block has the same shape but may differ from
others by a location and scale, implementing equation (18) */
void KDElocscale(
int *nn, /* Sample size */
int *mm, /* Number of components */
int *rr, /* Number of coordinates */
int *blockid, /* r-vector of block numbers */
double *mu, /* m by max(blockid) matrix of current mean estimates */
double *sigma, /* m by max(blockid) matrix of current stdev estimates */
double *x, /* n by r data matrix */
double *hh, /* scalar bandwidth */
double *z, /* n by m vector of normalized posteriors (or indicators in
stochastic case), normalized by "column" */
double *f /* n by m matrix of KDE products */
) {
int n=*nn, m=*mm, r=*rr;
int i, j, k, ii, kk, ell, ell2;
double sum1, sum2, tmp, h=*hh, u;
double const1 = -0.5 / (h * h);
double const2;
for(j=0; j<m; j++) {
/* In const2 expression, .3989... equals 1/(sqrt(2*pi)) */
const2 = 0.39894228040143267794/(h*sigma[j]*(double)r);
for (i=0; i<n; i++) {
f[i + j*n] = 1.0;
for (k=0; k<r; k++) {
ell = blockid[k]-1;
u = (x[i + k*n] - mu[j + m*ell])/sigma[j + m*ell];
sum1 = 0.0;
for (ii=0; ii<n; ii++) {
sum2 = 0.0;
for (kk=0; kk<r; kk++) {
ell2 = blockid[kk]-1;
tmp = (u - x[ii+kk*n] + mu[j + m*ell2])/sigma[j + m*ell2];
sum2 += exp(tmp * tmp * const1);
}
sum1 += z[ii + j*n] * sum2;
}
f[i + j*n] *= sum1 * const2;
}
}
}
}
//THE FOLLOWING IS NOT FINISHED!
//
///* Implement kernel density-estimation step for location-scale
// mixture model with a different shape for each component,
// equation (11) in Benaglia et al (2008) */
//void KDElocscale(
// int *nn, /* Sample size */
// int *mm, /* Number of components */
// int *rr, /* size of current block */
// double *mu, /* m-vector of current mean estimates */
// double *sigma, /* m-vector of current stdev estimates */
// double *x, /* data: vector of length n */
// double *hh, /* scalar bandwidth */
// double *z, /* nn*mm vector of normalized posteriors (or indicators in
// stochastic case), normalized by "column" */
// double *f /* nxm matrix of KDE products */
//) {
// int n=*nn, i, ii;
// int j, ell;//k, kk;
// double sum1, sum2, tmp, h=*hh, xik, u;
// double const1 = -0.5 / (h * h);
// double const2;
//
// for(j=0; j<m; j++) {
// /* In const2 expression, .3989... equals 1/(sqrt(2*pi)) */
// const2 = 0.39894228040143267794/(h*sigma[j]*(double)(*rr));
// for (i=0; i<n; i++) {
// for (k=0; k<r; k++) {
// u = (x[i + k*n] - mu[j])/sigma[j];
// f[i + j*n] = 1.0;
// sum1 = 0.0;
// for (ii=0; ii<n; ii++) {
// sum2 = 0.0;
// for (kk=0; kk<r; kk++) {
// tmp = (u - x[ii+kk*n] + mu[j])/sigma[j];
// sum2 += exp(tmp * tmp * const1);
// }
// sum1 += z[ii + j*n] * sum2;
// }
// f[i + j*n] *= sum1 * const2;
// }
// }
// }
//}
//