Add convex decreasing density estimator for π0 #56
Conversation
juliangehring
force-pushed
the
convex-decreasing-pi0
branch
from
60c134b
to
fa0e6fd
Compare
src/pi0-estimators.jl
Outdated
| end | ||
|
|
||
| function find_f_theta(x::Vector{Float64}, theta::Float64) | ||
| return @. 2 / theta^2 * (theta-x) * (x.<theta) |
There was a problem hiding this comment.
What does the @. notation mean?
There was a problem hiding this comment.
This applies the dot vectorisation to all functions on the right, and turns the line into 2 ./ theta^2 .* (theta.-x) .* (x.<theta). I find it a bit easier to read without all the extra dots.
There was a problem hiding this comment.
Cool thanks, that seems very useful! (Then there's one . too much in the code though).
There was a problem hiding this comment.
Thanks for spotting this, fixed in 2966dba.
src/pi0-estimators.jl
Outdated
| @. f = (1-ε)*f + ε*f_theta | ||
| @. f_p = @views gridsize*(f[idx_lower]-f[idx_upper])*px + f[idx_upper] | ||
| theta = dx * find_theta(t, p, f_p) | ||
| f_theta .= find_f_theta(x, theta) |
There was a problem hiding this comment.
For the sake of code readability, instead of calling this find_f_theta, can we maybe use Distributions.jl (TriangularDist) and use pdf?
There was a problem hiding this comment.
I'm struggling a bit with TriangularDist since the weights at x=0 always yield 0:
pdf.(TriangularDist(0.0, 1.0, 0.0), 0.0) (I would have expected it to yield the same as pdf.(TriangularDist(0.0, 1.0, 1.0), 1.0)). Is this a bug or a feature in Distributions?
There was a problem hiding this comment.
Oh interesting.. I'd say it is almost surely a bug with how they handle the boundary..
There was a problem hiding this comment.
Oh good, then it's not just me ;)
There was a problem hiding this comment.
I have given the function a more meaningful name in 2966dba. Since the TriangularDist in Distributions seems buggy, I'd keep our function for now.
src/pi0-estimators.jl
Outdated
| f_theta .= find_f_theta(x, theta) | ||
| pi0_new = f[end] | ||
| if abs(pi0_new - pi0_old) <= xtol | ||
| return pi0_new, thetas, true |
There was a problem hiding this comment.
I think the paper suggests a stopping rule based on sum((f_p.-f_theta_p)./f_p), do these criteria behave similarly?
There was a problem hiding this comment.
They mention the stopping rule sum((f_p.-f_theta_p)./f_p) < -δ in the algorithm, but don't use it in their reference implementation - they don't use any stopping rule there and simply stop after 100 iterations. Since I'm not sure how to choose a good value for δ in this case, I have used a general stopping rule based on the change of the π0 estimate as in the CensoredBUM estimator (#4).
src/pi0-estimators.jl
Outdated
| end | ||
| end | ||
| @. f = (1-ε)*f + ε*f_theta | ||
| @. f_p = @views gridsize*(f[idx_lower]-f[idx_upper])*px + f[idx_upper] |
There was a problem hiding this comment.
Can we get f_p exactly by updating the previous one instead of this linear interpolation?
There was a problem hiding this comment.
What do you have in mind for getting f_p exactly?
There was a problem hiding this comment.
If I see it correctly, f is the same as f_p, except evaluated on the sorted p-values rather than the grid. So e.g. above could update it with something like f_p = (1-ε)*f_p+ ε*f_theta_p (assuming properly initialized in the beginning). But maybe it's not worth doing this exactly also from computational perspective?
There was a problem hiding this comment.
I have run some simulations with both update rules, and I could never see a difference in the final π0 estimate. Your solution looks nicer and cleaner to me, I have added it in 2966dba.
|
Except the 3 comments, looks very good to me; this is very cool! With this + Grenander + isotonic, I wonder if it might be worth to move them to a |
|
@nignatiadis Do you think this needs more work, or are you fine with merging this feature? |
|
Hey Julian, sorry for delaying this. It looks good to me! |
|
No problem - thanks for having a look at it. |
Implements the convex decreasing density estimator for π0 (#3 ), based on section 4.3 in