inverse function of a spline
I have a function y = f(x) which is monotonically increasing (a cumulative distribution function) f is defined by piecewise polynomial interpolation, an interpolating spline on some points I would like to get the inverse function (ppf) x = f^{-1} (y) if the spline is of higher order than linear In the linear case it's trivial, because the inverse function is also just a piecewise linear interpolation. If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly? I don't know much about general properties of splines, and would appreciate any hints, so I can avoid numerical inversion (fsolve or similar) Thanks, Josef
On Fri, May 7, 2010 at 10:40 AM, <josef.pktd@gmail.com> wrote:
I have a function y = f(x) which is monotonically increasing (a cumulative distribution function) f is defined by piecewise polynomial interpolation, an interpolating spline on some points
I would like to get the inverse function (ppf) x = f^{-1} (y) if the spline is of higher order than linear
In the linear case it's trivial, because the inverse function is also just a piecewise linear interpolation.
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
I don't know much about general properties of splines, and would appreciate any hints, so I can avoid numerical inversion (fsolve or similar)
Since the curve is piecewise cubic the problem reduces to inverting a piece of a cubic, which inverse won't itself be a cubic in general. I think your best bet is interpolate the same points with x,y reversed, or resample using your spline and interpolate the new samples with x,y reversed. It won't be a exact inverse, but then, the original is probably not exact either. Chuck
On Fri, May 7, 2010 at 1:57 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Fri, May 7, 2010 at 10:40 AM, <josef.pktd@gmail.com> wrote:
I have a function y = f(x) which is monotonically increasing (a cumulative distribution function) f is defined by piecewise polynomial interpolation, an interpolating spline on some points
I would like to get the inverse function (ppf) x = f^{-1} (y) if the spline is of higher order than linear
In the linear case it's trivial, because the inverse function is also just a piecewise linear interpolation.
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
I don't know much about general properties of splines, and would appreciate any hints, so I can avoid numerical inversion (fsolve or similar)
Since the curve is piecewise cubic the problem reduces to inverting a piece of a cubic, which inverse won't itself be a cubic in general. I think your best bet is interpolate the same points with x,y reversed, or resample using your spline and interpolate the new samples with x,y reversed. It won't be a exact inverse, but then, the original is probably not exact either.
That's what I suspected, I was hoping for a trick (like one interpolator is the "natural" inverse of another one). resampling should give a good enough approximation. Without resampling, the error for round tripping x= f^{-1} ( f(x) ) might be too large to give consistent results. (Even if there are sampling and approximation errors, I still would prefer consistency.) Just a follow-up question on approximation: for the cdf (e.g. normal distribution) f:R->[0,1] f^{-1}:[0,1]->R Is it better to start with a spline on the inverse function (ppf), f^{-1}, because it has compact support, resample from it, and then create the cdf f from a resampled ppf; or the other way around, or it wouldn't really matter? Thanks for the information and hint, Josef
Chuck
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Fri, May 7, 2010 at 12:34 PM, <josef.pktd@gmail.com> wrote:
On Fri, May 7, 2010 at 1:57 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Fri, May 7, 2010 at 10:40 AM, <josef.pktd@gmail.com> wrote:
I have a function y = f(x) which is monotonically increasing (a cumulative distribution function) f is defined by piecewise polynomial interpolation, an interpolating spline on some points
I would like to get the inverse function (ppf) x = f^{-1} (y) if the spline is of higher order than linear
In the linear case it's trivial, because the inverse function is also just a piecewise linear interpolation.
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
I don't know much about general properties of splines, and would appreciate any hints, so I can avoid numerical inversion (fsolve or similar)
Since the curve is piecewise cubic the problem reduces to inverting a
piece
of a cubic, which inverse won't itself be a cubic in general. I think your best bet is interpolate the same points with x,y reversed, or resample using your spline and interpolate the new samples with x,y reversed. It won't be a exact inverse, but then, the original is probably not exact either.
That's what I suspected, I was hoping for a trick (like one interpolator is the "natural" inverse of another one).
resampling should give a good enough approximation. Without resampling, the error for round tripping x= f^{-1} ( f(x) ) might be too large to give consistent results. (Even if there are sampling and approximation errors, I still would prefer consistency.)
Just a follow-up question on approximation:
for the cdf (e.g. normal distribution) f:R->[0,1] f^{-1}:[0,1]->R
Is it better to start with a spline on the inverse function (ppf), f^{-1}, because it has compact support, resample from it, and then create the cdf f from a resampled ppf; or the other way around, or it wouldn't really matter?
Thanks for the information and hint,
I don't know the answer to that although I suspect starting from the inverse might be superior. The end points might be a problem though if the curve goes vertical. You might have to experiment a bit or use a spline in combination with other functions. There has probably been a small industry out there dealing with these sorts of problem but I don't know who they are.
Chuck
On 7 May 2010 12:40, <josef.pktd@gmail.com> wrote:
I have a function y = f(x) which is monotonically increasing (a cumulative distribution function) f is defined by piecewise polynomial interpolation, an interpolating spline on some points
I would like to get the inverse function (ppf) x = f^{-1} (y) if the spline is of higher order than linear
In the linear case it's trivial, because the inverse function is also just a piecewise linear interpolation.
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
I don't know much about general properties of splines, and would appreciate any hints, so I can avoid numerical inversion (fsolve or similar)
I should first say that even though your input points are monotonic, the spline is not guaranteed to be. (Though in practice if your sampled points have no sharp corners you're probably fine.) If this matters to you, there are algorithms for enforcing monotonicity of splines, some of which are simply procedures for jiggering the interpolation points just enough to avoid non-monotonicty and some of which are more clever. Sadly none are implemented in scipy. As Charles Harris pointed out, the inverse of a cubic is not a cubic, so the inverse function won't be a spline. But you can nevertheless efficiently evaluate it with scipy.interpolate.sproot, which a special-purpose numerical solver. I'm not sure whether it uses cubic root solvers or an optimized numerical solver with knowledge about the bounding properties of spline control points, but in any case it's quite efficient. It only finds zeros, but (check this) you should be able to shift a spline vertically by subtracting a constant from the coefficient array (c in t,c,k). Since you are constructing the spline in the first place, you should also think about whether you're evaluating f or f inverse more often and choose which one to be the spline appropriately. Anne
On Fri, May 7, 2010 at 4:24 PM, Anne Archibald <peridot.faceted@gmail.com> wrote:
On 7 May 2010 12:40, <josef.pktd@gmail.com> wrote:
I have a function y = f(x) which is monotonically increasing (a cumulative distribution function) f is defined by piecewise polynomial interpolation, an interpolating spline on some points
I would like to get the inverse function (ppf) x = f^{-1} (y) if the spline is of higher order than linear
In the linear case it's trivial, because the inverse function is also just a piecewise linear interpolation.
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
I don't know much about general properties of splines, and would appreciate any hints, so I can avoid numerical inversion (fsolve or similar)
I should first say that even though your input points are monotonic, the spline is not guaranteed to be. (Though in practice if your sampled points have no sharp corners you're probably fine.) If this matters to you, there are algorithms for enforcing monotonicity of splines, some of which are simply procedures for jiggering the interpolation points just enough to avoid non-monotonicty and some of which are more clever. Sadly none are implemented in scipy.
As Charles Harris pointed out, the inverse of a cubic is not a cubic, so the inverse function won't be a spline. But you can nevertheless efficiently evaluate it with scipy.interpolate.sproot, which a special-purpose numerical solver. I'm not sure whether it uses cubic root solvers or an optimized numerical solver with knowledge about the bounding properties of spline control points, but in any case it's quite efficient. It only finds zeros, but (check this) you should be able to shift a spline vertically by subtracting a constant from the coefficient array (c in t,c,k).
Since you are constructing the spline in the first place, you should also think about whether you're evaluating f or f inverse more often and choose which one to be the spline appropriately.
Thanks, I will try to figure out the sproot version. For now I'm stuck (and go somewhere else) because in the examples that I tried out, I get small non-monotonicities most of the time. The spline of the inverse function is backwards bending. I will stick with linear interpolation and kernel density estimation for the smooth case. BTW: I'm writing some histogram distribution and variants of empirical distribution classes that have the same methods as the ones in scipy.stats. Josef
Anne _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
Hi Josef,
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions. Nicky
On Fri, May 7, 2010 at 4:37 PM, nicky van foreest <vanforeest@gmail.com> wrote:
Hi Josef,
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions.
Now I know it's not, but I was designing the extension to the linear case on paper instead of in the interpreter, and got stuck on the wrong problem. Maybe I ask the question again when scipy has monotonic interpolators. Josef
Nicky _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
<josef.pktd <at> gmail.com> writes:
On Fri, May 7, 2010 at 4:37 PM, nicky van foreest <vanforeest <at> gmail.com>
wrote:
Hi Josef,
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions.
Now I know it's not, but I was designing the extension to the linear case on paper instead of in the interpreter, and got stuck on the wrong problem.
There's an algorithm for making constrained-to-be-monotonic spline interpolants (only in one dimension, though). The reference is Dougherty et al 1989 Mathematics of Computation, vol 52 no 186 pp 471-494 (April 1989). This is available on-line at www.jstor.org.
On Thu, Sep 29, 2011 at 12:37 PM, Jeff Brown <brownj@seattleu.edu> wrote:
<josef.pktd <at> gmail.com> writes:
On Fri, May 7, 2010 at 4:37 PM, nicky van foreest <vanforeest <at> gmail.com>
wrote:
Hi Josef,
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions.
Now I know it's not, but I was designing the extension to the linear case on paper instead of in the interpreter, and got stuck on the wrong problem.
There's an algorithm for making constrained-to-be-monotonic spline interpolants (only in one dimension, though). The reference is Dougherty et al 1989 Mathematics of Computation, vol 52 no 186 pp 471-494 (April 1989). This is available on-line at www.jstor.org.
Thanks for the reference. Maybe Ann's interpolators in scipy that take derivatives could be used for this. Shape preserving splines or piecewise polynomials would make a nice addition to scipy, but I'm only a potential user. I have dropped this for the moment, after taking a detour with (global) orthonormal polynomial approximation, where I also haven't solved the integration and function inversion problem yet (nice pdf but only brute force cdf and ppf). Josef
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Fri, Sep 30, 2011 at 12:37 PM, <josef.pktd@gmail.com> wrote:
On Thu, Sep 29, 2011 at 12:37 PM, Jeff Brown <brownj@seattleu.edu> wrote:
<josef.pktd <at> gmail.com> writes:
On Fri, May 7, 2010 at 4:37 PM, nicky van foreest <vanforeest <at>
gmail.com> wrote:
Hi Josef,
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions.
Now I know it's not, but I was designing the extension to the linear case on paper instead of in the interpreter, and got stuck on the wrong problem.
There's an algorithm for making constrained-to-be-monotonic spline interpolants (only in one dimension, though). The reference is Dougherty et al 1989 Mathematics of Computation, vol 52 no 186 pp 471-494 (April 1989). This is available on-line at www.jstor.org.
Thanks for the reference. Maybe Ann's interpolators in scipy that take derivatives could be used for this.
trying out how PiecewisePolynomial works, almost but not quite enough Josef spamming the world with messy examples
Shape preserving splines or piecewise polynomials would make a nice addition to scipy, but I'm only a potential user.
I have dropped this for the moment, after taking a detour with (global) orthonormal polynomial approximation, where I also haven't solved the integration and function inversion problem yet (nice pdf but only brute force cdf and ppf).
Josef
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
On Sat, Oct 1, 2011 at 8:52 AM, <josef.pktd@gmail.com> wrote:
On Fri, Sep 30, 2011 at 12:37 PM, <josef.pktd@gmail.com> wrote:
On Thu, Sep 29, 2011 at 12:37 PM, Jeff Brown <brownj@seattleu.edu> wrote:
<josef.pktd <at> gmail.com> writes:
On Fri, May 7, 2010 at 4:37 PM, nicky van foreest <vanforeest <at>
gmail.com> wrote:
Hi Josef,
If I have a cubic spline, or any other smooth interpolator in scipy, is there a way to get the inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions.
Now I know it's not, but I was designing the extension to the linear case on paper instead of in the interpreter, and got stuck on the wrong problem.
There's an algorithm for making constrained-to-be-monotonic spline interpolants (only in one dimension, though). The reference is Dougherty et al 1989 Mathematics of Computation, vol 52 no 186 pp 471-494 (April 1989). This is available on-line at www.jstor.org.
Thanks for the reference. Maybe Ann's interpolators in scipy that take derivatives could be used for this.
trying out how PiecewisePolynomial works, almost but not quite enough
IIRC, de Boor dealt with fitting distribution functions somewhere in his book "A Practical Guide to Splines". I don't recall whether or not he constrains things to positivity, but recalling one of the figures, I think that he was fitting histograms, perhaps their area. <snip> Chuck
On Sat, Oct 1, 2011 at 11:17 AM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Sat, Oct 1, 2011 at 8:52 AM, <josef.pktd@gmail.com> wrote:
On Fri, Sep 30, 2011 at 12:37 PM, <josef.pktd@gmail.com> wrote:
On Thu, Sep 29, 2011 at 12:37 PM, Jeff Brown <brownj@seattleu.edu> wrote:
<josef.pktd <at> gmail.com> writes:
On Fri, May 7, 2010 at 4:37 PM, nicky van foreest <vanforeest <at>
gmail.com> wrote:
Hi Josef,
> If I have a cubic spline, or any other smooth interpolator in scipy, > is there a way to get the > inverse function directly?
How can you ensure that the cubic spline approx is non-decreasing? I actually wonder whether using cubic splines is the best way to approximate distribution functions.
Now I know it's not, but I was designing the extension to the linear case on paper instead of in the interpreter, and got stuck on the wrong problem.
There's an algorithm for making constrained-to-be-monotonic spline interpolants (only in one dimension, though). The reference is Dougherty et al 1989 Mathematics of Computation, vol 52 no 186 pp 471-494 (April 1989). This is available on-line at www.jstor.org.
Thanks for the reference. Maybe Ann's interpolators in scipy that take derivatives could be used for this.
trying out how PiecewisePolynomial works, almost but not quite enough
IIRC, de Boor dealt with fitting distribution functions somewhere in his book "A Practical Guide to Splines". I don't recall whether or not he constrains things to positivity, but recalling one of the figures, I think that he was fitting histograms, perhaps their area.
looks nice, he matches the area in each bin. page 79 ff, I don't see any explicit non-negativity constraints. as for an efficient implementation in numpython: ? Josef
<snip>
Chuck
_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
An example of a possible usecase background I wrote a function for Lilliefors test of normality (Kolmogorov Smirnov test for normal distribution based on estimated mean and variance). (The reported power of the test is not great.) The p-values are non-standard and tabulated, by sample size and by probability (only 0.001 to 0.2) To safe work later on for other tests, I wrote a simple class, but I restrict to linear interpolation. Using two interp1d or using Rbf works fine, using interp2d gives me a warning - cannot add more knots - and the interpolated values are not correct. Since I never played much with interp2d, or interpolate.bisplrep or interpolate.BivariateSpline, human error is possible. Josef
participants (5)
-
Anne Archibald -
Charles R Harris -
Jeff Brown -
josef.pktd@gmail.com -
nicky van foreest