isolating parameters for Gibbs sampling

[Bob Carpenter]

Warning: the punchline's in the very last paragraph.

On 2/27/13 10:35 PM, Andrew Gelman wrote:
> I was thinking of doing the hack where we Stan would alter the parameters-to-be-sliced
> one at a time and examine which lines of the model change, thus it would build a
> "hash table" (if I'm using that term correctly) telling it which lines of code
> to compute for each part of the slice updating.

By "parameters to be sliced", Andrew's talking about
conditional sampling for unbounded discrete parameters (I now
am confused about whether to call this "Gibbs" given that
we don't have conjugacy relations so we don't sample the
conditional exactly).

This would be easy if for parameter theta[k], we could efficiently
evaluate the conditional

log p(theta[k] | theta[1...k-1], theta[k+1,...K], y)

up to a constant that doesn't depend on theta[k].

It's the same issue for exact multinomial sampling
for bound parameters or if we wanted to do Gibbs
sampling of continuous variables for some reason.

It's easy to compute the above because it's proportional to
the joint posterior

log p(theta | y )

which we already have up to a constant. But that's very
very inefficient.

I understand what Andrew wants to do in principle. (Yes,
it would be a hash table that maps from parameters to
sequences of lines of code. Abstractly, a hash table
implements a function from one type of object to another.
A list in R implements such a map from names to values.)

But it's a whole lot harder to literally do what Andrew's
asking than it sounds. If Stan were interpreted line-by-line
it would be more feasible.

For instance, consider this simple coding for regression:

for (n in 1:N) {
real mu;
mu <- alpha + beta * x[n];
y[n] ~ normal(mu,sigma);
}

If alpha changes, we need to evaluate this whole thing.
Ditto beta. But they only show up on the line with the
assignment to the local variable mu.

One angle we could use is our current propto calculations
which is what lets us compute the log joint efficiently up
to a constant. If we could somehow make theta[k] the only
"variable" from the perspective of our trait metaprograms,
we'd pretty much get the behavior we wanted for free.

- Bob

[Bob Carpenter]

Peter:

I'd start with the univariate continuous distributions,
which live in:

src/stan/prob/distributions/univariate/continuous

This work should be done on the main branch unless you'd
be more comfortable working on your own branch. The reason
I suggest the main branch is that unlike the fvar stuff,
these being incomplete won't break the model compilation.

If you do want to keep working in fvar, you should merge
the master branch into fvar first so that you're reasonably
up to date with stan/prob.

There's currently a single example implementation in:

src/stan/prob/distributions/univariate/continuous/normal.hpp
------------------------------------------------------------
template <class RNG>
inline double
normal_rng(double mu,
double sigma,
RNG& rng) {
using boost::variate_generator;
using boost::normal_distribution;
using stan::math::value_of;
variate_generator<RNG&, normal_distribution<> >
rng_unit_norm(rng, normal_distribution<>());
return mu + sigma * rng_unit_norm();
}

where RNG is the templated out type of the basic random number
generator. The use of a reference type ("RNG&") in the variate_generator
is critical in order to advance the shared random number generator we're
going to use.

This implementation samples from a Normal(0,1) distribution then
scales by multiplying by sigma (the scale parameter) and relocates
by adding mu (the location parameter).

Boost supplies many of the basic distributions and the Wikipedia
has pretty good descriptions of how to sample from the various
distribtuions.

I'm not sure how best to test these. The current test is rather
weak, to say the least:

src/test/prob/distributions/univariate/continuous/normal_tes.cpp
----------------------------------------------------------------
TEST(ProbDistributionsNormal, random) {
boost::random::mt19937 rng;
EXPECT_NO_THROW(stan::prob::normal_rng(2.0,1.0,rng));
}

Anyone have any suggestions? We could draw thousands of samples
and then make sure the sample-based estimates are within the
sampling variance at a very high threshold (so we don't create
a lot of false positive tests). We can also do things like
test that the output is finite and appropriately signed (e.g.,
gamma variates should be positive, uniform variates shold be
constrained between their endpoints, etc.)

- Bob

[Ben Goodrich]

On Thursday, February 28, 2013 5:34:22 PM UTC-5, Bob Carpenter wrote:

I'm not sure how best to test these.  The current test is rather
weak, to say the least:

src/test/prob/distributions/univariate/continuous/normal_tes.cpp
----------------------------------------------------------------
TEST(ProbDistributionsNormal, random) {
   boost::random::mt19937 rng;
   EXPECT_NO_THROW(stan::prob::normal_rng(2.0,1.0,rng));
}

Anyone have any suggestions?  

We could do a test of the null hypothesis that the pseudo-random values (with a given seed) came from the supposed distribution. A chi-square test

http://www.itl.nist.gov/div898/handbook/eda/section3/eda35f.htm

is probably the easiest to implement and using a more powerful test would just increase the chance of a false positive.

Ben

[Peter Li]

How would I create the rng templates for distributions such as double exponential or inv_chi_square, which are not included in boost/random?


Thanks,
Peter

On Mar 1, 2013, at 3:01 AM, Bob Carpenter wrote:

> [I'm moving this discussion back to stan-dev. I'm not
> even the best person to help with these issues, but I
> can answer this one.]
>
> On 2/28/13 7:03 PM, Peter Li wrote:
>> Or should it be
>>
>> return exp(mu + 0.5 * sigma * sigma) + (exp(sigma * sigma) - 1) * exp(2 * mu + sigma * sigma) * rng_unit_lognorm();
>>
>> On Feb 28, 2013, at 6:43 PM, Peter Li wrote:
>>
>>> Alright. So lognormal would look like the following?

>>>
>>> template <class RNG>
>>> inline double

>>> lognormal_rng(double mu,

>>> double sigma,
>>> RNG& rng) {
>>> using boost::variate_generator;

>>> using boost::lognormal_distribution;
>>> using stan::math::value_of;
>>> using std::exp;
>>> variate_generator<RNG&, lognormal_distribution<> >
>>> rng_unit_lognorm(rng, lognormal_distribution<>());
>>> return exp(mu + sigma * rng_unit_lognorm());
>>> }
>
> Not quite. Check out the Wikipedia page:
>
> http://en.wikipedia.org/wiki/Log-normal_distribution
>
> section "Generating log-normally distributed random variates".
>
> Your scaled and relocated exponential transform would be right
> for an underlying normal distribution (i.e., replace
> lognormal_distribution<> with normal_distribution<>).
>
> You can probably feed the parameters into the lognormal_distribution
> in Boost and not do the scaling yourself. It might make sense to
> do that in the regular normal distribution, too.
>
> http://www.boost.org/doc/libs/1_53_0/doc/html/boost/random/lognormal_distribution.html
>
> It looks like you can just do the following, which is
> probably a better idea in general than trying to code our
> own versions:

>
> template <class RNG>
> inline double

> lognormal_rng(double mu,

> double sigma,
> RNG& rng) {
> using boost::variate_generator;

> using boost::lognormal_distribution;
> variate_generator<RNG&, lognormal_distribution<> >
> lognormal_rng(rng, lognormal_distribution<>(mu,sigma));
> return lognormal_rng();
> }
>
> This actually cleans up a bunch of the includes. The
> value_of shouldn't be included in the normal case, either,
> if you could clean that up.
>
>
> - Bob

[Jiqiang Guo]

There is a formula for double exponential (laplace) at 

http://en.wikipedia.org/wiki/Laplace_distribution

for generating sample from uniform(-.5, .5) first and then doing a transformation. 

For inv_chi_square, we can generate chi_square and take the reciprocal.  Or inv_chi_squaure is a special case of inverse gamma. 

--

Jiqiang 

[Bob Carpenter]

You need to figure out how to sample from them.
Usually, you can just transform other distributions.

For instance, the double exponential is aka "Laplace
distribution", and a method of sampling from it is
described here:

http://en.wikipedia.org/wiki/Laplace_distribution

in the section "Generating random variables according to
the Laplace distribution"

The inverse chi square is an inverse transform. I
don't see a Wikipedia how to generate example, but I'm
sure you can find it somewhere. Or you can work it
out from scratch. I'd suggest reading the transformations
section of the manual (near the back), that explains the
basics.

- Bob

[Peter Li]

Thanks. So how would I implement the chi-squared goodness of fit test that Ben suggested?


Thanks,
Peter

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

[Bob Carpenter]

On 3/1/13 2:57 PM, Peter Li wrote:
> Thanks. So how would I implement the chi-squared goodness of fit test that Ben suggested?

How familiar are you with basic hypothesis testing?
It's not that complex mathematically, but can get confusing
conceptually with all the philosophy, and mathematical
details and side conditions. If looking at the refs below
isn't enough, we can implement one next week.

I'll also build in the code generator so we can plug these
in to the release after next for Stan.

Here's some background, and feel free to bring up questions.

http://en.wikipedia.org/wiki/Goodness_of_fit

Specifically, Ben's suggesting a chi-squared test

http://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test

For continuous distributions (and for unbounded
or large discrete distributions), you have to bin the
output using the cumulative distribution function.

Here's a nice walkthrough of an example for the normal.

http://courses.wcupa.edu/rbove/Berenson/10th%20ed%20CD-ROM%20topics/section12_5.pdf

The trick is then adjusting the chi-squared test for the
degrees of freedom, which will be just 1 for us, because we'll
know the parameters from which we were supposed to have
generated the samples (the other degrees of freedom come from
having to estimate those from the data).

The code will need inverse CDFs, for which you can use the
quantile functions from Boost, for example, the normal:

http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/normal_dist.html

I don't think this will work, by the way, for the Cauchy,
which doesn't have finite means or variances.

- Bob

[Peter Li]

This is what I've done:

TEST(ProbDistributionsNormal, isnormal) {
boost::random::mt19937 rng;

boost::math::normal_distribution<>dist (10.0, 2.0);

double loc[4];
for(int i = 1; i < 5; i++)
loc[i - 1] = quantile(dist, 0.2 * i);

int count = 0;
double a = 0;
int bin [5] = {0, 0, 0, 0, 0};

while(count < 1000)
{
a = stan::prob::normal_rng(10.0,2.0,rng);
if(a > loc[3])
bin[4]++;
else if(a < loc[3] && a > loc[2])
bin[3]++;
else if(a < loc[2] && a > loc[1])
bin[2]++;
else if(a < loc[1] && a > loc[0])
bin[1]++;
else
bin[0]++;
count++;
}

double chi = 0;
double expect [5] = {200, 200, 200, 200, 200};

for(int j = 0; j < 5; j++)
chi += ((bin[j] - expect[j]) * (bin[j] - expect[j]) / expect[j]);

EXPECT_TRUE(chi < 9.49);
}

I believe the degrees of freedom should be n-1 where n is the number of "bins". Is there a more efficient way to do this? Also, how do I provide a formula for the chi-square value, instead of just 9.49.


Thanks,
Peter

[Bob Carpenter]

This looks great.

More comments inline.

On 3/1/13 8:16 PM, Peter Li wrote:
> This is what I've done:
>
> TEST(ProbDistributionsNormal, isnormal) {
> boost::random::mt19937 rng;
>
> boost::math::normal_distribution<>dist (10.0, 2.0);
>
> double loc[4];
> for(int i = 1; i < 5; i++)
> loc[i - 1] = quantile(dist, 0.2 * i);
>
> int count = 0;
> double a = 0;

You should declare variable "a" closer to where
you use it (in the loop; see below).

> int bin [5] = {0, 0, 0, 0, 0};
>
> while(count < 1000)
> {
> a = stan::prob::normal_rng(10.0,2.0,rng);
> if(a > loc[3])
> bin[4]++;
> else if(a < loc[3] && a > loc[2])
> bin[3]++;
> else if(a < loc[2] && a > loc[1])
> bin[2]++;
> else if(a < loc[1] && a > loc[0])
> bin[1]++;
> else
> bin[0]++;
> count++;
> }

You can do this loop a little more compactly and generally
as follows.

while (count < 1000) {
double a = stan::prob::normal_rng(10.0,2.0,rng);
int i = 0;
while (i < 4 && a > loc[i]) ++i;
++bin[i];
}

Then you can pull out the total count (N=1000) and
total number of bins (K=5) and replace the hardcoded
1000 with N and the hardcoded 4 with (K-1).

>
> double chi = 0;
> double expect [5] = {200, 200, 200, 200, 200};
>
> for(int j = 0; j < 5; j++)
> chi += ((bin[j] - expect[j]) * (bin[j] - expect[j]) / expect[j]);
>
> EXPECT_TRUE(chi < 9.49);
> }
>
> I believe the degrees of freedom should be n-1 where n is the number of "bins". Is there a more efficient way to do this?

Right. You'd normally then add the number of
free parameters estimated in the goodness of
fit test, but here we know the parameters. So
it's just number of cells minus 1 (the minus 1 is
because they're constrained to add up, so once
you know K-1 of them, you know what the Kth is).

> Also, how do I provide a formula for the chi-square value, instead of just 9.49.

The value you want is the quantile (inverse CDF) for the
chi-square distribution for a quantile like 1 - 1e-6
where we won't get too many false positives in the tests.
To do this, you want to use the complementary quantile functions
because the regular quantiles are arithmetically unstable for values
close to 1.

See:

http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/nmp.html#math_toolkit.dist.dist_ref.nmp._anchor_id__math_dist_quantile_c___quantile_from_the_complement_of_the_probability_

Boy, things sure are hard to find in the boost doc!
And what's up with these paths?

It might be even better here to use a quantile value like
0.95 and then allow multiple tries so that you get a
false positive test rate per test of 1e-6 or so.

- Bob

[Peter Li]

Great. Thanks, but for some reason, when I try to run a test with something like:

make CC=clang++ O=0 -j2 test/prob/distributions/univariate/continuous/normal

nothing happens.

[Bob Carpenter]

On 3/1/13 10:07 PM, Peter Li wrote:

> Great. Thanks, but for some reason, when I try to run a test with something like:
>
> make CC=clang++ O=0 -j2 test/prob/distributions/univariate/continuous/normal
>
> nothing happens.

When you say "nothing happens", what does the shell print out?

Can you do a clean first and see what happens?

Are you on the master branch or in fvar?

Here's what I get working on the master branch.

~/stan(master)$ git pull
Already up-to-date.

~/stan(master)$ make clean-all
rm -f
rm -f
rm -f -r doc/api
cd src/docs/stan-reference; rm -f *.brf *.aux *.bbl *.blg *.log *.toc *.pdf *.out *.idx *.ilg *.ind *.cb *.cb2 *.upa
rm -f -r models src/stan/model/model_header.hpp.gch src/stan/model/model_header.hpp.pch src/stan/model/model_header.hpp.d
rm -f -r test bin doc
rm -f
rm -f src/test/gm/model_specs/compiled/*.cpp src/test/gm/model_specs/compiled/*.o
src/test/gm/model_specs/compiled/declarations src/test/gm/model_specs/compiled/dims
src/test/gm/model_specs/compiled/function_signatures1 src/test/gm/model_specs/compiled/function_signatures6
src/test/gm/model_specs/compiled/function_signatures7 src/test/gm/model_specs/compiled/function_signatures_bernoulli
src/test/gm/model_specs/compiled/function_signatures_beta
src/test/gm/model_specs/compiled/function_signatures_beta_binomial
src/test/gm/model_specs/compiled/function_signatures_binomial
src/test/gm/model_specs/compiled/function_signatures_categorical
src/test/gm/model_specs/compiled/function_signatures_cauchy
src/test/gm/model_specs/compiled/function_signatures_chi_square
src/test/gm/model_specs/compiled/function_signatures_dirichlet
src/test/gm/model_specs/compiled/function_signatures_double_exponential
src/test/gm/model_specs/compiled/function_signatures_exponential
src/test/gm/model_specs/compiled/function_signatures_gamma
src/test/gm/model_specs/compiled/function_signatures_hypergeometric
src/test/gm/model_specs/compiled/function_signatures_inv_chi_square
src/test/gm/model_specs/compiled/function_signatures_inv_gamma
src/test/gm/model_specs/compiled/function_signatures_logistic
src/test/gm/model_specs/compiled/function_signatures_lognormal
src/test/gm/model_specs/compiled/function_signatures_multinomial
src/test/gm/model_specs/compiled/function_signatures_neg_binomial
src/test/gm/model_specs/compiled/function_signatures_normal
src/test/gm/model_specs/compiled/function_signatures_ordered_logistic
src/test/gm/model_specs/compiled/function_signatures_pareto src/test/gm/model_specs/compiled/function_signatures_poisson
src/test/gm/model_specs/compiled/function_signatures_scaled_inv_chi_square
src/test/gm/model_specs/compiled/function_signatures_student_t_0
src/test/gm/model_specs/compiled/function_signatures_student_t_1
src/test/gm/model_specs/compiled/function_signatures_student_t_2
src/test/gm/model_specs/compiled/function_signatures_student_t_3
src/test/gm/model_specs/compiled/function_signatures_uniform
src/test/gm/model_specs/compiled/function_signatures_weibull
src/test/gm/model_specs/compiled/illegal_generated_quantities src/test/gm/model_specs/compiled/illegal_transformed_data
src/test/gm/model_specs/compiled/illegal_transformed_parameters src/test/gm/model_specs/compiled/mat_assign
src/test/gm/model_specs/compiled/print_chars
cd src/test/agrad/distributions/univariate/continuous; rm -f *_generated_test.cpp
cd src/test/agrad/distributions/univariate/discrete; rm -f *_generated_test.cpp
cd src/test/agrad/distributions/multivariate/continuous; rm -f *_generated_test.cpp
cd src/test/agrad/distributions/multivariate/discrete; rm -f *_generated_test.cpp

~/stan(master)$ make CC=clang++ O=0 -j4 test/prob/distributions/univariate/continuous/normal
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -Wc++11-extensions -I lib/gtest_1.6.0/include -I lib/gtest_1.6.0
lib/gtest_1.6.0/src/gtest-all.cc -o test/gtest.o
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/agrad/agrad.o src/stan/agrad/agrad.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/math/matrix.o src/stan/math/matrix.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/agrad/matrix.o src/stan/agrad/matrix.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/expression_grammar_inst.o
src/stan/gm/grammars/expression_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/program_grammar_inst.o
src/stan/gm/grammars/program_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/statement_2_grammar_inst.o
src/stan/gm/grammars/statement_2_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/statement_grammar_inst.o
src/stan/gm/grammars/statement_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/term_grammar_inst.o
src/stan/gm/grammars/term_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/var_decls_grammar_inst.o
src/stan/gm/grammars/var_decls_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/grammars/whitespace_grammar_inst.o
src/stan/gm/grammars/whitespace_grammar_inst.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -o bin/stan/gm/ast_def.o src/stan/gm/ast_def.cpp
ar rv test/libgtest.a test/gtest.o
ar: creating archive test/libgtest.a
a - test/gtest.o
ar -rs bin/libstan.a bin/stan/agrad/agrad.o bin/stan/math/matrix.o bin/stan/agrad/matrix.o
ar: creating archive bin/libstan.a
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -c -O0 -Wc++11-extensions -I lib/gtest_1.6.0/include -I lib/gtest_1.6.0
src/test/prob/distributions/univariate/continuous/normal_test.cpp -o test/prob/distributions/univariate/continuous/normal.o
ar -rs bin/libstanc.a bin/stan/gm/grammars/expression_grammar_inst.o bin/stan/gm/grammars/program_grammar_inst.o
bin/stan/gm/grammars/statement_2_grammar_inst.o bin/stan/gm/grammars/statement_grammar_inst.o
bin/stan/gm/grammars/term_grammar_inst.o bin/stan/gm/grammars/var_decls_grammar_inst.o
bin/stan/gm/grammars/whitespace_grammar_inst.o bin/stan/gm/ast_def.o
ar: creating archive bin/libstanc.a
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE
-Wno-unused-function -Wno-tautological-compare -O0 lib/gtest_1.6.0/src/gtest_main.cc
test/prob/distributions/univariate/continuous/normal.o -Wc++11-extensions -I lib/gtest_1.6.0/include -I lib/gtest_1.6.0
-o test/prob/distributions/univariate/continuous/normal test/libgtest.a -Lbin -lstan -Lbin -lstanc
test/prob/distributions/univariate/continuous/normal
--gtest_output="xml:test/prob/distributions/univariate/continuous/normal.xml"
Running main() from gtest_main.cc
[==========] Running 2 tests from 1 test case.
[----------] Global test environment set-up.
[----------] 2 tests from ProbDistributionsNormal
[ RUN ] ProbDistributionsNormal.random
[ OK ] ProbDistributionsNormal.random (0 ms)
[ RUN ] ProbDistributionsNormal.isnormal
[ OK ] ProbDistributionsNormal.isnormal (2 ms)
[----------] 2 tests from ProbDistributionsNormal (2 ms total)

[----------] Global test environment tear-down
[==========] 2 tests from 1 test case ran. (2 ms total)
[ PASSED ] 2 tests.
~/stan(master)$


- Bob

P.S. Daniel gave me some fancy shell stuff for my .profile
that prints out which branch I'm on if I'm in a git directory.
If you're interested, I can forward.

[Peter Li]

Problem solved. Thanks.

Peter

[Peter Li]

On Mar 1, 2013, at 9:44 PM, Bob Carpenter wrote:

The value you want is the quantile (inverse CDF) for the
chi-square distribution for a quantile like 1 - 1e-6
where we won't get too many false positives in the tests.
To do this, you want to use the complementary quantile functions
because the regular quantiles are arithmetically unstable for values
close to 1.

See:

http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/nmp.html#math_toolkit.dist.dist_ref.nmp._anchor_id__math_dist_quantile_c___quantile_from_the_complement_of_the_probability_

Boy, things sure are hard to find in the boost doc!
And what's up with these paths?

It might be even better here to use a quantile value like
0.95 and then allow multiple tries so that you get a
false positive test rate per test of 1e-6 or so.

Would something like this work:

quantile(complement(mydist, 1e-6)))

Thanks,

Peter

[Bob Carpenter]

I haven't tried it, but that's what it looks like from
the doc. mydist should be the chi-square distribution
with the appropriate degrees of freedom.

To make sure it's doing what you expect, you can test
vs. R (qchisq function; the lower.tail arg
is the ordinary vs. complement flag). You could also
test vs. the ordinary quantile, where you'd expect the
quantile for 0.9 to be equal to the complementary quantile
for 0.1.

- Bob

[Daniel Lee]

Just a heads up, one of the random number generator tests have been failing.

ProbDistributionsBeta.doesFit

http://d1m1s1b1.stat.columbia.edu:8080/job/Stan%20Unit%20Tests/526/testReport/

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Peter Li]

Thanks. I didn't know that that one got pushed. I'll fix it as soon as I can.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Andrew Gelman]

See comment below.

On Feb 28, 2013, at 2:08 PM, Bob Carpenter wrote:

> I understand what Andrew wants to do in principle. (Yes,
> it would be a hash table that maps from parameters to
> sequences of lines of code. Abstractly, a hash table
> implements a function from one type of object to another.
> A list in R implements such a map from names to values.)
>
> But it's a whole lot harder to literally do what Andrew's
> asking than it sounds. If Stan were interpreted line-by-line
> it would be more feasible.
>
> For instance, consider this simple coding for regression:
>
> for (n in 1:N) {
> real mu;
> mu <- alpha + beta * x[n];
> y[n] ~ normal(mu,sigma);
> }
>
> If alpha changes, we need to evaluate this whole thing.
> Ditto beta. But they only show up on the line with the
> assignment to the local variable mu.

Bob: The above is not quite what I had in mind. My proposal is not to _parse_ the code for alpha but rather to change alpha and then _evaluate_ the code and see which lines give different results. Only the lines that give different results would need to be recomputed during the slice sampling for alpha. Of course, this will fail under certain settings (such as a function that's equal to alpha for alpha>0 and equal to 0 for alpha<0) but I think it could work pretty well for most models. Also, we could re-check the eval every 10th or 100th iteration or so, just to check that the "active" lines of code haven't changed.

[Bob Carpenter]

Maybe we can discuss this at the next meeting. I don't
see how it (Andrew's suggestion below) can work, but he
keeps suggesting it, so I apparently haven't done a very good
job of explaining why I don't think it'll work.

Even if Andrew's suggestion could be made to work, it'll
be very hard to code, because we'd need to literally instrument
every line of generated code somehow, index it to store values,
etc. We really want to look at expressions, not lines. And
the lines and expressions have contexts.

Having said all that, I don't like the idea of using something
that's not guaranteed to work. We have no way to debug it if
it goes wrong. Maybe a limited language subset could work?

My point below is that even if we find the lines that change,
they depend on other lines. Andrew is perhaps imagining that
a loop is implemented like in BUGS or JAGS by blowing out the graph,
so that:

for (n in 1:2)
y[n] ~ normal(0,1);

was just

y[1] ~ normal(0,1);
y[2] ~ normal(0,1);

But that's just not how it works. So we'd have to index the lines
that changed by the loop indices. These can depend on arbitrary
calculations elsewhere, so it's hard to see how to do this.

- Bob

[Andrew Gelman]

Sure, let's discuss at the meeting. I did not think that Stan blows out the graph; I understand that Stan lines are executable and that vectorziation is done at once. But I still don't see why we can't do things line by line. Also I'm not so worried about it not always working, as I think we can have it do occasional checks as often as we want it to.

Right now I'd say Pr(Bob is right and my idea is hopeless) = .8, and since I'm probably overconfident, I guess that should be .9. But that last .1 makes it worth a discussion at the meeting, I think!
A