Mixture model challenge

[Bill Harris]

I'm excited to be getting involved with stan finally. I've done a few
simple models, mostly duplicating what's in the documentation, and now
I'm trying a mixture model and having difficulty.

This is essentially the model from the manual with a bit of fake data:

------------------------------------------------------------------------

## Generate data
require(rstan)
lambda <- 0.6
nsamp <- 1000
n1 <- floor(lambda * nsamp)
n2 <- nsamp - n1
mu1 <- 900
mu2 <- 1100
sigma1 <- 20
sigma2 <- 40
data <- c(rnorm(n1,mu1,sigma1), rnorm(n2,mu2,sigma2))
sdata <- data/1000
ndata <- rescale(data)
mdata <- data.frame(
idx=1:nsamp,
data=data,
sdata=sdata,
ndata=ndata)

## Set up model
mixmodel <- '
data {
int<lower=1> K; // # components
int<lower=1> N; // # points
real<lower=0> y[N]; // data
}
parameters {
simplex[K] theta; // mixing proportions
real<lower=0> mu[K]; // location
real<lower=0,upper=3> sigma[K]; // scale
}
model {
real ps[K];
for (k in 1:K) {
mu[k] ~ normal(0,3)T[0,];
}
for (n in 1:N) {
for (k in 1:K) {
ps[k] <- log(theta[k])
+ normal_log(y[n],mu[k],sigma[k]);
}
lp__ <- lp__ + log_sum_exp(ps);
}
}
'
mixdata <- list(K = 2, N = 1000, y = mdata$ndata)
fit.mix <- stan(model_code = mixmodel, data = mixdata,
iter = 10000, chains = 4, verbose = FALSE)
------------------------------------------------------------------------

When I run the model, I get

,----
| > print(fit.mix)
| Inference for Stan model: mixmodel.
| 4 chains: each with iter=1000; warmup=500; thin=1; 1000 iterations saved.
|
| mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
| theta[1] 0.5 0.1 0.1 0.4 0.4 0.5 0.6 0.6 2 7.7
| theta[2] 0.5 0.1 0.1 0.4 0.4 0.5 0.6 0.6 2 7.7
| mu[1] 0.4 0.2 0.3 0.1 0.2 0.4 0.7 0.8 2 90.7
| mu[2] 0.4 0.2 0.3 0.1 0.1 0.4 0.7 0.8 2 77.2
| sigma[1] 0.1 0.0 0.1 0.1 0.1 0.1 0.1 0.2 2 27.3
| sigma[2] 0.1 0.0 0.0 0.1 0.1 0.1 0.1 0.1 2000 10.9
| lp__ 463.4 27.7 39.2 395.6 457.3 485.7 486.9 487.9 2 33.0
|
| Samples were drawn using NUTS2 at Sat Mar 16 15:23:44 2013.
| For each parameter, n_eff is a crude measure of effective sample size,
| and Rhat is the potential scale reduction factor on split chains (at
| convergence, Rhat=1).
`----

That seems pretty far off from the fake data. I did add a bit of a
prior on mu[k], FWIW. Yesterday I spent a fair amount of time trying to
fit real data and got no better.

Is there something obvious that might lead to better convergence?

Thanks,

Bill
--
Bill Harris
Facilitated Systems
http://makingsense.facilitatedsystems.com/

[Ben Goodrich]

On Saturday, March 16, 2013 6:27:00 PM UTC-4, Bill Harris wrote:

I'm excited to be getting involved with stan finally.  I've done a few
simple models, mostly duplicating what's in the documentation, and now
I'm trying a mixture model and having difficulty.

A colleague was having a similar problem with a mixture model this week. Can you verify whether there are multiple modes? There are some R functions to facilitate this. If there are multiple modes, the traceplot() function should have as many lines as you have chains that each flatten out but at different places. Also, pairs() of a stanfit-object with optional argument chain_id = c(i,j) will contrast the bivariate kernel densities for the i-th and j-th chains with the lower-diagonal and upper-diagonal of the plot respectively.

Ben

[Bill Harris]

Ben,

I should have checked traceplot. Yes, there are multiple modes, and
there were in my real data model yesterday, too. theta 1 and 2 both
settle into 2 modes around 400 iterations, mu 1 stablizes at 2 modes
around 450, mu 2 at around 350, and the like.

You should be able to duplicate that easily; the code I attached
previously is the entire code I loaded and ran (except for loading arm,
ggplot2, and scales, which I didn't need for the example).

My earlier attempts had less informative priors and much larger Rhat
estimates.

I think I misunderstood your other suggestion; I got

,----
| > pairs(fit.mix,chain_id = c(1,2))
| There were 50 or more warnings (use warnings() to see the first 50)
| > warnings()
| Warning messages:
| 1: In plot.window(...) : "chain_id" is not a graphical parameter
| 2: In plot.xy(xy, type, ...) : "chain_id" is not a graphical parameter
| .
| .
| .
`----

pairs(fit.mix) gave lots of plots with two clouds and, often, a seeming
outlier in another corner.

I presume we don't post graphics here, but I'll be glad to attach the
graphs if that helps and is allowed.

[Ben Goodrich]

On Saturday, March 16, 2013 8:02:31 PM UTC-4, Bill Harris wrote:

| > pairs(fit.mix,chain_id = c(1,2))
| There were 50 or more warnings (use warnings() to see the first 50)
| > warnings()
| Warning messages:
| 1: In plot.window(...) : "chain_id" is not a graphical parameter
| 2: In plot.xy(xy, type, ...) : "chain_id" is not a graphical parameter

Do you have rstan v1.2.0? Those warnings sound like an earlier version of rstan.
 

I presume we don't post graphics here, but I'll be glad to attach the
graphs if that helps and is allowed.

Graphs are fine. I know that Google Groups accepts .png or .pdf files.

Anyway, if there are multiple modes, then the Rhat and n_eff values that rstan produces are suspect. But if there are multiple modes, then you have a more challenging inferential problem anyway.

Ben
 

[Bill Harris]

Ben Goodrich <goodri...@gmail.com> writes:

> On Saturday, March 16, 2013 8:02:31 PM UTC-4, Bill Harris wrote:

> Do you have rstan v1.2.0? Those warnings sound like an earlier version of
> rstan.

I have v 1.1.1. I'll have to upgrade and get back to you.

> Graphs are fine. I know that Google Groups accepts .png or .pdf files.

Okay. Here's the traceplot.

> Anyway, if there are multiple modes, then the Rhat and n_eff values that
> rstan produces are suspect. But if there are multiple modes, then you have
> a more challenging inferential problem anyway.

I think the problem is inferential, not computational. I gave up trying
to model real data with a similar model Friday and went back to an EM
approach (mixtools), but I'd really like to do this with stan. I'm
somewhat new to BDA, and so I'm sure I'm doing some strange things (or
not doing some good things). If anyone has suggestions to improve the
model, that would be great.

[Bill Harris]

Bill Harris <bill_...@facilitatedsystems.com> writes:

> Ben Goodrich <goodri...@gmail.com> writes:
>> Anyway, if there are multiple modes, then the Rhat and n_eff values that
>> rstan produces are suspect. But if there are multiple modes, then you have
>> a more challenging inferential problem anyway.

I modified the model slightly and got better mixing on everything but
theta[k].

------------------------------------------------------------------------
## Generate data

require(YaleToolkit)
require(arm)
require(rstan)
require(ggplot2)
require(scales)

lambda <- 0.6
nsamp <- 1000
n1 <- floor(lambda * nsamp)
n2 <- nsamp - n1
mu1 <- 900
mu2 <- 1100
sigma1 <- 20
sigma2 <- 40

data <- c(rnorm(n1,mu1,sigma1), rnorm(n2,mu2,sigma2))
sdata <- data/1000
ndata <- rescale(data)
mdata <- data.frame(
idx=1:nsamp,
data=data,
sdata=sdata,
ndata=ndata)

rm(data)
rm(sdata)
rm(ndata)

## Set up model

mixmodel <- '
data {
int<lower=1> K; // # components
int<lower=1> N; // # points
real<lower=0> y[N]; // data
}
parameters {
simplex[K] theta; // mixing proportions
real<lower=0> mu[K]; // location
real<lower=0,upper=3> sigma[K]; // scale
}
model {
real ps[K];
for (k in 1:K) {
mu[k] ~ normal(0,3)T[0,];
}
for (n in 1:N) {
for (k in 1:K) {
ps[k] <- log(theta[k])
+ normal_log(y[n],mu[k],sigma[k]);
}
lp__ <- lp__ + log_sum_exp(ps);
}
}
'

mixdata <- list(K = 2, N = 1000, y = mdata$sdata)

fit.mix <- stan(model_code = mixmodel, data = mixdata,

iter = 1000, chains = 4, verbose = FALSE)

------------------------------------------------------------------------

,----
| > print(fit.mix)
| Inference for Stan model: mixmodel.
| 4 chains: each with iter=1000; warmup=500; thin=1; 1000 iterations saved.
|
| mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

| theta[1] 0.3 0.3 0.4 0.0 0.0 0.0 0.3 1.0 2 402.0
| theta[2] 0.7 0.3 0.4 0.0 0.7 1.0 1.0 1.0 2 402.0
| mu[1] 1.9 0.4 1.6 0.1 1.0 1.2 2.6 6.1 14 1.1
| mu[2] 1.3 0.5 1.1 0.4 1.0 1.0 1.0 4.9 6 1.2
| sigma[1] 1.1 0.4 0.9 0.1 0.1 0.9 1.9 2.9 5 1.3
| sigma[2] 0.4 0.4 0.7 0.1 0.1 0.1 0.1 2.5 3 1.8
| lp__ 842.3 0.1 1.6 838.4 841.5 842.7 843.6 844.6 445 1.0
`----

[Ben Goodrich]

On Sunday, March 17, 2013 1:56:14 AM UTC-4, Bill Harris wrote:

What I get looks like

,----
| > source("mixturemodels.R")
|
| TRANSLATING MODEL 'mixmodel' FROM Stan CODE TO C++ CODE NOW.
| Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
|   failed to parse Stan model 'mixmodel' with error message:
| EXPECTATION FAILURE LOCATION: file=input; line=16, column=5
|
|     theta[k] ~ dirichlet(alpha);
|     ^-- here
|
|
| DIAGNOSTIC(S) FROM PARSER:
| no matches for function name="dirichlet_log"
|     arg 0 type=real
|     arg 1 type=vector
| unknown distribution=dirichlet
| Parser expecting: <statement>
`----

but that looks a lot like the example on p. 79, as best as I can tell.  

What am I missing?

Just

theta ~ dirichlet(alpha);

no indexing of theta.

Ben
 

[Bill Harris]

Ben Goodrich <goodri...@gmail.com> writes:

> theta ~ dirichlet(alpha);
>
> no indexing of theta.

Thanks, Ben. Now I see: in the example in the manual, the indexing of
theta and phi occur over different indices than the dimensions of alpha
and beta.

The model compiles and runs now. I just have to figure out how to get
sufficient mixing on theta to attain convergence.

[Bill Harris]

Bill Harris <bill_...@facilitatedsystems.com> writes:

> The model compiles and runs now. I just have to figure out how to get
> sufficient mixing on theta to attain convergence.

Now I've got a model initialization or formulaton question, not a Stan
usage question. I've created the following mixture model and am trying
to fit the fake data that the script generates. This script should
replicate the entire calculation.

vector<lower=0>[K] alpha; // prior

real<lower=0> y[N]; // data
}
parameters {
simplex[K] theta; // mixing proportions
real<lower=0> mu[K]; // location
real<lower=0,upper=3> sigma[K]; // scale
}
model {
real ps[K];

theta ~ dirichlet(alpha);

for (k in 1:K)
mu[k] ~ normal(0,3)T[0,];
for (n in 1:N) {
for (k in 1:K) {
ps[k] <- log(theta[k])
+ normal_log(y[n],mu[k],sigma[k]);
}
lp__ <- lp__ + log_sum_exp(ps);
}
}
'

mixdata <- list(K = 2, N = 1000, alpha = 1*c(0.6,0.4), y = mdata$sdata)

fit.mix <- stan(model_code = mixmodel, data = mixdata,

iter = 400, chains = 4, verbose = FALSE)
## iter = 1000, chains = 4, verbose = FALSE)

------------------------------------------------------------------------

So far, I've been unsuccessful in getting the chains to mix, and theta
makes the biggest problem.

I've tried the advice in chapter 18 of BDA to make the parameters of
Dir() be proportional to the relative sizes of the thetas and their
magnitude be proportional to the strength of the prior distribution.
With the current formulation, the parameters match the chosen values of
theta in the fake data, and the coefficient can be used to scale them up
or down.

I've tried the coefficient equal to 0.1 (making the distribution
Dir(0.06, 0.04)), 1, 10, and 100. I've also set alpha = c(1,1) and
other combinations. In no case can I recall getting the chain for theta
to mix, and sometimes the others don't, either.

With the coefficient equal to 5 (alpha = c(3,2)), I get

,----
| > print(fit.mix)
| Inference for Stan model: mixmodel.
| 4 chains: each with iter=1000; warmup=500; thin=1; 1000 iterations saved.
|
| mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

| theta[1] 0.4 0.3 0.4 0.0 0.1 0.3 0.7 1.0 2 52.4
| theta[2] 0.6 0.3 0.4 0.0 0.3 0.7 0.9 1.0 2 52.4
| mu[1] 1.2 0.3 0.9 0.5 0.9 0.9 1.0 4.2 7 1.2
| mu[2] 1.3 0.4 1.0 0.5 1.0 1.0 1.1 4.7 6 1.2
| sigma[1] 0.3 0.0 0.7 0.0 0.0 0.0 0.1 2.6 2000 1.6
| sigma[2] 0.4 0.4 0.7 0.0 0.0 0.1 0.1 2.7 3 1.7
| lp__ 1064.2 201.9 285.7 828.9 834.3 947.0 1171.6 1535.5 2 191.9
`----

Here's a traceplot.

[Bob Carpenter]

You probably want to put a prior on the scales sigma.

I was about to point you to the doc, when I realized
that the model there doesn't have a prior on sigma.
I'll go add it, because it usually helps with convergence
(even if not strictly necessary for the prior to be
proper).

Do you mean you have a problem getting the chains
of a mixture model to mix with each other? That's
typical because mixture models are not identified. Do
you get the same answer but with reversed
indices? The traceplots of theta[1] and theta[2]
dont' have enough resolutions for me to tell.

Also, you can replace all of this:

> for (k in 1:K)
> mu[k] ~ normal(0,3)T[0,];

with just

mu ~ normal(0,3);

becuase the normal's vectorized and you don't need the
truncation if both location and scale are constants and
the variate mu is properly constrained.

Your data's well separated, so I don't imagine there's
going to be a problem finding the modes.

Why are you doing this:

> mu1 <- 900
> mu2 <- 1100
> sigma1 <- 20
> sigma2 <- 40
>
> data <- c(rnorm(n1,mu1,sigma1), rnorm(n2,mu2,sigma2))
> sdata <- data/1000
> ndata <- rescale(data)

rather than just taking mu = 0.9, mu2 = 1.1 and
sigma1 = 0.02 and sigma2 = 0.04?

You should be able to get away with a Dirichlet(2,2)
or even Dirichlet(1,1) distribution, but I'd stay away
from Dirichlet(0.06, 0.04) for a problem where you expect
both items to have some count. If you don't include a prior
at all, it's equivalent to Dirichlet(1,1).

- Bob

> Does anyone have a suggestion for getting a mixture model such as this
> to converge? I didn't think it would be particularly challenging.
>
> Thanks,
>
> Bill
>
>
>

[Bob Carpenter]

Maybe I should go to bed :-) You have an implicit
uniform prior over the declared support of sigma, which
is proper here because the support is bounded.

Initialization will be at the halfway point, which is
still high (1.5) for your scaled down parameters with
sigma = 0.02 and 0.04. Stan can usually figure this out,
but it will take extra time to adapt compared to
initializations that are closer to their true values.

- Bob

On 3/18/13 1:11 AM, Bill Harris wrote:

[Bob Carpenter]

OK, just as a sanity check, I went back and made
sure that the model in the manual fits properly.

Here's the model:

normal_mixture_k.stan
---------------------------
data {
int<lower=1> K;
int<lower=1> N;
real y[N];
}
parameters {
simplex[K] theta;
real mu[K];
real<lower=0,upper=10> sigma[K];
}
model {
real ps[K];
mu ~ normal(0,10);

for (n in 1:N) {
for (k in 1:K)
ps[k] <- log(theta[k])
+ normal_log(y[n],mu[k],sigma[k]);
lp__ <- lp__ + log_sum_exp(ps);
}
}

-------------------------

Here's the R code I used for simulation and to run
it, which is like yours, but without the rescaling
or data-frame wrangling:

------------------------------------------
K <- 3;

N1 <- 500;
N2 <- 300;
N3 <- 200;
N <- N1 + N2 + N3;
mu1 <- -2.0;
mu2 <- 0.5;
mu3 <- 3.5;
sigma1 <- 1.0;
sigma2 <- 0.5;
sigma3 <- 2;

y <- c(rnorm(N1,mu1,sigma1),rnorm(N2,mu2,sigma2),rnorm(N3,mu3,sigma3));

library("rstan");
fit <- stan(model="normal_mixture_k.stan", iter=500, chains=1, init=0);
---------------------------------


And here's the print for the fit in R:


--------------------------------
> print(fit)
Inference for Stan model: normal_mixture_k.
1 chains: each with iter=500; warmup=250; thin=1; 500 iterations saved.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

theta[1] 0.5 0.0 0.0 0.4 0.5 0.5 0.5 0.5 144 1
theta[2] 0.2 0.0 0.0 0.2 0.2 0.2 0.2 0.3 104 1
theta[3] 0.3 0.0 0.0 0.3 0.3 0.3 0.3 0.4 76 1
mu[1] -2.1 0.0 0.1 -2.2 -2.1 -2.1 -2.0 -1.9 109 1
mu[2] 3.4 0.0 0.4 2.6 3.2 3.4 3.6 4.0 90 1
mu[3] 0.4 0.0 0.1 0.3 0.4 0.4 0.5 0.5 81 1
sigma[1] 1.0 0.0 0.1 0.9 0.9 1.0 1.0 1.1 101 1
sigma[2] 2.1 0.0 0.2 1.8 2.0 2.1 2.3 2.6 88 1
sigma[3] 0.6 0.0 0.1 0.5 0.6 0.6 0.6 0.7 102 1
lp__ -2186.7 0.3 2.0 -2192.0 -2187.6 -2186.3 -2185.4 -2184.0 49 1

Samples were drawn using NUTS2 at Mon Mar 18 03:31:18 2013.

For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

---------------------------------


The fitted mu[1]/sigma[1] correspond to the simulated mu1/sigma1,
the fitted mu[2]/sigma[2] correspond to the simulated mu3/sigma3,
and the fitted mu[3]/sigma[3] correspond to simulated mu2/sigma2.
And the mixing proportions are estimated correctly.

Note that I only ran a single chain. If you run multiple chains,
you'll find the indexes scrambling. This model is not identifiable.
In fact, none of the parameters is identifiable. The only thing
that will converge across chains is lp__.

If you need to assess convergence, you need to somehow re-identify
the model. One way to do this is to sort.

You'd need to sort one of the parameters, like mu, then
apply the same sort order to the other parameters. Then you
should be able to measure the convergence of the sorted parameters.
You want to choose a parameter that will sort properly, too.
If all the scaleses are the same, sorting on sigma won't work,
and if all the mixture proportions are the same, then sorting
on theta won't work.

We really need to add something like JAGS's sort in order to do
this easily within Stan. It should be simple to add because
there are no gradients to compute for sort.

Ho long it takes to converge will depend on how good the
starting points are. Random starting points seem to converge
in somewhere between a few dozen and few hundred iterations
depending on how close they are, and init=0 initializations
seem to converge reliably in a few dozen iterations.

- Bob

[Bill Harris]

Thanks, Bob, for all the suggestions.

In this note, you asked (essentially) why I scaled the data the way I did.  My real data will be more like the 900 and 1100, but I thought that scaling it might make the model work more easily.  Rescale() seemed like the best way to go, but dividing by 1000 seemed similar while giving me numbers that sounded familiar.

I'm still quite new to Dirichlet distributions, so your practical advice helps.  The reminders about vectorization help, too.

More later.

Bill

[Bob Carpenter]

On 3/18/13 7:52 AM, Bill Harris wrote:
> Thanks, Bob, for all the suggestions.

> In this note, you asked (essentially) why I scaled the data the way I did. My real data will be more like the 900 and
> 1100, but I thought that scaling it might make the model work more easily. Rescale() seemed like the best way to go,
> but dividing by 1000 seemed similar while giving me numbers that sounded familiar.

Yes, that will work better. Stan can figure out
the scaling during warmup, but it can slow down
convergence if the initialization is poor (not near
where the true parameter values are).

Also, it often happens with setting "default" priors.

> I'm still quite new to Dirichlet distributions, so your practical advice helps. The reminders about vectorization help,
> too.

You might want to think about Beta(theta|a,b) first,
because Dirichlet is just the multivariate version.
It's a bit different because it's univariate, so you
only get the first component of the two, but you
can visualize the various distributions more easily.

The basic idea is that your prior should be fat (or
at least flat) where you expect the parameter values to
be, not the other way around, as you see in a Beta(0.01,0.01)
if you're expecting a value that's not near 0 or 1.

- Bob

[Bill Harris]

Bob,

I can't make your model run on my Windows box the way you sent it, and I'm not near a Linux box right now to test it there.

I modified your .R file to include the model definition a string called bobcmodel:

--------------------------------------------------------------------------------------------------------------------

N <- N1 + N2 + N3; 

mu1 <- -2.0; 

mu2 <- 0.5; 

mu3 <- 3.5; 

sigma1 <- 1.0; 

sigma2 <- 0.5; 

sigma3 <- 2; 

y <- c(rnorm(N1,mu1,sigma1),rnorm(N2,mu2,sigma2),rnorm(N3,mu3,sigma3)); 

bobcmodel <- '

data { 

   int<lower=1> K; 

   int<lower=1> N; 

   real y[N]; 

} 

parameters { 

   simplex[K] theta; 

   real mu[K]; 

   real<lower=0,upper=10> sigma[K]; 

} 

model { 

   real ps[K]; 

   mu ~ normal(0,10); 

   for (n in 1:N) { 

     for (k in 1:K) 

       ps[k] <- log(theta[k]) 

                + normal_log(y[n],mu[k],sigma[k]); 

     lp__ <- lp__ + log_sum_exp(ps); 

   } 

}

'

library("rstan"); 

## fit <- stan(model="normal_mixture_k.stan", iter=500, chains=1, init=0);

## fit <- stan(model_code = bobcmodel, iter=500, chains=1, init=0);

## fit <- stan(model_code = bobcmodel, data = list(K=K,N=N,y=y),iter=500, chains=1, init=0);

## fit <- stan(model = "normal_mixture_k.stan", data = list(K=K,N=N,y=y),iter=500, chains=1, init=0);

## fit <- stan(file="normal_mixture_k.stan", iter=500, chains=1, init=0);

fit <- stan(file = "normal_mixture_k.stan", data = list(K=K,N=N,y=y),iter=500, chains=1, init=0);

--------------------------------------------------------------------------------------------------------------------  

If I run it as you sent it, using the first stan() call, I get

> source("normal_mixture_k.R")

Error in stan(model = "normal_mixture_k.stan", iter = 500, chains = 1,  (from normal_mixture_k.R#40) : 

  argument 1 matches multiple formal arguments

If I run the second stan() call (the one that references the model code rather than the model file), I get

> source("normal_mixture_k.R")

TRANSLATING MODEL 'bobcmodel' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'bobcmodel' NOW.

cygwin warning:

  MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf

  Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf

  CYGWIN environment variable option "nodosfilewarning" turns off this warning.

  Consult the user's guide for more details about POSIX paths:

    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames

C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]

Error in eval(expr, envir, enclos) : 

  variable does not exist; processing stage=data initialization; variable name=K; base type=int

failed to create the sampler; sampling not done

> K

[1] 3

If I run it with the third stan() call, it works:

> source("normal_mixture_k.R")

TRANSLATING MODEL 'bobcmodel' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'bobcmodel' NOW.

cygwin warning:

  MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf

  Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf

  CYGWIN environment variable option "nodosfilewarning" turns off this warning.

  Consult the user's guide for more details about POSIX paths:

    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames

C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]

SAMPLING FOR MODEL 'bobcmodel' NOW (CHAIN 1).

Iteration:   1 / 500 [  0%]  (Adapting)

Iteration:  50 / 500 [ 10%]  (Adapting)

If I try the fourth stan() model, I get

> source("normal_mixture_k.R")

Error in stan(model = "normal_mixture_k.stan", data = list(K = K, N = N,  (from normal_mixture_k.R#46) : 

  argument 1 matches multiple formal arguments

The fifth fails, too:

> source("normal_mixture_k.R")

TRANSLATING MODEL 'normal_mixture_k' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'normal_mixture_k' NOW.

cygwin warning:

  MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf

  Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf

  CYGWIN environment variable option "nodosfilewarning" turns off this warning.

  Consult the user's guide for more details about POSIX paths:

    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames

C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]

Error in eval(expr, envir, enclos) : 

  variable does not exist; processing stage=data initialization; variable name=K; base type=int

failed to create the sampler; sampling not done

The sixth version works:

> source("normal_mixture_k.R")

TRANSLATING MODEL 'normal_mixture_k' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'normal_mixture_k' NOW.

cygwin warning:

  MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf

  Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf

  CYGWIN environment variable option "nodosfilewarning" turns off this warning.

  Consult the user's guide for more details about POSIX paths:

    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames

C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]

SAMPLING FOR MODEL 'normal_mixture_k' NOW (CHAIN 1).

Iteration:   1 / 500 [  0%]  (Adapting)

Iteration:  50 / 500 [ 10%]  (Adapting)

In other words, I seem to need to pass in the data as an explicit list, not as globals, as you seemed to be able to do.  Furthermore, the argument of the file name on my Windows box is called "file", not "model".  I don't know if those are differences between Windows and other OSes or not.  The error message (especially the "argument 1 matches multiple formal arguments") weren't clear to me for a while.

Thanks for this, though.  Now I've got something that might help me with my inference challenge.

Bill

[Bob Carpenter]

Responses inline.

On 3/18/13 12:33 PM, Bill Harris wrote:
> Bob,
>
> I can't make your model run on my Windows box the way you sent it, and I'm not near a Linux box right now to test it there.

There shouldn't be a difference between the
way RStan runs under Windows or under Linux other
than the paths you use to refer to files.

From the stan() doc in R:

file
A character string file name or a connection that R supports
containing the text of a model specification in the Stan modeling
language; a model may also be specified directly as a character
string using parameter model_code or through a previous fit using
parameter fit. When fit is specified, parameter file is ignored.

model_code
A character string either containing the model definition
or the name of a character string object in the workspace.
This parameter is used only if parameter file is not specified.
When fit is specified, the model compiled previously is used so
specifying model_code is ignored.

So its' file=model_file or model_code=model_string.

You need version 3 or 6, which correctly specify the
model and the data.

(Sorry if I said "model" at some point -- I keep doing
that in my own code and having to go back and look at the
doc.)

I think this might be easier if it were overloaded, where
we just say model=blah_blah_blah and RStan first sees
if that referes to a file, and if not interprets it
as a string.

> --------------------------------------------------------------------------------------------------------------------
>
> If I run it as you sent it, using the first stan() call, I get
>
> > source("normal_mixture_k.R")
> Error in stan(model = "normal_mixture_k.stan", iter = 500, chains = 1, (from normal_mixture_k.R#40) :
> argument 1 matches multiple formal arguments

Right -- "model" isn't a defined key.
Apparently R lets you use prefixes of argument
names, which is why this error comes up.

> If I run the second stan() call (the one that references the model code rather than the model file), I get
>
> > source("normal_mixture_k.R")
>
> TRANSLATING MODEL 'bobcmodel' FROM Stan CODE TO C++ CODE NOW.
> COMPILING THE C++ CODE FOR MODEL 'bobcmodel' NOW.
> cygwin warning:
> MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf
> Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf
> CYGWIN environment variable option "nodosfilewarning" turns off this warning.
> Consult the user's guide for more details about POSIX paths:
> http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
> C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but
> not used [-Wunused-function]
> Error in eval(expr, envir, enclos) :
> variable does not exist; processing stage=data initialization; variable name=K; base type=int
> failed to create the sampler; sampling not done

This also shouldn't run because you haven't specified the
data. You either need an array of variable names to read
from the environment, or a list of variable definitions.

> > K
> [1] 3
>
> If I run it with the third stan() call, it works:
>
> > source("normal_mixture_k.R")
>
> TRANSLATING MODEL 'bobcmodel' FROM Stan CODE TO C++ CODE NOW.
> COMPILING THE C++ CODE FOR MODEL 'bobcmodel' NOW.
> cygwin warning:
> MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf
> Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf
> CYGWIN environment variable option "nodosfilewarning" turns off this warning.
> Consult the user's guide for more details about POSIX paths:
> http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
> C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but
> not used [-Wunused-function]
> SAMPLING FOR MODEL 'bobcmodel' NOW (CHAIN 1).
> Iteration: 1 / 500 [ 0%] (Adapting)
> Iteration: 50 / 500 [ 10%] (Adapting)
>
> If I try the fourth stan() model, I get
>
> > source("normal_mixture_k.R")
> Error in stan(model = "normal_mixture_k.stan", data = list(K = K, N = N, (from normal_mixture_k.R#46) :
> argument 1 matches multiple formal arguments

Same problem as 1.

> The fifth fails, too:
>
> > source("normal_mixture_k.R")

> TRANSLATING MODEL 'normal_mixture_k' FROM Stan CODE TO C++ CODE NOW.
> COMPILING THE C++ CODE FOR MODEL 'normal_mixture_k' NOW.
> cygwin warning:
> MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf
> Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf
> CYGWIN environment variable option "nodosfilewarning" turns off this warning.
> Consult the user's guide for more details about POSIX paths:
> http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
> C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but
> not used [-Wunused-function]
> Error in eval(expr, envir, enclos) :
> variable does not exist; processing stage=data initialization; variable name=K; base type=int
> failed to create the sampler; sampling not done

You need the data here.

> The sixth version works:
>
> > source("normal_mixture_k.R")
>
> TRANSLATING MODEL 'normal_mixture_k' FROM Stan CODE TO C++ CODE NOW.
> COMPILING THE C++ CODE FOR MODEL 'normal_mixture_k' NOW.
> cygwin warning:
> MS-DOS style path detected: c:/R/R-215~1.3/etc/x64/Makeconf
> Preferred POSIX equivalent is: /cygdrive/c/R/R-215~1.3/etc/x64/Makeconf
> CYGWIN environment variable option "nodosfilewarning" turns off this warning.
> Consult the user's guide for more details about POSIX paths:
> http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
> C:/R/library/rstan/include//stansrc/stan/agrad/agrad.hpp:2190:17: warning: 'void stan::agrad::free_memory()' defined but
> not used [-Wunused-function]
> SAMPLING FOR MODEL 'normal_mixture_k' NOW (CHAIN 1).
> Iteration: 1 / 500 [ 0%] (Adapting)
> Iteration: 50 / 500 [ 10%] (Adapting)
>
> In other words, I seem to need to pass in the data as an explicit list, not as globals, as you seemed to be able to do.

If the variables are defined in the global environment, you
should be able to pass in data using the argument

data=c("K","N","y")

You can't just use c(K,N,y), because we need the names of
the variables as well as the data itself.

> Furthermore, the argument of the file name on my Windows box is called "file", not "model".

That's true everywhere -- sorry if I said "model" at some point.

> I don't know if those are
> differences between Windows and other OSes or not. The error message (especially the "argument 1 matches multiple
> formal arguments") weren't clear to me for a while.

No.

> Thanks for this, though. Now I've got something that might help me with my inference challenge.

You're welcome. Let us know if you have other issues.

- Bob

[Bill Harris]

Thanks, Bob.  It seems to be working.  

Bill