Supervised Classification

[Karl Möller]

Hi

I am new to Bayesian modelling and MCMC / Stan so please forgive all rookie mistakes... 

I am trying to construct a minimal working example of a Bayesian classifier for labeled data that should work for my requirements. Let us assume that 1D training samples X = {x_i} from 

every class (label) c_i are generated by some normal distribution with unknown \mu_i and \sigma_i specific to class c_i. So for a new test pattern x_new we know that:

P(c_i | x_new, X) = P(x_new | c_i, X) P(c_i) / (\sum_j (P(x_new | c_j, X) P(c_j))  

where we can further write P(x_new | c_i, X) = \int d \theta P(x_new | c_i, \theta_i) P(\theta_i | c_i, X) where \theta_i = (\mu_i, \sigma_i)

So to get distribution P(\theta_i|c_i, X) is pretty trivial using the following Stan code:

data {

    int<lower=0> N; // num of training set samples 

    int<lower=1> M; // num of classes

    real X[N];  //training data

    int<lower=1, upper=M> C[N]; // class for sample n

   real x_new; //new test pattern to classify

}

transformed data {}

parameters {

    real mu[M]; 

    real sigma[M]; 

    real x_pred[M];

}

transformed parameters {}

model {

    //likelihood

    for (n in 1:N) 

        X[n] ~ normal(mu[C[n]], sigma[C[n]]) 

       

    //predictive posterior

    for (m in 1:M)

       x_pred[m] ~ normal(mu[m], sigma[m]) 

        

}

generated quantities {

}

I probably don't need the predictive posterior, but just used it as a sanity check.

What I do not understand how to do in Stan is the marginalisation (integral/sum) over \theta to get the expected likelihood (and evidence

follows trivially) and then find P(c_i | x_new). What I can do is 

generated quantities{

    real x_new_log_likelihood[M];

    for (m in 1:M)

        x_new_log_likelihood[m] <- normal_log(x_new, mu[m], sigma[m])

}

And then in Python calculated the mean of 'x_new_log_likelihood' for every class (actually doing the required transformations because log \sum is not equal to \sum log) and plug in some priors for P(c_i) and just calculate

P(c_i | x_new) in Python directly. Is there a better way I can do all of this in Stan directly? I hope my explanation is sufficient. I have looked at the 

Naive Bayes example in the manual and I think it should be the same as what I want, but the classification for a new test pattern is not part of the example if I understand

it correctly. 

Regards

Karl

[Bob Carpenter]

Stan will do this all for you pretty much automatically.
The manual has some examples of how to do the predictive
distributions within Stan (that is, you use MCMC for
the integral).

As to the normalized category probabilities, you just
need to normalize:

vector[K] lp;
simplex[K] log_theta;

log_theta = lp - log_sum_exp(lp);

This assumes lp[k] = p(y | k) * p(k),
where p(y|k) is probably of observation given category k
and p(k) is the marginal probability of categories.

Then log_theta[k] is log Pr[C = k | y].

exp(theta) = exp(lp) / sum(exp(lp))

- Bob

> On Jul 20, 2016, at 11:27 AM, Karl Möller <ka...@excors.net> wrote:
>
> Not doing too well for my 1st post to the group, sorry. Correction:

>
> generated quantities{
>
> real x_new_log_likelihood[M];
> for (m in 1:M)
> x_new_log_likelihood[m] <- normal_log(x_new, mu[m], sigma[m])
>
> }
>

> And everywhere the word "score" pops up, I mean X (stripped out this from something I am trying). I hope these all the errors...
>
> Regards
> Karl

>
> On Wednesday, July 20, 2016 at 5:11:02 PM UTC+2, Karl Möller wrote:
> Hi
>
> I am new to Bayesian modelling and MCMC / Stan so please forgive all rookie mistakes...
>
> I am trying to construct a minimal working example of a Bayesian classifier for labeled data that should work for my requirements. Let us assume that 1D training samples X = {x_i} from

> every class (label) c_i are generated by some normal distribution with unknown \mu_i and \sigma_i specific to class c_i. So for a new text pattern x_new we know that:

>
> P(c_i | x_new, X) = P(x_new | c_i, X) P(c_i) / (\sum_j (P(x_new | c_j, X) P(c_j))
>
> where we can further write P(x_new | c_i, X) = \int d \theta P(x_new | c_i, \theta_i) P(\theta_i | c_i, X) where \theta_i = (\mu_i, \sigma_i)
>
> So to get distribution P(\theta_i|c_i, X) is pretty trivial using the following Stan code:
>
> data {
> int<lower=0> N; // num of training set samples
> int<lower=1> M; // num of classes
> real X[N]; //training data

> int<lower=1, upper=M> C[N]; // class for score n

> real x_new; //new test pattern to classify
>
> }
> transformed data {}
> parameters {
> real mu[M];
> real sigma[M];
> real x_pred[M];
> }
> transformed parameters {}
> model {
> //likelihood
> for (n in 1:N)

> scores[n] ~ normal(mu[C[n]], sigma[C[n]])

>
> //predictive posterior
> for (m in 1:M)
> x_pred[m] ~ normal(mu[m], sigma[m])
>
> }
> generated quantities {
> }
>
> I probably don't need the predictive posterior, but just used it as a sanity check.
> What I do not understand how to do in Stan is the marginalisation (integral/sum) over \theta to get the expected likelihood (and evidence
> follows trivially) and then find P(c_i | x_new). What I can do is
>
> generated quantities{
>
> real x_new_log_likelihood[M];
> for (m in 1:M)

> new_score_log_likelihood[m] <- normal_log(x_new)
>
> }
>
> And then in Python calculated the mean of 'new_score_log_likelihood' for every class (actually doing the required transformations because log \sum is not equal to \sum log) and plug in some priors for P(c_i) and just calculate

> P(c_i | x_new) in Python directly. Is there a better way I can do all of this in Stan directly? I hope my explanation is sufficient. I have looked at the
> Naive Bayes example in the manual and I think it should be the same as what I want, but the classification for a new test pattern is not part of the example if I understand
> it correctly.
>
> Regards
> Karl
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Karl Möller]

Hi Bob

Thanks for the reply. I looked at some of the examples in the manual, but I am still a bit lost. I am not

sure how to parametrise p(y | k)  (where y is the new test sample) in the model block, or am I missing something and 

this should be a generated quantity? For something like linear regression you can parametrise the mean of a normal, as in plug in

the new test sample x_new:  y_new ~ normal(x_new * beta, sigma); and declare y_new as a parameter.

But I can't do something similar in my case?

Sorry if it's really a stupid question, I have been trawling the manual and the mailing list for a while, maybe I am just missing the right keywords.

Regards

Karl

[Karl Möller]

Hi

Is this what you had in mind:

generated quantities{

  vector[M] lp; 

  vector[M] log_theta;  //should not be a simplex, log probabilities don't sum to 1

  simplex[M] theta; // P(C = m | x_new)

  

  for (m in 1:M)

      lp[m] <- 1/M + normal_log(x_new, mu[m], sigma[m]); //just assuming the class priors are 1/M

      

  for (m in 1:M)

      log_theta[m] <- lp[m] - log_sum_exp(lp);

      

  theta <- exp(log_theta);

    

}

I don't see how this can be done as part of the model block (do I even want that?) So here I get a probability distribution P(C = m | x_new) and then the expected value of

that distribution with a loss function (which I can't do in Stan, but in Python or R) gives me a single number?

Regards

Karl

[Karl Möller]

Mmm no, P(C = m | x_new) should be a single value right?

[Bob Carpenter]

That's right other than that it should be -log(M)
rather than 1/M, because everythng's on the log scale.

theta is a vector. And you can do the same thing
for more predictions.

Basically, it's just Bayes's rule coded in Stan.
I thought we had a "naive Bayes" classifier example, but
we don't. I'll add one:

https://github.com/stan-dev/stan/issues/1238#issuecomment-234316954

- Bob

P.S. As of Stan 2.10, the syntax is

normal_lpdf(x_new | ...)

rather than

normal_log(x_new, ...)

[the "lpdf" is for "log probability density function"].

Also, we're now using "=" instead of "<-" for assignment.

We'll still support the old syntax for quite a while.

[Karl Möller]

Ah yes, -log(M) obviously. More predictions I I think understand, what I am not fully understand is why I get a vector if I only predict for a single scalar. Or I think I do, but maybe 

it is not exactly what I want. If I calculate P(C=k | x_new) mathematically, that should be a single value right? As far as I understand, the generated_quantities 

block gets executed on every draw from the posterior p(mu, sigma), so after a specific mu_i, sigma_i is sampled, a value for theta is generated. So after Stan runs,

I get as many different thetas as there are different mu_i, sigma_i. So if I want my single value, I should take the mean of the thetas in Python after Stan has run?

Am I run or am I missing something? Also thanks for telling me about the new syntax, I have been looking at the PyStan manual which is still 2.9.

Regards

Karl

[Bob Carpenter]

Pr[C = k | x_new, x, y] is a single value, but then it's a vector
when you consider all k values.

The final output is a sample from the posterior consisting of
a sequence of draws. You can get a single expectation by
averaging, and you can also get uncertainty in that estimate by
looking at quantiles. It's standard posterior predictive inference,
which averages over the posterior.

Pr[C = k | x_new, x, y]

= INTEGRAL_theta Pr[C = k | x_new, theta] p(theta | x, y) d.theta

where theta are your parameters.

- Bob

[Karl Möller]

Ah yes thanks it makes sense to me now, I guess I was stuck thinking of it the way I originally wrote the probability expansion (which seems equivalent).

I wanted to marginalise out theta when calculating Pr[x = x_new | x, theta] and then calculate the class assignment posteriors Pr[C = k | x_new, x] using that result.