Question regarding the use of lookup table for an unavailable distribution

[Albert Lee]

Hi. 

I am trying to fit a model that has landau distribution. Since Landau distribution is not available in Stan, I decided to get away with this using the lookup table approach (I previously asked the related question in this post ).

I implemented just that using all sorts of hacks as you can see in the attached stan code (I can imagine the images of stan devs rolling their eyes looking at my code...).

With this code, I wanted to test if I can recover some parameters using the simulated data.

Simulation is as follows:

X ~ Landau() // generate 1000 data points from Landau distribution

Y = b0 + b1 * X // transform it to Y using b0 = 13, b1=5

The object is to recover b0 and b1 using only Y.

You can download the data file(landau_logpdf_lkup.stan.R) at 

https://drive.google.com/file/d/0B0ZDURib4NwAQmVEeThmcVZrYXM/view?usp=sharing

and I have attached the stan code(landau.stan)

The stan code compiles just fine, and the function "compute_landau_logpdf" I implemented returns the correct logpdf(see Expected_landau_logpdf_using_mathematica.png)

But using "optimize" method in cmdStan, I cannot recover b0 and b1.

Here's the output of running cmdstan

   ===> ./landau optimize data file=landau_logpdf_lkup.stan.R

 method = optimize

   optimize

     algorithm = lbfgs (Default)

       lbfgs

         init_alpha = 0.001 (Default)

         tol_obj = 9.9999999999999998e-13 (Default)

         tol_rel_obj = 10000 (Default)

         tol_grad = 1e-08 (Default)

         tol_rel_grad = 10000000 (Default)

         tol_param = 1e-08 (Default)

         history_size = 5 (Default)

     iter = 2000 (Default)

     save_iterations = 0 (Default)

 id = 0 (Default)

 data

   file = landau_logpdf_lkup.stan.R

 init = 2 (Default)

 random

   seed = 1001198607

 output

   file = output.csv (Default)

   diagnostic_file =  (Default)

   refresh = 100 (Default)

initial log joint probability = -4818.34

    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes

       1      -4818.34             0       316.886    0.005548       0.001       49

Optimization terminated with error:

  Line search failed to achieve a sufficient decrease, no more progress can be made

Expected logpdf for 3.23: -2.679

Expected logpdf for -5: -54.55

Expected logpdf for 10001: -18.42

And here's the output for output.csv

# stan_version_major = 2

# stan_version_minor = 9

# stan_version_patch = 0

# model = landau_model

# method = optimize

#   optimize

#     algorithm = lbfgs (Default)

#       lbfgs

#         init_alpha = 0.001 (Default)

#         tol_obj = 9.9999999999999998e-13 (Default)

#         tol_rel_obj = 10000 (Default)

#         tol_grad = 1e-08 (Default)

#         tol_rel_grad = 10000000 (Default)

#         tol_param = 1e-08 (Default)

#         history_size = 5 (Default)

#     iter = 2000 (Default)

#     save_iterations = 0 (Default)

# id = 0 (Default)

# data

#   file = landau_logpdf_lkup.stan.R

# init = 2 (Default)

# random

#   seed = 1001198607

# output

#   file = output.csv (Default)

#   diagnostic_file =  (Default)

#   refresh = 100 (Default)

lp__,b1,b0

-4818.34,13.4974,15.4925

My questions are

• Is there a fundamental flaw in using the lookup table approach?
• Or if the logpdf landau function is implemented correctly(which is big IF..but at least three test cases return correct results), what's the potential reason that optimization didn't work?
• Is there some crucial piece of knowledge I am missing?
• Is there any other suggestions?

Thank you!!

[Bob Carpenter]

It needs to yield the right derivatives, not just the right
values. What you're calling the AS-INTEGER HACK is fundamentally
broken.

There's no way to recover b0 and b1 using only Y --- you need
X, too!

The Wikipedia page

https://en.wikipedia.org/wiki/Landau_distribution

has a simple approximation, which is easily coded in Stan.

Also, optimization might not work in principle if whatever
you're optimizing doesn't have exactly one mode.

Also, we strongly discourage interval-constrained parameters
(your b1 and b0).

log(1/b1) is just -log(b1)

You might also want to independently test how you're generating
the data.

- Bob

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