-
Notifications
You must be signed in to change notification settings - Fork 30
Crosstab to PMF
b-k edited this page Apr 14, 2013
·
1 revision
A version of this used to be in the Apophenia library, but it got cut.
#include <apop.h>
/** Convert a crosstab to a PMF model, which you can use to make random draws, as a prior
for Bayesian updating (via \ref apop_update), and so on.
\param d An input crosstab, expressed either in the matrix of an apop_data set.
\return A PMF model which has a single line for each nonzero line of the crosstab.
The model points to a data set generated here, which you will need to
apop_data_free whenever you apop_model_free this model.
*/
apop_model *apop_crosstab_to_pmf (apop_data *d){
Apop_stopif(!d, return NULL, 0, "NULL input.");
Apop_stopif(!d->matrix, return NULL, 0, "You gave me an input set with no matrix data.");
int tsize = d->matrix->size1 * d->matrix->size2;
apop_data *outdata =apop_data_alloc(tsize, 2);
outdata->weights = gsl_vector_alloc(tsize);
size_t ctr = 0;
double x;
for (size_t i=0; i < d->matrix->size1; i++)
for (size_t j=0; j < d->matrix->size2; j++)
if ((x = gsl_matrix_get(d->matrix, i, j))) {
apop_data_set(outdata, ctr, 0, i);
apop_data_set(outdata, ctr, 1, j);
gsl_vector_set(outdata->weights, ctr++, x);
}
return apop_estimate(outdata, apop_pmf);
}
int main(){ //test script with magic square.
apop_model *pmf = apop_crosstab_to_pmf(
apop_data_fill(apop_data_alloc(3,3),
8, 1, 6,
3, 5, 7,
4, 9, 2));
int draw_ct=1e5;
apop_data *draws= apop_data_alloc(draw_ct, 2);
gsl_rng *r = apop_rng_alloc(213);
for (int i=0; i < draw_ct; i++){
Apop_row(draws, i, onerow);
apop_draw(onerow->data, r, pmf);
}
apop_data *sum = apop_data_summarize(draws);
assert(fabs(apop_data_get(sum, 0, .colname="mean")-1) < 1e-3);
assert(fabs(apop_data_get(sum, 1, .colname="mean")-1) < 1e-3);
}