`fast_pdf` defaults can create likelihood approximation singularities

[mshvartsman]

I ran into what I think is unintended behavior in the density function in fast_pdf.

Specifically, fast_1d_kde defaults to only covering the range of the input data, and np.interp defaults to returning the edge value if queried past the edge. Both are set to their defaults in fast_pdf.

Suppose that your model can get close to a point mass somewhere over the support of the data (e.g. a DDM with a very low threshold). Then the density at that one point is pretty high. fast_1d_kde then defaults to only computing the density over that very narrow range and then np.interp defaults to returning that high density value for all the rest of the data (which lives outside the bounds of the grid the KDE function returned). Then a reasonable optimizer would keep tightening that point mass and getting ever higher likelihoods.

Below's a minimal (extreme) example, though I've originally seen this come up in a real model.

from RunDEMC.density import fast_pdf
import numpy as np 
model = np.array([1.00001,1,1,1,1], dtype=np.float64)
data =  np.array([1,2,3,4,5,6,7,8,9,10],dtype=np.float64)
print(fast_pdf(model, data, 200))  # probably not intended
# [ 98667.44131307  25479.89240456  25479.89240456  25479.89240456  25479.89240456 25479.89240456  25479.89240456  25479.89240456 25479.89240456  25479.89240456]
print(fast_pdf(model, data, 200, (0, 12)))  # looks better
# [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]

An easy fix is either default the extrema to the range of (min(data,model), max(data,model) or call np.interp with left=0, right=0. I'm happy to submit a pull request with either option.

[psederberg]

Thanks for looking into this! Please submit a PR with whichever option you
prefer. Note that I haven't actually been using fast_pdf, instead I've been
using kdensity because I could verify that it gave identical answers to R's
density functions.

Have you compared the speed of the two? It might be the case that kdensity
isn't much slower, but I'd definitely be interested in knowing either way.

On Thu, Mar 17, 2016 at 10:14 PM, Michael Shvartsman <
notifications@github.com> wrote:

I ran into what I think is unintended behavior in the density function in
fast_pdf.

Specifically, fast_1d_kde defaults to only covering the range of the
input data, and np.interp defaults to returning the edge value if queried
past the edge. Both are set to their defaults in fast_pdf.

Suppose that your model can get close to a point mass somewhere over the
support of the data (e.g. a DDM with a very low threshold). Then the
density at that one point is pretty high. fast_1d_kde then defaults to
only computing the density over that very narrow range and then np.interp
defaults to returning that high density value for all the rest of the data
(which lives outside the bounds of the grid the KDE function returned).
Then a reasonable optimizer would keep tightening that point mass and
getting ever higher likelihoods.

Below's a minimal (extreme) example, though I've originally seen this come
up in a real model.

from RunDEMC.density import fast_pdfimport numpy as np
model = np.array([1.00001,1,1,1,1], dtype=np.float64)
data = np.array([1,2,3,4,5,6,7,8,9,10],dtype=np.float64)print(fast_pdf(model, data, 200)) # probably not intendedprint(fast_pdf(model, data, 200, (0, 12))) # looks better

An easy fix is either default the extrema to the range of
(min(data,model), max(data,model) or call np.interp with left=0, right=0.
I'm happy to submit a pull request with either option.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
#1

[psederberg]

If you don't want to bother with the PR, I'm happy to just make the change
given it's a single line :)

I think I'll change the interp call with left and right...

On Thu, Mar 17, 2016 at 10:27 PM, Per B. Sederberg psederberg@gmail.com
wrote:

Thanks for looking into this! Please submit a PR with whichever option you
prefer. Note that I haven't actually been using fast_pdf, instead I've been
using kdensity because I could verify that it gave identical answers to R's
density functions.

Have you compared the speed of the two? It might be the case that kdensity
isn't much slower, but I'd definitely be interested in knowing either way.

On Thu, Mar 17, 2016 at 10:14 PM, Michael Shvartsman <
notifications@github.com> wrote:

I ran into what I think is unintended behavior in the density function in
fast_pdf.

Specifically, fast_1d_kde defaults to only covering the range of the
input data, and np.interp defaults to returning the edge value if
queried past the edge. Both are set to their defaults in fast_pdf.

Suppose that your model can get close to a point mass somewhere over the
support of the data (e.g. a DDM with a very low threshold). Then the
density at that one point is pretty high. fast_1d_kde then defaults to
only computing the density over that very narrow range and then np.interp
defaults to returning that high density value for all the rest of the data
(which lives outside the bounds of the grid the KDE function returned).
Then a reasonable optimizer would keep tightening that point mass and
getting ever higher likelihoods.

Below's a minimal (extreme) example, though I've originally seen this
come up in a real model.

from RunDEMC.density import fast_pdfimport numpy as np
model = np.array([1.00001,1,1,1,1], dtype=np.float64)
data = np.array([1,2,3,4,5,6,7,8,9,10],dtype=np.float64)print(fast_pdf(model, data, 200)) # probably not intendedprint(fast_pdf(model, data, 200, (0, 12))) # looks better

An easy fix is either default the extrema to the range of
(min(data,model), max(data,model) or call np.interp with left=0, right=0.
I'm happy to submit a pull request with either option.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
#1

[mshvartsman]

My bottleneck at the moment is the simulation, so I haven't compared speed of the density estimation yet (though I know that the naive way of doing it was expensive in my last implementation).

Agreed that the interp call sounds more elegant (otherwise you have to pay for min/max calls on each call to fast_pdf). But on reflection: is there a case (e.g. in very long-tailed model distributions) where you want the density to fall off more gradually? I think usually it shouldn't be an issue because you should have enough model data to make your bandwidth small enough that the density is effectively 0 outside the model's simulated support. But for small enough, widely spread-out data (e.g. when the error distribution is wide but error rate is low) you might underestimate the density?