speaker-confusion-model/code/models/posteriors.stan

// TODO
// use speech rates to set priors on truth_vocs
data {
    int<lower=1> n_classes; // number of classes

    // analysis data block
    int<lower=1> n_recs;
    int<lower=1> n_children;

    array[n_recs] int<lower=1> children;
    array[n_recs] real<lower=1> age;
    array[n_recs] int<lower=-1> siblings;
    array[n_children] int<lower=1> corpus;

    real<lower=0> recs_duration;

    // speaker confusion data block
    int<lower=1> n_clips;   // number of clips
    int<lower=1> n_groups; // number of groups
    int<lower=1> n_corpora;

    int<lower=0> n_validation;

    // actual speech rates
    int<lower=1> n_rates;
    int<lower=1> n_speech_rate_children;

    array [n_rates,n_classes] int<lower=0> speech_rates;
    array [n_rates] int group_corpus;
    array [n_rates] real<lower=0> durations;
    array [n_rates] real<lower=0> speech_rate_age;
    array [n_rates] int<lower=-1> speech_rate_siblings;
    array [n_rates] int<lower=1,upper=n_speech_rate_children> speech_rate_child;

    // parallel processing
    int<lower=1> threads;
}

transformed data {
    array[n_speech_rate_children] int<lower=1> speech_rate_child_corpus;

    array[n_children] int<lower=-1> child_siblings;
    array[n_speech_rate_children] int<lower=-1> speech_rate_child_siblings;
    int no_siblings = 0;
    int has_siblings = 0;

    for (k in 1:n_rates) {
        speech_rate_child_corpus[speech_rate_child[k]] = group_corpus[k];
    }

    for (k in 1:n_recs) {
        child_siblings[children[k]] = siblings[k];
    }

    for (c in 1:n_children) {
        if (child_siblings[c] == 0) {
            no_siblings += 1;
        }
        else if (child_siblings[c] > 0) {
            has_siblings += 1;
        }
    }

    for (k in 1:n_rates) {
        speech_rate_child_siblings[speech_rate_child[k]] = speech_rate_siblings[k];
    }
}

parameters {
    #include "blocks/behavior_model_parameters.stan"
    #include "blocks/human_annotations_parameters.stan"
}

model {
    //actual model
    #include "blocks/behavior_model_priors_uncentered.stan"
    #include "blocks/human_annotations_uncentered.stan"  
}

generated quantities {
    matrix<lower=0>[n_children,n_classes-1] mu_child_level;
    vector [n_children] child_dev_age;
    matrix<lower=0> [n_recs, n_classes] truth_vocs;

    for (c in 1:n_children) {
        child_dev_age[c] = normal_rng(alpha_dev, sigma_dev);
        // if there is sibling data
        if (child_siblings[c]>=0) {
            int distrib = child_siblings[c]>0?2:1;

            mu_child_level[c,1] = gamma_rng(
                alpha_corpus_level[distrib,1,corpus[c]],
                (alpha_corpus_level[distrib,1,corpus[c]]/(mu_corpus_level[1,corpus[c]]*exp(
                   child_siblings[c]>0?beta_sib_och:0 
                )))
            ); 
            mu_child_level[c,2:] = to_vector(gamma_rng(
                alpha_corpus_level[distrib,2:,corpus[c]],
                (alpha_corpus_level[distrib,2:,corpus[c]]./mu_corpus_level[2:,corpus[c]]*exp(
                   child_siblings[c]>0?beta_sib_adu/10:0 
                ))
            ))'; //'
        }
        // otherwise
        else {
            // assuming no sibling
            if (bernoulli_rng(p_sib)) {
                mu_child_level[c,1] = gamma_rng(alpha_corpus_level[2,1,corpus[c]], alpha_corpus_level[2,1,corpus[c]]/(mu_corpus_level[1,corpus[c]]*exp(beta_sib_och)));
                mu_child_level[c,2] = gamma_rng(alpha_corpus_level[2,2,corpus[c]], alpha_corpus_level[2,2,corpus[c]]/(mu_corpus_level[2,corpus[c]]*exp(beta_sib_adu/10)));
                mu_child_level[c,3] = gamma_rng(alpha_corpus_level[2,3,corpus[c]], alpha_corpus_level[2,3,corpus[c]]/(mu_corpus_level[3,corpus[c]]*exp(beta_sib_adu/10)));
            }
            else {
                mu_child_level[c,1] = gamma_rng(alpha_corpus_level[1,1,corpus[c]], alpha_corpus_level[1,1,corpus[c]]/(mu_corpus_level[1,corpus[c]]));
                mu_child_level[c,2] = gamma_rng(alpha_corpus_level[1,2,corpus[c]], alpha_corpus_level[1,2,corpus[c]]/(mu_corpus_level[2,corpus[c]]));
                mu_child_level[c,3] = gamma_rng(alpha_corpus_level[1,3,corpus[c]], alpha_corpus_level[1,3,corpus[c]]/(mu_corpus_level[3,corpus[c]]));
            }
        }
    }

    for (k in 1:n_recs) {
        real chi_mu = mu_pop_level[1]*exp(
            child_dev_age[children[k]]*age[k]/12.0/10.0+beta_dev*(mu_child_level[children[k],2]+mu_child_level[children[k],3]-mu_pop_level[3]-mu_pop_level[4])*age[k]/12.0/10.0
        );
        truth_vocs[k,1] = 1000*recs_duration*gamma_rng(alpha_child_level[1], alpha_child_level[1]/chi_mu);
        truth_vocs[k,2:] = 1000*recs_duration*to_vector(gamma_rng(alpha_child_level[2:], alpha_child_level[2:]./mu_child_level[children[k],:]'))'; //'
    }
}