Mixed data - parameter simplex

[James Savage]

Quick question - 

Let's say I have the first 4 values of a six-dimensional simplex, and I want to a) make estimates of alpha below, and b) generate predictions (par1 .. par3), which are parameters of the model. The restriction of course is that the parameters + data = 1. 

[data1, data2, data3, data4, par1, par2, par3]' ~ dirichlet(alpha)

Any suggestions on how to implement this? 

Currently, I'm building the simplex in the transformed parameters block, then calculating the log absolute determinant of the inverse transform (I think?). Is this just the log of the second expression on p487 of the 2.11 manual? To my mind, the model would then be something like

model {

// prior on alpha

// Expression which I suppose gives some structure to par1...par3?

target += my_log_abs_jac_inv_trans_calculated_elsewher;

//likelihood

my_simplex ~ dirichlet(alpha);

}

Does that seem right? I'm doing this at the moment and it's not working. So I'm trying to work out whether my approach to implementing is wrong, or whether I'm calculating the wrong log abs jac... 

Cheers, 

Jim

[Daniel Lee]

Hmm... it seems like you should have something like:

data {

  vector[4] x;

}

parameters {

  simplex[2] alpha;

}

transformed parameters {

  simplex[6] beta;

  beta[1:4] = x;

  beta[5] = alpha[1] / (1 - sum(x));

  beta[6] = alpha[2] / (1 - sum(x));

}

Since x is data, I don't think you need a Jacobian here. That's beta[5] and beta[6] are just alpha[1] and alpha[2] divided by a constant.

Does that make sense?

(Of course, you could simplify in this case by using real<lower = 0, upper = 1> alpha and using alpha and (1 - alpha)).

Daniel

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

Nice formulation! You could also code this in one go as:

beta[5:6] = alpha / (1 - sum(x));

I'm still getting used to this slicing magic myself.

The transform of alpha to beta[5:6] is linear in the data,
so no need for a Jacobian.

- Bob

> 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.
>
>

> --
> 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.

[James Savage]

Very neat indeed, and yes Bob- magical. 

Now my problem is a bit more difficult. My data actually look more like

Y = 

[.4, .2, NA, NA;

 .5, .1, .1, NA;

 .2, NA, NA, NA]

where each row is a simplex by design, but where different rows have different

number of missing rows (that will be imputed using the parameters above). 

Each row has at least one missing value (to be clear, I'm not giving Stan missing 

values; I'm giving it an accompanying position matrix for the missing observations. 

So I'm not sure I can use the elegant solution above. 

Cheers, 

Jim

[Daniel Lee]

You can't do it easily in Stan without the ragged data structure, but you could code that up if you wanted. 

I think it'd be hard to write a general program to handle all sorts of missingness, so multiple data sets with one program, but for one data set you can write it out in Stan. 

Daniel

[James Savage]

Thanks Daniel - 

Would coding it up involve something like the approach in my original question? 

Jim

[Bob Carpenter]

You just need to build up each row of Y separately.
There's no good way to generalize here because each simplex
has to be declared separately to get the size right.

You could pass in the data as a matrix as long as you
don't mind some unused values, but you'd need something
else like

int row_data_sizes[M]; // num data elements in row m

to say how many rows each has, and given that each simplex
is separate, you'd want to break the build-up into a function
to avoid code duplication.

- Bob

[Michael Betancourt]

Seems like it would be much easier to do some preprocessing
to order all of the NAs to the left or the right and split up the data
into those simplices with one missing element, two, … etc. Then
you can built the Stan program in a straightforward fashion.

[James Savage]

Mike - you've just saved me a lot of time and heartache. Thanks. 

[Michael Betancourt]

:thumbsup: