Wrapped cauchy density implementation

[benlambertStan]

Summary:

Add wrapped cauchy density to Stan math library.

Description:

Implement a wrapped cauchy distribution (a useful distribution in circular statistics) as described in ``An extended family of circular distributions related to wrapped Cauchy distributions via Brownian motion'', S. Kato and M.C. Jones, Bernoulli 19(1), 2013, 154–171. This distribution is characterised by a location parameter -pi()<=mu<=pi(), and a scale parameter 0<=rho<=1, with support from -pi()<=y<=pi().

Current Output:

This function doesn't exist in the Stan math library.

Expected Output:

Five functions: wrapped_cauchy_log, wrapped_cauchy_cdf, wrapped_cauchy_cdf_log, wrapped_cauchy_ccdf_log, wrapped_cauchy_rng.

Current Version:

v2.12.0

[bob-carpenter]

Ideally, there would also be a _lcdf and _lccdf function to complement the _lpdf and _rng.

We have to decide whether to restrict input to (-pi, pi) or to treat it the way we treated von Mises and just work mod pi and force users to pick normalization interval.

For definition of PDF, see: https://en.wikipedia.org/wiki/Wrapped_Cauchy_distribution

But it didn't have a CDF formula or an indication of how to generate random variates.

[benlambertStan]

Hi Bob,

Ok, I have now cloned the Stan/math repository and have added five
functions to my local copy: wrapped_cauchy_log, wrapped_cauchy_cdf,
wrapped_cauchy_cdf_log, wrapped_cauchy_ccdf_log, wrapped_cauchy_rng, which
appear to work ok.

I know that the latest Stan functions use _lpdf, but I was unsure whether
this is done at the 'function name' level, or handled later on? This is
because I found the prob function declarations in
math-develop\stan\math\prim\scal\prob
were suffixed _log.

I also had a few questions though, due to my inexperience with both Stan
development and git, I'm afraid:

• Should wrapped_cauchy_cdf return the partial derivatives of the log cdf
  or the cdf? At the moment, I have copied the approach taken in the existent
  _cdf densities, where the log cdf partial derivatives are multiplied by the
  cdf value. Is this correct?
• I have added tests for the four probability densities to: test\prob.
  These take the form of example parameter values and the corresponding value
  of the log prob. I couldn't see any direct tests of the partial derivatives
  however. Where are these implemented?
• I have added a test of the _rng capability of the wrapped_cauchy_rng
  function to: test\unit\math\prim\scal\prob. Do I need to do anything else
  here?
• Do I need to add anything else to my clone of stan/math? For example, do
  I need to add a section in the manual somewhere?
• Once I am happy with my functions and tests I understand that I need to
  make sure that the following tests work:

make manual
make doxygen
make test-headers
./runTests.py src/test/unit
./runTests.py src/test/integration
./runTests.py src/test/prob

If these work ok (I am still trying to figure out how to do these on a
Windows machine) am I ok to submit a pull request?

Best,

Ben

On Wed, Sep 21, 2016 at 9:05 PM, Bob Carpenter notifications@github.com
wrote:

Ideally, there would also be a _lcdf and _lccdf function to complement
the _lpdf and _rng.

We have to decide whether to restrict input to (-pi, pi) or to treat it
the way we treated von Mises and just work mod pi and force users to pick
normalization interval.

For definition of PDF, see: https://en.wikipedia.org/wiki/
Wrapped_Cauchy_distribution

But it didn't have a CDF formula or an indication of how to generate
random variates.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
#396 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AMbAy32caTfLrKXdjC3abNzgKBH7036nks5qsY33gaJpZM4KCaf0
.

[bob-carpenter]

Nice. See inline.

On Oct 2, 2016, at 6:46 AM, benlambertStan notifications@github.com wrote:

Hi Bob,

Ok, I have now cloned the Stan/math repository and have added five
functions to my local copy: wrapped_cauchy_log, wrapped_cauchy_cdf,
wrapped_cauchy_cdf_log, wrapped_cauchy_ccdf_log, wrapped_cauchy_rng, which
appear to work ok.

I know that the latest Stan functions use _lpdf, but I was unsure whether
this is done at the 'function name' level, or handled later on? This is
because I found the prob function declarations in
math-develop\stan\math\prim\scal\prob
were suffixed _log.

Yes, we haven't changed the math library yet, so keep using
_log for and _cdf_log and _ccdf_log for now. We'll change them all
at once and then change the code generator to match and it'll simplify
a lot of the internal naming hacks in the parser/generator.

I also had a few questions though, due to my inexperience with both Stan
development and git, I'm afraid:

• Should wrapped_cauchy_cdf return the partial derivatives of the log cdf
  or the cdf?

Always the derivative of the function at hand. Having said that, we don't
really need the _cdf going forward. We do everything on the log scale.

At the moment, I have copied the approach taken in the existent
_cdf densities, where the log cdf partial derivatives are multiplied by the
cdf value. Is this correct?

Yes, but only because it cancels out to give the you derivative of
cauchy_cdf.

d/dx log cauchy_cdf(...x...)

= 1 / cauchy_cdf() * d/dx cauchy_cdf(...x...)

so multiplying this by the cdf value gives you the derivative of the cdf.

• I have added tests for the four probability densities to: test\prob.
  These take the form of example parameter values and the corresponding value
  of the log prob. I couldn't see any direct tests of the partial derivatives
  however. Where are these implemented?

Everything's done with finite differences in the test framework
I believe. Daniel?

• I have added a test of the _rng capability of the wrapped_cauchy_rng
  function to: test\unit\math\prim\scal\prob. Do I need to do anything else
  here?

No, it should only accept primitive arguments.

• Do I need to add anything else to my clone of stan/math? For example, do
  I need to add a section in the manual somewhere?

No, that'll go into stan-dev/stan when it gets added to the
set of known functions for the parser and documented.

It will need to get added to the appropriate top-level headers for
inclusion.

• Once I am happy with my functions and tests I understand that I need to
  make sure that the following tests work:

make manual
make doxygen
make test-headers
./runTests.py src/test/unit
./runTests.py src/test/integration
./runTests.py src/test/prob

If these work ok (I am still trying to figure out how to do these on a
Windows machine) am I ok to submit a pull request?

Be sure to do this:

make cpplint

You can skip the prob and integration tests and let those
run on the server. The prob tests take a huge amount of
disk space and time. Daniel can probably explain how to just
run one---I always forget.

• Bob

[spinkney]

@benlambertStan I could not find the wrapped cauchy code in your repo. Do you still have it?