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

mrc-4324: support for taking expectations of stochastic functions #296

Merged
merged 7 commits into from
Jun 26, 2023

Conversation

richfitz
Copy link
Member

@richfitz richfitz commented Jun 13, 2023

Required for differentiating models that contain stochastic distributions. Normally we cope with this in dust by putting distributions into "deterministic mode" which applies the same logic. Here, we do the transformation of a call to a stochastic function (e.g., rnorm(mu, sd)) with its expectation (mu) symbolically so that we can later differentiate the resulting expressions.

No visible change from any of this to the user, this is just support code to make the eventual autodiff support manageable

@codecov
Copy link

codecov bot commented Jun 13, 2023

Codecov Report

Patch coverage: 100.00% and no project coverage change.

Comparison is base (e069323) 100.00% compared to head (2ef5f87) 100.00%.

❗ Current head 2ef5f87 differs from pull request most recent head 791103d. Consider uploading reports for the commit 791103d to get more accurate results

Additional details and impacted files
@@            Coverage Diff            @@
##            master      #296   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           46        47    +1     
  Lines         5502      5509    +7     
=========================================
+ Hits          5502      5509    +7     
Impacted Files Coverage Δ
R/differentiate-support.R 100.00% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

@richfitz richfitz marked this pull request as ready for review June 13, 2023 17:50
Copy link
Contributor

@M-Kusumgar M-Kusumgar left a comment

Choose a reason for hiding this comment

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

Nicely done, was fun to just look through the list of distributions and either try and remember their mean or google it 🙂

Comment on lines +50 to +52
rgeom = function(expr) {
substitute((1 - p) / p, list(p = expr[[2]]))
},
Copy link
Contributor

Choose a reason for hiding this comment

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

I might be being a bit dumb, rgeom takes the probability of success right, so should this not be 1/p?

Copy link
Contributor

Choose a reason for hiding this comment

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

i realise you're making these functions yourself so could tell people that you need to put the probability of failure in the rgeom function but would be nice to have consistent behaviour with the rgeom from R function

Copy link
Member Author

Choose a reason for hiding this comment

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

This is consistent with rgeom (see tests); the actual random variable here is the number of observed failures before a success is observed:

    x, q: vector of quantiles representing the number of failures in a
          sequence of Bernoulli trials before success occurs.

so in the edge case of pr = 1, we have 0 / 1 = 0 and rgeom(n, 1) is 0 (repeated n times).

Or am I misunderstanding?

Copy link
Contributor

Choose a reason for hiding this comment

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

ahhh apologies, I got confused with what the R function took, I thought it took a random variable, X, which was the number of Bernoulli trials to get one success, just tested it out numerically too! rgeom is definitely consistent with your expectation!

Comment on lines +82 to +84
rwilcox = function(expr) {
substitute(m * n / 2, list(m = expr[[2]], n = expr[[3]]))
},
Copy link
Contributor

Choose a reason for hiding this comment

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

For anyone interested, result is here in the "Normal approximation and tie correction" section

Comment on lines +85 to +87
rsignrank = function(expr) {
substitute(n * (n + 1) / 4, list(n = expr[[2]]))
})
Copy link
Contributor

Choose a reason for hiding this comment

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

Again for those interested, result is here in the "Computing the null distribution" section


## Discrete distrbutions are somewhat harder
expectation_discrete <- function(fd, fq, pars, tol = 1e-12) {
end <- do.call(fq, c(list(p = 1 - tol), unname(pars)))
Copy link
Contributor

Choose a reason for hiding this comment

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

again maybe this is very standard but think this could use a tiny comment to explain fq which i believe is the quantile function which will find an end value such that our expectation is an average over 1 - tol area under the graph, probably a way to word it better than that!

Copy link
Contributor

Choose a reason for hiding this comment

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

Agreed with this, and maybe a link to the definitions. To make this really obvious if we come back to this in the future.

Copy link
Contributor

@r-ash r-ash left a comment

Choose a reason for hiding this comment

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

This looks good, one comment about the TODO

substitute(1 / rate, list(rate = expr[[2]]))
},
rf = function(expr) {
## TODO: only valid for df2 > 2!
Copy link
Contributor

Choose a reason for hiding this comment

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

Do you want to add some asserts to check these edge cases?


## Discrete distrbutions are somewhat harder
expectation_discrete <- function(fd, fq, pars, tol = 1e-12) {
end <- do.call(fq, c(list(p = 1 - tol), unname(pars)))
Copy link
Contributor

Choose a reason for hiding this comment

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

Agreed with this, and maybe a link to the definitions. To make this really obvious if we come back to this in the future.

Copy link
Contributor

@weshinsley weshinsley left a comment

Choose a reason for hiding this comment

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

Ok - with discussion with the team, I think I understand this one, and looks good - version bump when we get to merging it

@@ -1,6 +1,6 @@
Package: odin
Title: ODE Generation and Integration
Version: 1.5.0
Version: 1.5.1
Copy link
Contributor

Choose a reason for hiding this comment

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

Bump after PR #295 is merged

@@ -813,7 +813,7 @@ test_that_odin("overlapping graph", {
list(y * p, p + p2)
}
cmp <- deSolve::ode(1, tt, f, NULL)
tol <- variable_tolerance(mod, js = 1e-6)
tol <- variable_tolerance(mod, js = 3e-6)
Copy link
Contributor

Choose a reason for hiding this comment

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

This will go away when #295 is merged

@richfitz richfitz merged commit 1653828 into master Jun 26, 2023
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.

4 participants