Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/issue 2682 add 7 parameter ddm pdf #2822

Merged

Conversation

Franzi2114
Copy link
Collaborator

@Franzi2114 Franzi2114 commented Oct 7, 2022

Summary

We provide an implementation of the LPDF drift diffusion model with up to 7 parameters. There will be two new functions available in the stan::math namespace: wiener_full_lpdf and wiener_full_prec_lpdf. The only difference between both functions is the control of the precision in the computation of the partial derivatives. We did not deal with the CDF functions. With the new funtion it is possible to compute the 4-, 5-, 6-, 7- parameter drift diffusion models. Therefore, it is also an alternative to the wiener_lpdf by @kendalfoster.

The PR consists of the following new files:
Main file:
-stan/math/prim/prob/wiener_full_lpdf.hpp.
Helper files:
-stan/math/prim/fun/wiener5_lpdf.hpp: log density and partial derivatives for 5 parameter model
-stan/math/prim/fun/wiener7_lpdf.hpp: log density and partial derivatives for 7 parameter model
-stan/math/prim/functor/hcubature.hpp: new multidimensional integrator, adapted from Steven Johnson
-test/unit/math/prim/prob/wiener_full_test.cpp: own unit tests
-test/unit/math/prim/prob/wiener_full_prec_test.cpp: own unit tests for function with precision control

The following files are adapted in order to be able to implement a function with 8 parameters:
-stan/math/prim/functor/operands_and_partials.hpp
-stan/math/rev/functor/operands_and_partials.hpp

Discourse threads:
Testing for 7-parameter function with known partials -> Therefore, we only let scalar input be accepted and no vectorization.
Is there a multidimensional integrate-function? -> Therefore, we implemented hcubature.
How to document correctly?/How to add a new function’s documentation to the Stan Function Reference? -> Up to now, we have the doxygen documentation.
Adding a function with LARGE known partial derivatives -> We introduce 3 helper files.
HMC: Does precision of derivatives affect validity of results?

Issue:
Extension operands_and_partials to 8 edges #2698

Tests

stan/math/test/unit/math/prim/prob/wiener_full_test.cpp and wiener_full_prec_test.cpp
In this discourse thread we asked on how to test the new function. Since we have more parameters than all other functions, we wrote own unit tests to test all exceptions for invalid input, and tested whether valid scalar types give the correct result. Results are compared with values from the R-package WienR.

Side Effects

We had to adapt the operands_and_partials routine such that it accepts 8 parameters.
We had to implement a new, multi-dimensional integrator: hcubature.

Release notes

Two new functions in the stan::math namespace availabe: wiener_full_lpdf and wiener_full_prec_lpdf.
Along with some helper functions in stan::math::internal.
New integrator: hcubature.
Adaptation of operands_and_partials to 8 parameters.

Checklist

@spinkney
Copy link
Collaborator

I believe you'll want @SteveBronder, @andrjohns, or @rok-cesnovar to have a closer look. I think we'll want to discuss what's the best name for these.

I'm not an expert in these but I'm curious how the 5 parameter one different than the wiener_lpdf already in stan? It's 4 parameters but then we separate in a model whether the outcome hit the lower or upper bound.

@Franzi2114
Copy link
Collaborator Author

I'm not an expert in these but I'm curious how the 5 parameter one different than the wiener_lpdf already in stan? It's 4 parameters but then we separate in a model whether the outcome hit the lower or upper bound.

The 5 parameter model incorporates the variability in the drift rate parameter. This is not possible in wiener_lpdf. You can set this variability to 0, then you obtain the 4 parameter model which has the same output as wiener_lpdf, but since we have all partial derivatives implemented, the whole procedure iscomputationally much more efficient. We also need the wiener5 functions for the calculation of the wiener7 functions.

@andrjohns
Copy link
Collaborator

I believe you'll want @SteveBronder, @andrjohns, or @rok-cesnovar to have a closer look. I think we'll want to discuss what's the best name for these.

I'm not an expert in these but I'm curious how the 5 parameter one different than the wiener_lpdf already in stan? It's 4 parameters but then we separate in a model whether the outcome hit the lower or upper bound.

I can review this tomorrow

@spinkney
Copy link
Collaborator

The 5 parameter model incorporates the variability in the drift rate parameter. This is not possible in wiener_lpdf. You can set this variability to 0, then you obtain the 4 parameter model which has the same output as wiener_lpdf, but since we have all partial derivatives implemented, the whole procedure is computationally much more efficient. We also need the wiener5 functions for the calculation of the wiener7 functions.

Very cool! When we expose this to Stan we can just use wiener_lpdf and have different signatures for each type. Do we even want to keep the current wiener_lpdf if this one is superior?

If you get up the energy to add the lcdfs and the rng functions, I'd love to have those in. I have a model that needs the lcdf and lccdf I've been wanting to test. I did code up the 4 parameter _rng function (#2801) in Stan but haven't had the time to put it in.

@valentin12
Copy link
Contributor

valentin12 commented Oct 20, 2022

Very cool! When we expose this to Stan we can just use wiener_lpdf and have different signatures for each type. Do we even want to keep the current wiener_lpdf if this one is superior?

Thanks a lot for your review and your ideas. Putting it together with wiener_lpdf sounds good to me. I think as well that replacing the current wiener_lpdf implementation is worth considering. @Franzi2114 or I could do a comparison to be able to put a number on the performance difference between the implementations.

@Franzi2114
Copy link
Collaborator Author

add the lcdfs and the rng functions,

Yeah, thanks @spinkney! That's true. When setting all variabilities to 0 and calling the function with wiener_full_lpdf(a,v,w,t0,0,0,0) the four parameter model will be computed. Therefore, the wiener_lpdf() gives no advantages anymore.

We have not yet implemented the lcdf- and rng-functions. We have to check how much effort it is to bring them to Stan. But we will check this and may contribute them later in a separate PR. First, we wanted to add the lpdf-functions.

@spinkney
Copy link
Collaborator

We have not yet implemented the lcdf- and rng-functions. We have to check how much effort it is to bring them to Stan. But we will check this and may contribute them later in a separate PR. First, we wanted to add the lpdf-functions.

Yes, I agree! It's best to do these in separate prs.

@andrjohns
Copy link
Collaborator

@Franzi2114 from my understanding this PR comprises three main contributions:

  • Extending operands_and_partials to 8 parameters
  • Introducing hcubature integration routine
  • Introducing wiener_full_lpdf and wiener_full_prec_lpdf

Can you split this into three PRs rather than all together? That will make it much easier to ensure that each component is fully tested, and also gives a bit of a cleaner development history for the repo. Thanks!

@Franzi2114
Copy link
Collaborator Author

Franzi2114 commented Oct 24, 2022

Can you split this into three PRs rather than all together?

Hi @andrjohns, thanks for your input! Yes, we can split our PR into three PRs. I will refer to the same issue in the branchnames (Feature/issue 2682 operands and partials 8, and Feature/issue2682 hcubature).

Up to now, we don't have separate tests for hcubature. The integration is checked in the tests for wiener_full_lpdf(). So we will add separate tests and do a third PR.

PR for operands_and_partials: #2833
PR for hcubature #2838

@Franzi2114
Copy link
Collaborator Author

Hey @andrjohns, I am confused. I only merged the develop branch (where the new operands_and_partials_8 routine is integrated now) into this branch and the continuous-integration test fails at a point I do not understand. Did I do anything wrong?

@andrjohns
Copy link
Collaborator

@serban-nicusor-toptal would you mind restarting the Threading tests for this Jenkins run?

Is it possible to give my account permissions to restart runs as well? So you don't have to get bothered with pings everytime

@Franzi2114
Copy link
Collaborator Author

https://github.com/stan-dev/math/pull/2822/files#diff-e6ffa5dc854b843b3ee3c3c28f8eae2f436c2df2b1ca299cca1fa5982e390cf8

If you copy paste the one in develop to the one in this PR I think that should be fine

Ah, sure. I didn't feel addressed from your comment above, as I didn't touch the Jenkinsfile at all. I now copied it from the develop branch.

@stan-buildbot
Copy link
Contributor


Name Old Result New Result Ratio Performance change( 1 - new / old )
arma/arma.stan 0.19 0.2 0.95 -5.19% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan 0.01 0.01 1.03 2.91% faster
gp_regr/gen_gp_data.stan 0.02 0.03 0.71 -40.27% slower
gp_regr/gp_regr.stan 0.11 0.11 0.98 -1.96% slower
sir/sir.stan 81.3 84.22 0.97 -3.59% slower
irt_2pl/irt_2pl.stan 4.72 4.36 1.08 7.52% faster
eight_schools/eight_schools.stan 0.05 0.06 0.97 -3.52% slower
pkpd/sim_one_comp_mm_elim_abs.stan 0.37 0.27 1.38 27.31% faster
pkpd/one_comp_mm_elim_abs.stan 20.74 19.4 1.07 6.48% faster
garch/garch.stan 0.47 0.49 0.96 -4.54% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan 3.47 3.14 1.1 9.36% faster
arK/arK.stan 2.09 1.71 1.22 17.95% faster
gp_pois_regr/gp_pois_regr.stan 2.73 2.66 1.02 2.34% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan 10.96 9.85 1.11 10.08% faster
performance.compilation 201.49 215.5 0.93 -6.95% slower
Mean result: 1.0321114703119478

Jenkins Console Log
Blue Ocean
Commit hash: f70fb84db501edcbf822691c02664c43b76287b2


Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal

CPU:
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
Address sizes: 46 bits physical, 48 bits virtual
CPU(s): 80
On-line CPU(s) list: 0-79
Thread(s) per core: 2
Core(s) per socket: 20
Socket(s): 2
NUMA node(s): 2
Vendor ID: GenuineIntel
CPU family: 6
Model: 85
Model name: Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz
Stepping: 4
CPU MHz: 2400.000
CPU max MHz: 3700.0000
CPU min MHz: 1000.0000
BogoMIPS: 4800.00
Virtualization: VT-x
L1d cache: 1.3 MiB
L1i cache: 1.3 MiB
L2 cache: 40 MiB
L3 cache: 55 MiB
NUMA node0 CPU(s): 0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78
NUMA node1 CPU(s): 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79
Vulnerability Gather data sampling: Mitigation; Microcode
Vulnerability Itlb multihit: KVM: Mitigation: VMX disabled
Vulnerability L1tf: Mitigation; PTE Inversion; VMX conditional cache flushes, SMT vulnerable
Vulnerability Mds: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Meltdown: Mitigation; PTI
Vulnerability Mmio stale data: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Retbleed: Mitigation; IBRS
Vulnerability Spec rstack overflow: Not affected
Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1: Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2: Mitigation; IBRS, IBPB conditional, STIBP conditional, RSB filling, PBRSB-eIBRS Not affected
Vulnerability Srbds: Not affected
Vulnerability Tsx async abort: Mitigation; Clear CPU buffers; SMT vulnerable
Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb cat_l3 cdp_l3 invpcid_single pti intel_ppin ssbd mba ibrs ibpb stibp tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm cqm mpx rdt_a avx512f avx512dq rdseed adx smap clflushopt clwb intel_pt avx512cd avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local dtherm ida arat pln pts hwp hwp_act_window hwp_epp hwp_pkg_req pku ospke md_clear flush_l1d arch_capabilities

G++:
g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:
clang version 10.0.0-4ubuntu1
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

Copy link
Collaborator

@andrjohns andrjohns left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure why there are Jenkinsfile changes still showing in the diff, since I've git-reverted the file to match develop (and the changes in the diff are just the current state of develop). Bizarre.

Anyway, saw that there was an outstanding review comment from me, so marking this as approved. Looks like this is all good to go!

Feel free to merge once the current test run finishes up. Thanks for sticking through this marathon PR, your contribution (and patience) is very much appreciated!

@stan-buildbot
Copy link
Contributor


Name Old Result New Result Ratio Performance change( 1 - new / old )
arma/arma.stan 0.19 0.18 1.05 4.85% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan 0.01 0.01 1.04 3.77% faster
gp_regr/gen_gp_data.stan 0.02 0.02 1.08 7.07% faster
gp_regr/gp_regr.stan 0.11 0.1 1.09 8.08% faster
sir/sir.stan 77.44 76.06 1.02 1.78% faster
irt_2pl/irt_2pl.stan 3.9 3.72 1.05 4.78% faster
eight_schools/eight_schools.stan 0.05 0.05 0.98 -2.29% slower
pkpd/sim_one_comp_mm_elim_abs.stan 0.25 0.24 1.04 4.3% faster
pkpd/one_comp_mm_elim_abs.stan 18.16 17.53 1.04 3.48% faster
garch/garch.stan 0.46 0.43 1.07 6.41% faster
low_dim_gauss_mix/low_dim_gauss_mix.stan 2.79 2.69 1.03 3.36% faster
arK/arK.stan 1.63 1.57 1.04 3.56% faster
gp_pois_regr/gp_pois_regr.stan 2.49 2.5 0.99 -0.78% slower
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan 9.06 9.03 1.0 0.32% faster
performance.compilation 169.86 174.26 0.97 -2.59% slower
Mean result: 1.0327686550929835

Jenkins Console Log
Blue Ocean
Commit hash: f70fb84db501edcbf822691c02664c43b76287b2


Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal

CPU:
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
Address sizes: 46 bits physical, 48 bits virtual
CPU(s): 80
On-line CPU(s) list: 0-79
Thread(s) per core: 2
Core(s) per socket: 20
Socket(s): 2
NUMA node(s): 2
Vendor ID: GenuineIntel
CPU family: 6
Model: 85
Model name: Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz
Stepping: 4
CPU MHz: 2400.000
CPU max MHz: 3700.0000
CPU min MHz: 1000.0000
BogoMIPS: 4800.00
Virtualization: VT-x
L1d cache: 1.3 MiB
L1i cache: 1.3 MiB
L2 cache: 40 MiB
L3 cache: 55 MiB
NUMA node0 CPU(s): 0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78
NUMA node1 CPU(s): 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79
Vulnerability Gather data sampling: Mitigation; Microcode
Vulnerability Itlb multihit: KVM: Mitigation: VMX disabled
Vulnerability L1tf: Mitigation; PTE Inversion; VMX conditional cache flushes, SMT vulnerable
Vulnerability Mds: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Meltdown: Mitigation; PTI
Vulnerability Mmio stale data: Mitigation; Clear CPU buffers; SMT vulnerable
Vulnerability Retbleed: Mitigation; IBRS
Vulnerability Spec rstack overflow: Not affected
Vulnerability Spec store bypass: Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1: Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2: Mitigation; IBRS, IBPB conditional, STIBP conditional, RSB filling, PBRSB-eIBRS Not affected
Vulnerability Srbds: Not affected
Vulnerability Tsx async abort: Mitigation; Clear CPU buffers; SMT vulnerable
Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid dca sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb cat_l3 cdp_l3 invpcid_single pti intel_ppin ssbd mba ibrs ibpb stibp tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm cqm mpx rdt_a avx512f avx512dq rdseed adx smap clflushopt clwb intel_pt avx512cd avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local dtherm ida arat pln pts hwp hwp_act_window hwp_epp hwp_pkg_req pku ospke md_clear flush_l1d arch_capabilities

G++:
g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:
clang version 10.0.0-4ubuntu1
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

@SteveBronder
Copy link
Collaborator

@Franzi2114 feel free to click the merge button at your convenience :)

Your contribution is very appreciated! When you have a moment would you mind writing the Stan docs for the new 5 and 7 parameter wiener distributions in the page below?

https://github.com/stan-dev/docs/blob/master/src/functions-reference/positive_lower-bounded_distributions.qmd#L173

@Franzi2114
Copy link
Collaborator Author

@Franzi2114 feel free to click the merge button at your convenience :)

Great! Wow, I am so happy that the PR is eventually completed! ^^

However, I cannot press the merge button as I am not authorized to do this. So, could one of you merge the PR?

@spinkney
Copy link
Collaborator

spinkney commented Mar 8, 2024

I think we can add @Franzi2114 to let him do the honors himself

@SteveBronder
Copy link
Collaborator

Oh! I'm terribly sorry! idk when we changed that? I think it's nice for contributors to click merge on their own PRs. Let me see tmrw how to change that

@SteveBronder
Copy link
Collaborator

@Franzi2114 I just made you a member of the Stan math repo :) You may have to check your email to turn it on. But you should have merge permissions now!

@Franzi2114
Copy link
Collaborator Author

Hm ... I accepted the invitation and are now member of the repo. I reloaded this page and signed out and in again. But still the message is 'You are not authorized to merge this pull request.' Seems more tricky to do the last step than expected ^^

@SteveBronder
Copy link
Collaborator

@Franzi2114 sorry try one more time. Thought I gave you the right access but it was one less than you needed

@Franzi2114
Copy link
Collaborator Author

Franzi2114 commented Mar 12, 2024

Now there popped up another box with a rocket saying 'This branch has not been deployed', but I still cannot merge the PR.

@SteveBronder
Copy link
Collaborator

My apologies the whole thing is confusing. I just sent you a member request for the Stan dev team which should give you sufficient privileges (check your email). If after clicking the email you still cannot merge the PR then I'll merge it.

I find this really weird though and want to figure out why this is happening. I really think it's important and satisfying to let collaborators merge their own pull requests (first PR is a big deal!)

@SteveBronder
Copy link
Collaborator

@serban-nicusor-toptal do you know how to make it so outside collaborators can still click the merge button on their own PRs?

@serban-nicusor-toptal
Copy link
Contributor

serban-nicusor-toptal commented Mar 13, 2024

Good morning @SteveBronder @Franzi2114, I think everything is setup correctly. In order for Franzi to get proper access, he will need to accept the invitation to join the Math org team (this team is mapped in the branch protections). Once that's done, the team will grant you the required permissions.
image

@Franzi2114 Franzi2114 merged commit b21e184 into stan-dev:develop Mar 13, 2024
8 checks passed
@Franzi2114
Copy link
Collaborator Author

YEEEEEEEEEEYYYYYYYYYY!!!!!!!!!! It worked!!!!!!!!!!!!!!!!!!!

@Franzi2114
Copy link
Collaborator Author

Wuhu! Thanks =D
Wow, that was really a long journey, but I think I learned a lot.

Thank you so much for all your proof reading, input, corrections and all the time you spent on this PR!

I already prepared the CDF functions for the 7-parameter diffusion model. I will make the PR in the next days.

@WardBrian
Copy link
Member

Congrats on the huge effort!

We can expose the basic scalar version in the compiler very soon (most likely next week!). A more advanced typechecking which lets us expose the vectorized version is something we’ve started discussing as well

@spinkney
Copy link
Collaborator

Congrats! Looking forward to testing this out and seeing the cdfs. Are you planning to implement the rng functions too?

@Franzi2114 Franzi2114 deleted the feature/issue-2682-Add-7-parameter-DDM-PDF branch March 15, 2024 08:34
@Franzi2114
Copy link
Collaborator Author

rng functions too?

First, I do the CDF and CCDF functions. But compared to the CDF and the PDF, the rng should be quickly done. So I'll probably implement this as well.

@SteveBronder
Copy link
Collaborator

In terms of use I think the rngs would be really nice so folks can do prediction etc.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet