Distribution for modeling a scale

[Ross Boylan]

I have a variable that takes on values in the range 1-4. I was thinking
of modeling this as a beta distribution, but since there are a lot of
cases at the extremes that seems problematic. I'm not sure the beta
function in stan will even accept such inputs and produce estimates. I
could fake it by pulling the values in a bit, but I'm not sure that's a
good idea.

I'd appreciate any suggestions on how to fit such outcomes. I don't
think it's a big deal if the fit doesn't reproduce the detailed features
of the distribution for my purposes, though if this lack of fit leads to
estimation problems that would matter.

"Detailed" features refers to that fact that most of the observations
are exactly integers, and all the values are evenly divisible by 1/3.

My main focus is on how covariates shift the distribution. One approach
might be to treat the outcome as categorical or ordinal, but that seems
likely to create a lot more parameters.

The distribution also kind of looks like a truncated normal if one uses
1 unit size bins.

If you're curious, here's the basic data:

> table(tfat$FrndPartRace)

1 1.33333333333333 1.66666666666667 2
162 55 83 156
2.33333333333333 2.66666666666667 3 3.33333333333333
71 73 93 19
3.66666666666667 4
5 29

Thanks.
Ross

[Andrew Gelman]

HI, my quick hack would be to take the logit of the variable after transforming to a 0.5 to 4.5 scale, that is, work with logit((y - 0.5)/4). This will make it continuous and will spread out the extremes a bit. Of course it won't give correct predictions but otherwise it could do what you need.
If it works (or if it doesn't), please let me know!
A

> --
> 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.
> For more options, visit https://groups.google.com/groups/opt_out.

[Ross Boylan]

On 1/31/2014 3:33 PM, Andrew Gelman wrote:
> HI, my quick hack would be to take the logit of the variable after transforming to a 0.5 to 4.5 scale, that is, work with logit((y - 0.5)/4). This will make it continuous and will spread out the extremes a bit. Of course it won't give correct predictions but otherwise it could do what you need.
> If it works (or if it doesn't), please let me know!
> A

I get this error:

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model 'stan_model' with error message:
EXPECTATION FAILURE LOCATION: file=input; line=14, column=21

z <- logit((y-0.5)/4);
^-- here


DIAGNOSTIC(S) FROM PARSER:
no matches for function name="subtract"
arg 0 type=real[1]
arg 1 type=real
available function signatures for subtract:
0. subtract(vector, vector) : vector
1. subtract(row vector, row vector) : row vector
2. subtract(matrix, matrix) : matrix
3. subtract(vector, real) : vector
4. subtract(row vector, real) : row vector
5. subtract(matrix, real) : matrix
6. subtract(real, vector) : vector
7. subtract(real, row vector) : row vector
8. subtract(real, matrix) : matrix
expression is ill formed
Parser expecting: <eps>


The relevant block is:
transformed data {
real z[N];
z <- logit((y-0.5)/4);
}

The error seems to indicate stan can't subtract a constant from a
vector, but that seems unlikely.

Ross

P.S. I've been running a truncated normal model, but it had poor n_eff
(3 for 200 iterations) in tests and is taking a long time in a
production run.

[Ross Boylan]

On 1/31/2014 4:59 PM, Ross Boylan wrote:
> On 1/31/2014 3:33 PM, Andrew Gelman wrote:
>> HI, my quick hack would be to take the logit of the variable after
>> transforming to a 0.5 to 4.5 scale, that is, work with logit((y -
>> 0.5)/4). This will make it continuous and will spread out the
>> extremes a bit. Of course it won't give correct predictions but
>> otherwise it could do what you need.
>> If it works (or if it doesn't), please let me know!
>> A
>

I programmed around the vectorization problem, but got this on a test
run (200 iterations):

> system.time(rb <- stan(model_code = stan_model, data=stanDat, seed=seed,
+ iter=200, chains=1))

TRANSLATING MODEL 'stan_model' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'stan_model' NOW.
SAMPLING FOR MODEL 'stan_model' NOW (CHAIN 1).
Iteration: 1 / 200 [ 0%] (Warmup)
Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Location parameter[0] is -nan:0, but must be finite!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Location parameter[0] is -nan:0, but must be finite!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Iteration: 200 / 200 [100%] (Sampling)
Elapsed Time: 5.49 seconds (Warm-up)
3.62 seconds (Sampling)
9.11 seconds (Total)

user system elapsed
43.507 1.096 49.298
> summary(summary(rb)$summary[, c('n_eff', 'Rhat')])
n_eff Rhat
Min. : 2.737 Min. :0.9899
1st Qu.: 5.525 1st Qu.:1.0215
Median : 8.201 Median :1.1171
Mean : 9.599 Mean :1.2534
3rd Qu.:11.767 3rd Qu.:1.3808
Max. :83.189 Max. :3.0188


[BTW: becuase is misspelled in the error message.]

Here's the model:

data {
int<lower=0> N; // number of respondents
int<lower=0> J; //number of covariates
int<lower=0> NSD; // number of classes of beta for sd
int<lower=1,upper=NSD> SDClass[J]; // beta[j, k] has sigma_beta[SDClass[j]]
int<lower=0> NCLUSTER; //number of seeds
int iCluster[N]; // cluster id
matrix[N, J] x; // covariates
real<lower=1,upper=4> y[N]; // outcome

}
transformed data {
real z[N];

for (i in 1:N)
z[i] <- logit((y[i]-0.5)/4);
}
parameters {
vector[J] beta_raw; // standard covariate coefficients
vector[NCLUSTER] alpha_raw; // cluster effects
real<lower=0> sigma;
real<lower=0> sigma_alpha;
real<lower=0> sigma_beta[NSD];
}
transformed parameters {
vector[J] beta;
vector[NCLUSTER] alpha;
for(ir in 1:J) {
beta[ir] <- beta_raw[ir]*sigma_beta[SDClass[ir]];
}
alpha <- alpha_raw*sigma_alpha;
}
model {
vector[N] beta_x;
sigma ~ cauchy(0, 5);
sigma_alpha ~ cauchy(0, 5);
sigma_beta ~ cauchy(0, 5);
beta_raw[1] ~ uniform(-5, 5);
for (j in 2:J) {
beta_raw[j] ~ normal(0, 1);
}
alpha_raw ~ normal(0, 1);
for (n in 1:N) {
beta_x[n] <- (x[n] * beta)+alpha[iCluster[n]];
}
z ~ normal(beta_x, sigma);
}

I haven't updated stan for at least a few weeks.
stan_version() gives 2.0.1

The other model I'm running had something similar:

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Location parameter is -nan:0, but must be finite!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

I may have some programming error, though I'm not seeing it.



Ross

[Ross Boylan]

For a more production-sized run I still got only 2 errors, this time
about scale of 0:

> system.time(sflist <- mclapply(1:4, function(i) stan(fit=rb, data=stanDat, iter=4000, chains=1,
+ seed=seed, chain_id=i, refresh=-1,
+ warmup=500, thin=4),
+ mc.cores=4)
+ )

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

user system elapsed
156.973 0.984 59.846
> myfit <- sflist2stanfit(sflist)
> summary(summary(myfit)$summary[, c('n_eff', 'Rhat')])
n_eff Rhat
Min. : 65.8 Min. :0.9996
1st Qu.: 407.3 1st Qu.:1.0031
Median : 731.1 Median :1.0058
Mean : 749.3 Mean :1.0078
3rd Qu.:1059.7 3rd Qu.:1.0107
Max. :1761.0 Max. :1.0664

[Bob Carpenter]

You get the error listed below because arithmetic is
not supported on arrays, only vectors. Your y is
a one-dimensional array (real[1]) and 0.5 is real, and
real[1] - real is not allowed. You can see the list of
signatures that are allowed in the warning message.

Any suggestions on how to make that diagnostic message
clearer? I think what would probably help is to back-translate
from the function name ("-" is replaced with the function "subtract").

- Bob

[Bob Carpenter]

On 2/1/14, 2:39 AM, Ross Boylan wrote:
> For a more production-sized run I still got only 2 errors, this time about scale of 0:

^^^^^^
...

> Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
> Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
> If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the
> sampler is fine,
> but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
>

If you're talking about the "Informational Message" reports, those are NOT
errors, they are warnings.

Did you read the part that says

"If this warning occurs sporadically, such as for highly constrained variable types like covariance
matrices, then the sampler is fine,"

and if so, could you let us know how we can make it clearer that
this is just an informational message and not an error?

You're not the only one who's interpreted these messages as errors.
I think part of the confusion is that the embedded message says "error",
whereas we've tried to keep the referencing part clear that it's a warning.

But we want to keep them because if the rejections happen repeatedly,
then it can be a problem for sampling efficiency (up to the point where
if they happen in every step of every iteration the sampler won't move
at all).

I'm also wondering if users are familiar with the notion of Metropolis
rejection or if that's too much MCMC jargon.

- Bob

>
> > system.time(sflist <- mclapply(1:4, function(i) stan(fit=rb, data=stanDat, iter=4000, chains=1,
> + seed=seed, chain_id=i, refresh=-1,
> + warmup=500, thin=4),
> + mc.cores=4)
> + )
>

> Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
> Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
> If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the
> sampler is fine,
> but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
>
>
>
>

[Bob Carpenter]

It's also multimodal around the even integers.

If you really want to capture this, I don't see any way to do it
other than an ordinal logit which could have the cutpoints vary in
size.

- Bob

>> table(tfat$FrndPartRace)
>
> 1 1.33333333333333 1.66666666666667 2
> 162 55 83 156
> 2.33333333333333 2.66666666666667 3 3.33333333333333
> 71 73 93 19
> 3.66666666666667 4
> 5 29

[Ross Boylan]

On Sat, 2014-02-01 at 11:08 +0100, Bob Carpenter wrote:
> On 2/1/14, 2:39 AM, Ross Boylan wrote:
> > For a more production-sized run I still got only 2 errors, this time about scale of 0:
> ^^^^^^
> ...
> > Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
> > Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
> > If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the
> > sampler is fine,
> > but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
> >
>
> If you're talking about the "Informational Message" reports, those are NOT
> errors, they are warnings.
>
> Did you read the part that says
>
> "If this warning occurs sporadically, such as for highly constrained variable types like covariance
> matrices, then the sampler is fine,"
>
> and if so, could you let us know how we can make it clearer that
> this is just an informational message and not an error?

I read it; that's why I said "only", because the message indicated that
2 wasn't really a problem. It still registered as an error with me
because I had bounded the parameter at 0 and so it didn't seem it should
generate 0 as a candidate value. Since 0 <= 0 that expectation is
probably unreasonable, although I thought, as a matter of implementation
if not principle, that variables would never take on the value at the
bounds.

>
> You're not the only one who's interpreted these messages as errors.
> I think part of the confusion is that the embedded message says "error",
> whereas we've tried to keep the referencing part clear that it's a warning.

The "error" does contribute to the sense there's been an error, even if
the rest of the message is quite explicit.

>
> But we want to keep them because if the rejections happen repeatedly,
> then it can be a problem for sampling efficiency (up to the point where
> if they happen in every step of every iteration the sampler won't move
> at all).
>
> I'm also wondering if users are familiar with the notion of Metropolis
> rejection or if that's too much MCMC jargon.

That seems pretty basic for a user of the package, though I suppose if
it takes the world by storm there will be users who don't know much
about how it gets the results.

Ross

[Ross Boylan]

On Sat, 2014-02-01 at 10:58 +0100, Bob Carpenter wrote:
> You get the error listed below because arithmetic is
> not supported on arrays, only vectors. Your y is
> a one-dimensional array (real[1]) and 0.5 is real, and
> real[1] - real is not allowed. You can see the list of
> signatures that are allowed in the warning message.
>
> Any suggestions on how to make that diagnostic message
> clearer? I think what would probably help is to back-translate
> from the function name ("-" is replaced with the function "subtract").
>

I think the meaning of real[1] is less clear than the fact that
"subtract"="-". In this case the message could say "real array" or
perhaps "real array[1]" or even "1-dimensional real array", but that
might be too verbose for general use.

My problem was not that I didn't understand the message, but that I
didn't believe it :) That is, I thought combining arrays and scalars
must be supported, and so either I had made an error, the system had
made an error, or I had misunderstood the message.

I think it would be nice if scalars and arrays could be combined. I
can't think of other examples where they can't (e.g., R, python numeric,
C++ math libraries other than stan, SAS IML).

Does this mean that, unless I have a good reason, I should prefer
vector[n] x to real x[n]?

[Bob Carpenter]

Thanks for the response!

We'd like to have lower=0 be exclusive, but floating point
just doesn't cooperate. Instead, the exp(x) inverse transform
underflows to 0.

Then, when the value gets used as say, a scale for a normal
density, the density function throws an error. That gets caught
further up the chain and turned into a Metropolis rejection and
the error message wrapped in that verbose and confusing message.

I suppose we could consider adding epsilon to the underflows and
doing the same for gradients, but that's much harder for constraints
like covariance matrices, where our big problems are with this
message and with inverse transforms.

- Bob

[Andrew Gelman]

Bob
I don't think the exact values really matter, at least that was the impression I got from Ross's email. And with a judicious logit transform, you get rid of the modes at the ends because they get spread out. Hence my suggestion of logit with boundaries at .5 and 4.5.
A

[Ross Boylan]

On Sat, 2014-02-01 at 19:54 +0100, Andrew Gelman wrote:
> Bob
> I don't think the exact values really matter, at least that was the impression I got from Ross's email. And with a judicious logit transform, you get rid of the modes at the ends because they get spread out. Hence my suggestion of logit with boundaries at .5 and 4.5.
> A

I don't think the exact values matter, mostly because in the larger
system the outcome variable here becomes a covariate in downstream
predictions. And the specifications of those predictions use simple
terms like beta*FrndPartRace, albeit in nonlinear models.

I don't see how the transform changes the fact the distribution is
discrete and multimodal.

The measurements are of an underlying concept that is continuous, though
probably only ordinal. The observed measurement come from averaging
three items coded 1-4 (disagree - agree scale).

My alternate model as a truncated normal finally finished. The results
were not pretty. So thanks to Andrew for his suggestion. In case
anyone is interested, here are the truncated normal results:

> system.time(sflist <- mclapply(1:6, function(i) stan(fit=rb, data=stanDat, iter=4000, chains=1,
+ control=list(adapt_delta=0.99),

+ seed=seed, chain_id=i, refresh=-1,
+ warmup=500, thin=4),

+ mc.cores=6)

+ )

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:

Error in function stan::prob::normal_log(N4stan5agrad3varE): Location parameter is -nan:0, but must be finite!

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

user system elapsed

37652.281 1.804 13591.803 # 3.8 hours wall-time
Warning messages:
1: In selectChildren(ac, 1) : error 'Interrupted system call' in select
2: In selectChildren(ac, 1) : error 'Interrupted system call' in select
3: In selectChildren(ac, 1) : error 'Interrupted system call' in select
4: In selectChildren(ac, 1) : error 'Interrupted system call' in select
5: In selectChildren(ac, 1) : error 'Interrupted system call' in select

> b5 <- sflist2stanfit(sflist)
> #alternative from Sean O'Riordain
> summary(summary(b5)$summary[, c('n_eff', 'Rhat')])
n_eff Rhat
Min. : 3.017 Min. : 1.006
1st Qu.: 5.279 1st Qu.: 1.155
Median : 8.019 Median : 1.296
Mean : 24.437 Mean : 1.644
3rd Qu.: 17.629 3rd Qu.: 1.561
Max. :478.178 Max. :13.870

I'm not really sure what to make of the interrupted system call
messages, which have not occurred in other, similar, calls. Presumably
they are coming from R's parallel library and the underlying OS, rather
than stan. Perhaps they have something to do with this being a first
time experiment in which I ran more threads than physical cores (4),
though fewer than the virtual, hyperthreaded, cores.

At any rate, the computation did return 6 full sets of answers. I'm a
little suprised n_eff isn't a bit better given the individual chains:
> lapply(sflist, function(x) summary(summary(x)$summary[,c("n_eff", "Rhat")]))
[[1]]
n_eff Rhat
Min. : 2.402 Min. :0.9989
1st Qu.: 4.330 1st Qu.:1.0421
Median : 6.476 Median :1.2018
Mean : 8.464 Mean :1.3598
3rd Qu.:10.006 3rd Qu.:1.4946
Max. :62.820 Max. :4.2685

[[2]]
n_eff Rhat
Min. : 31.33 Min. :0.9989
1st Qu.:113.66 1st Qu.:0.9996
Median :147.23 Median :1.0021
Mean :155.63 Mean :1.0059
3rd Qu.:192.01 3rd Qu.:1.0080
Max. :466.75 Max. :1.0594

[[3]]
n_eff Rhat
Min. : 7.142 Min. :0.9989
1st Qu.: 26.530 1st Qu.:1.0023
Median : 37.334 Median :1.0132
Mean : 40.613 Mean :1.0330
3rd Qu.: 50.733 3rd Qu.:1.0436
Max. :138.818 Max. :1.3402

[[4]]
n_eff Rhat
Min. : 2.344 Min. :0.9989
1st Qu.: 3.312 1st Qu.:1.1151
Median : 4.419 Median :1.4453
Mean : 5.357 Mean :1.5919
3rd Qu.: 6.519 3rd Qu.:1.8981
Max. :22.484 Max. :4.0677

[[5]]
n_eff Rhat
Min. : 2.359 Min. :0.9989
1st Qu.: 3.049 1st Qu.:1.1413
Median : 4.082 Median :1.5455
Mean : 5.091 Mean :1.6899
3rd Qu.: 5.993 3rd Qu.:2.0985
Max. :22.764 Max. :4.1930

[[6]]
n_eff Rhat
Min. : 23.94 Min. :0.9989
1st Qu.: 83.62 1st Qu.:0.9999
Median :112.07 Median :1.0038
Mean :114.83 Mean :1.0095
3rd Qu.:141.93 3rd Qu.:1.0121
Max. :271.34 Max. :1.1148

I'd expect n_eff to be more like the sum of the individual chain
n_eff's, or at least > than the largest single value. Since chain 2 has
min n_eff of 31, how can all the chains together have smaller minimum
n_eff than chain 2 alone?

Model:

data {
int<lower=0> N; // number of respondents
int<lower=0> J; //number of covariates
int<lower=0> NSD; // number of classes of beta for sd
int<lower=1,upper=NSD> SDClass[J]; // beta[j, k] has sigma_beta[SDClass[j]]
int<lower=0> NCLUSTER; //number of seeds
int iCluster[N]; // cluster id
matrix[N, J] x; // covariates
real<lower=1,upper=4> y[N]; // outcome
}

parameters {
vector[J] beta_raw; // standard covariate coefficients
vector[NCLUSTER] alpha_raw; // cluster effects
real<lower=0> sigma;
real<lower=0> sigma_alpha;
real<lower=0> sigma_beta[NSD];
}
transformed parameters {
vector[J] beta;
vector[NCLUSTER] alpha;
for(ir in 1:J) {
beta[ir] <- beta_raw[ir]*sigma_beta[SDClass[ir]];
}
alpha <- alpha_raw*sigma_alpha;
}
model {
vector[N] beta_x;
sigma ~ cauchy(0, 5);
sigma_alpha ~ cauchy(0, 5);
sigma_beta ~ cauchy(0, 5);
beta_raw[1] ~ uniform(-5, 5);
for (j in 2:J) {
beta_raw[j] ~ normal(0, 1);
}
alpha_raw ~ normal(0, 1);
for (n in 1:N) {
beta_x[n] <- (x[n] * beta)+alpha[iCluster[n]];

//truncated distributions don't vectorize, probably
y[n] ~ normal(beta_x[n], sigma) T[1, 4];
}

Ross

[Andrew Gelman]

This is getting away from Stan so maybe it's not worth continuing this discussion on the models, but just briefly:
1. As Jennifer and I talk about in our book, it can be fine to fit a continuous model to discrete data, as long as the goal is not to predict individual values.
2. When you talk about multi modality, there are 2 things going on. First, the data are discrete and lumpy so have lots of modes all over. Second, there are lots of data near the beginning and end of the scale, so even a smoothed distribution of the data has 2 or 3 modes. Tranforming to the logit scale will spread out the modes at the ends, indeed in some cases can spread out such modes so much that they are no longer modes.
A

[Bob Carpenter]

> On Feb 1, 2014, at 10:48 PM, Ross Boylan wrote:

...

>> I don't see how the transform changes the fact the distribution is
>> discrete and multimodal.

Andrew answered this one. Transforms affect modes. It's
why maximum likelihood estimates are sensitive to parameterization.

>> Warning messages:
>> 1: In selectChildren(ac, 1) : error 'Interrupted system call' in select

...

>> I'm not really sure what to make of the interrupted system call
>> messages, which have not occurred in other, similar, calls.Presumably
>> they are coming from R's parallel library and the underlying OS, rather
>> than stan. Perhaps they have something to do with this being a first
>> time experiment in which I ran more threads than physical cores (4),
>> though fewer than the virtual, hyperthreaded, cores.

No idea, but they're not coming from Stan.

>> At any rate, the computation did return 6 full sets of answers. I'm a
>> little suprised n_eff isn't a bit better given the individual chains:

...

>> I'd expect n_eff to be more like the sum of the individual chain
>> n_eff's, or at least > than the largest single value. Since chain 2 has
>> min n_eff of 31, how can all the chains together have smaller minimum
>> n_eff than chain 2 alone?

That's what you'd get with more tradition n_eff calculations,
but Stan's is more conservative, using both within-chain and
cross-chain measures. You can see the definition in the Stan manual
and in the new edition of Gelman et al.'s Bayesian Data Analysis.

What you'll see is that higher R-hat will cause lower effective sample
sizes. As R-hat goes to 1, the n_eff we calculate will converge to
calculating each chain independently and averaging.

- Bob

[Bob Carpenter]

On 2/1/14, 7:31 PM, Ross Boylan wrote:
> On Sat, 2014-02-01 at 10:58 +0100, Bob Carpenter wrote:
>> You get the error listed below because arithmetic is
>> not supported on arrays, only vectors. Your y is
>> a one-dimensional array (real[1]) and 0.5 is real, and
>> real[1] - real is not allowed. You can see the list of
>> signatures that are allowed in the warning message.
>>
>> Any suggestions on how to make that diagnostic message
>> clearer? I think what would probably help is to back-translate
>> from the function name ("-" is replaced with the function "subtract").
>>
> I think the meaning of real[1] is less clear than the fact that
> "subtract"="-". In this case the message could say "real array" or
> perhaps "real array[1]" or even "1-dimensional real array", but that
> might be too verbose for general use.

That's a good idea. I don't want to be too verbose, but we can
do better than real[1]. I thought the
1 would be confusing. Also, real[2] looks like an array of size
2, not a 2-dimensional array.

I created an issue to fix this in the future:

https://github.com/stan-dev/stan/issues/529

> My problem was not that I didn't understand the message, but that I
> didn't believe it :) That is, I thought combining arrays and scalars
> must be supported, and so either I had made an error, the system had
> made an error, or I had misunderstood the message.
>
> I think it would be nice if scalars and arrays could be combined. I
> can't think of other examples where they can't (e.g., R, python numeric,
> C++ math libraries other than stan, SAS IML).

Eigen's an example, and it's the C++ math lib we use :-)
You can wrap an array into a matrix or vice-versa, but you can't do
matrix arithmetic on arrays.

The reason we can't just automatically coerce back and forth is
that vectors are either row vectors (1 x K) or column vectors (K x 1),
so there's not an unambiguous way to cast an array without context.

We do allow them to mix in contexts like vectorized probability functions
where the use is unambiguous.

> Does this mean that, unless I have a good reason, I should prefer
> vector[n] x to real x[n]?

There's a chapter in the manual on when to prefer one over the other.
There's really no reason to have 1-dimensional real arrays other than
for generality of notation and similarity to BUGS. There are substantive
differences between 2D arrays and matrices. There are no integer
vectors or matrices, and only arrays can have 3+ dimensions.

- Bob

[Ross Boylan]

I had another thought about how to handle the situation below. It
didn't occur to me mostly because it's probably more work to implement.

If stan counted the errors (oops: I guess I should say "exceptional
conditions") internally and only reported a problem when they seemed
"high" that would avoid the useless noise and questions, but still alert
people when it looked as if there might be a real problem.

Ross

[Bob Carpenter]

That's been suggested several times, too. Or to just report
after a fixed percentage of warmup iterations, where most of them
occur.

It is a lot more trouble to implement and seems to me rather ad hoc
for this particular message. I'd rather just control our logging
levels.

Not clear what we'll actually do. It probably depends most on who
on the dev team volunteers to implement it and how much time they have.

- Bob

[Ross Boylan]

On 1/31/2014 3:33 PM, Andrew Gelman wrote:

HI, my quick hack would be to take the logit of the variable after transforming to a 0.5 to 4.5 scale, that is, work with logit((y - 0.5)/4).  This will make it continuous and will spread out the extremes a bit.  Of course it won't give correct predictions but otherwise it could do what you need.
If it works (or if it doesn't), please let me know!
A

The resulting predicted distribution, which I have finally examined, doesn't resemble the original very well.  Almost all the mass is near the extreme values.

I interpreted "won't give correct predictions" as meaning

1. the range will be more extreme
2. values will be predicted that don't exist in the original (original is always a multiple of 1/3)
3. exact shape, including clustering around integers, won't match.

But maybe you meant something stronger.  At any rate, I think I need to hit a little closer to the mark.  Of most relevance for my purposes, it seems likely that predicted effects will cause a bigger shift in the outcome--namely a flip from one extreme to the other--than a better model would show.

Here are the two distributions:

> quantile(simvals, probs=seq(0, 1, .05))  # simulated from the model
       0%        5%       10%       15%       20%       25%       30%       35% 
0.5000000 0.5002916 0.5022871 0.5089260 0.5259889 0.5634938 0.6390487 0.7803809 
      40%       45%       50%       55%       60%       65%       70%       75% 
1.0294206 1.4209338 1.9705903 2.6217905 3.2569433 3.7685059 4.1151658 4.3187676 
      80%       85%       90%       95%      100% 
4.4245009 4.4738727 4.4933017 4.4991461 4.5000000 
> quantile(tfat$FrndPartRace, probs=seq(0, 1, .05))  # original data
      0%       5%      10%      15%      20%      25%      30%      35% 
1.000000 1.000000 1.000000 1.000000 1.000000 1.333333 1.666667 1.666667 
     40%      45%      50%      55%      60%      65%      70%      75% 
1.666667 2.000000 2.000000 2.000000 2.000000 2.333333 2.333333 2.666667 
     80%      85%      90%      95%     100% 
2.666667 3.000000 3.000000 3.333333 4.000000 

Ross

[Ross Boylan]

I think I have good enough results using what I think of as a 1-parameter beta model, where the mean is the one parameter.
I took the original 1-4 scale and divided by 5, leaving the range from .2 through .8.  This was to avoid putting a lot of mass at the end points, since the beta doesn't like that (the manual seems  to indicate 0 and 1 are not legal values).  Then I fit a beta model parameterized by mean and "size" (courtesy of wikipedia) where the regression model applied only to the mean.  So there was a 2nd beta parameter, but it was just a single parameter rather than a fuller model.  Finally, when simulating I pull values outside of [1, 4] to the endpoints, which gives me a bit of a pile up at the end points, just like the original data.

Small syntax question: I didn't declare some variables in the model block (I just noticed it); is that a problem (see below)?

Here's the code:

data {
   int<lower=0> N; // number of respondents
   int<lower=0> J;  //number of covariates
   int<lower=0> NSD;  // number of classes of beta for sd
   int<lower=1,upper=NSD> SDClass[J];  // beta[j, k] has sigma_beta[SDClass[j]]
   int<lower=0> NCLUSTER;  //number of seeds
   int iCluster[N]; // cluster id
   matrix[N, J] x;  // covariates
   vector<lower=1,upper=4>[N] y;  // outcome
}
transformed data {
   vector<lower=0,upper=1>[N] z;
   z <- y/5.0;
}
parameters {
    vector[J] beta_raw; // standard covariate coefficients
    vector[NCLUSTER] alpha_raw; // cluster effects
    real<lower=0> nu; // 'size' parameter
    real<lower=0> sigma_alpha;
    real<lower=0> sigma_beta[NSD];
}
transformed parameters {
    vector[J] beta;
    vector[NCLUSTER] alpha;
    for(ir in 1:J) {
       beta[ir] <- beta_raw[ir]*sigma_beta[SDClass[ir]];
    }
    alpha <- alpha_raw*sigma_alpha;
}
model {
    vector[N] mu;
    vector[N] parm1;
    vector[N] parm2;
    nu ~ cauchy(0, 5); //maybe it should start above 0
    sigma_alpha ~ cauchy(0, 5);
    sigma_beta ~ cauchy(0, 5);
    beta_raw[1] ~ uniform(-5, 5);
    for (j in 2:J) {
      beta_raw[j] ~ normal(0, 1);
    }
    alpha_raw ~ normal(0, 1);
    for (n in 1:N) {
        mu[n] <- inv_logit((x[n] * beta)+alpha[iCluster[n]]);
    }
    parm1 <- mu*nu;
    parm2 <- (1-mu)*nu;
    z ~ beta(parm1, parm2);
}

parm1 and parm2 don't have declared data type, which seems to contradict the first sentence of section 21 of the manual.  Is that a problem?

Ross

[Bob Carpenter]

On 2/7/14, 7:48 PM, Ross Boylan wrote:
> I think I have good enough results using what I think of as a 1-parameter beta model, where the mean is the one parameter.
> I took the original 1-4 scale and divided by 5, leaving the range from .2 through .8. This was to avoid putting a lot
> of mass at the end points, since the beta doesn't like that (the manual seems to indicate 0 and 1 are not legal
> values). Then I fit a beta model parameterized by mean and "size" (courtesy of wikipedia)

Andrew discusses this parameterization in chapter 5 of BDA and we also have
an example in the reparameterization section of the Stan manual.

> where the regression model
> applied only to the mean.

Though we didn't use a regression for the mean on the logit scale, which
is a nice trick.

I'm not sure I understand. You declare them as local
variables in the model:

> vector[N] parm1;
> vector[N] parm2;

- Bob

[Ross Boylan]

On 2/7/2014 10:53 AM, Bob Carpenter wrote:
>> parm1 and parm2 don't have declared data type, which seems to
>> contradict the first sentence of section 21 of the
>> manual. Is that a problem?
>
> I'm not sure I understand. You declare them as local
> variables in the model:
>
> > vector[N] parm1;
> > vector[N] parm2;

You're right: no mystery.
Ross

[Ross Boylan]

On 2/1/2014 1:48 PM, Ross Boylan wrote:
> > system.time(sflist <- mclapply(1:6, function(i) stan(fit=rb, data=stanDat, iter=4000, chains=1,
> + control=list(adapt_delta=0.99),
> + seed=seed, chain_id=i, refresh=-1,
> + warmup=500, thin=4),
> + mc.cores=6)
> + )
>

> Warning messages:
> 1: In selectChildren(ac, 1) : error 'Interrupted system call' in select
> 2: In selectChildren(ac, 1) : error 'Interrupted system call' in select
> 3: In selectChildren(ac, 1) : error 'Interrupted system call' in select
> 4: In selectChildren(ac, 1) : error 'Interrupted system call' in select
> 5: In selectChildren(ac, 1) : error 'Interrupted system call' in select
>
>

I've been getting more of these selectChildren warnings and looked into it.
Apparently they require all of
1) using R (specifically the parallel library mclapply)
2) using Rmpi
3) using an old version of openmpi.

See
http://stackoverflow.com/questions/14149361/rmpi-mclapply-in-selectchildrenac-1-error-interrupted-system-call-in-s
for details (and the link there to the openmpi list).

Apparently loading the Rmpi library resets SIGCHLD handling (or did with
older open mpi's) which produces the warnings.

They appear harmless. And, to repeat, r/stan is a bystander to this
infelicity, not a cause.

Ross Boylan