restricting interpolating functions to be positive
dantimatter
I would like to construct an interpolation that can only take on
positive values. The data that I'm trying to interpolate is all
positive, but for some reason the Interpolation[] function 'wants' to
connect the dots by dipping below zero. Is there a way to force
positivity? Or should I simply adjust the InterpolationOrder until it
works?
Thanks,
Dan
Mark McClure
Here's a (too?) simple approach:
data = {{0, 2}, {1, 0.2}, {2, 0.2}, {3, 2}};
f = Interpolation[data];
g[x_] := Max[{f[x], 0}];
Plot[{f[x], g[x]}, {x, 1, 2},
PlotStyle -> {Thin, Thick}]
You'll lose some smoothness but you'd certainly lose smoothness by
decreasing the InterpolationOrder as well.
Mark McClure
dh
Hi,
Interpolation by default fits cubic polynomials. These tend to
over/undershoot the points given. If you can live with linear
interpolation, you may specify InterpolationOrder->0. Then the points
are simply connected by straight lines.
Daniel
Valeri Astanoff
Good day,
Here is a DIY way to do it (just a toy example) :
In[1]:= data = {{0, 2}, {1, 1}, {2, .2}, {3, .01}, {4, 1},
{5, 0.9}, {6, 0.05}, {7, 0.1}, {8, 1}};
xm = data[[-1, 1]];
[ we're gonna get negative interpolation with data close to zero ]
In[3]:= f1 = Interpolation[data, InterpolationOrder -> 1];
f2 = Interpolation[data];
[ f1 is used just to locate the minimums, an additional test
should be used below to only select minimums close to zero ]
In[5]:= dx = 0.1;
[ dx is to be adjusted accordingly ]
In[6]:= data2 = Table[{{x}, f2[x],
If[ f1'[x - dx] < 0 && f1'[x + dx] > 0 , 0, f2'[x]]}, {x, 0, xm}]
During evaluation of In[6]:= InterpolatingFunction::dmval: Input \
value {-0.1} lies outside the range of data in the interpolating \
function. Extrapolation will be used. >>
Out[6]= {{{0}, 2., -0.963333}, {{1}, 1., -0.968333}, {{2},
0.2, -0.563333}, {{3}, 0.01, 0}, {{4}, 1., 0.823333}, {{5},
0.9, -0.531667}, {{6}, 0.05, 0}, {{7}, 0.1, 0.483333}, {{8}, 1.,
1.30833}}
[ Notice derivatives at 3 and 6 set to zero ]
In[7]:= f3 = Interpolation[data2];
In[8]:= Plot[f1[x], {x, 0, xm}]
[...]
In[9]:= Plot[f2[x], {x, 0, xm}]
[...]
In[10]:= Plot[f3[x], {x, 0, xm}]
[...]
Seems okay ( sure not for all cases ! )
--
Valeri Astanoff
Bill Rowe
>I would like to construct an interpolation that can only take on
>positive values. The data that I'm trying to interpolate is all
>positive, but for some reason the Interpolation[] function 'wants'
>to connect the dots by dipping below zero.
By default, Interpolation uses piecewise cubic polynomials to
fit your data. The coefficients of these polynomials are found
by selecting a sufficient number of consecutive (4 for cubics)
and solving for the needed coefficients. The dips below zero are
a consequence of using a 3 order polynomial and the specific
values of the data you are interpolating.
>Is there a way to force positivity? Or should I simply adjust the
>InterpolationOrder until it works?
You really don't have much option here. If you reduce the
interpolating order to 1, and all of the data is positive you
are guaranteed all values of the resulting interpolation will be
positive. But this also means derivatives will not be continuous
which may or may not be a problem for you. Lack of continuity in
the first and second derivatives means a plot of the
interpolating function will not be smooth. No other
interpolation order is certain to only provide positive values
given positive data. Higher order interpolating polynomials are
more likely to dip below zero.
If you want continuous first and second derivatives (a smooth
interpolation) and only positive values, you either need to use
some other basis for interpolating your data than polynomials
such as a set of exponential terms or give up the requirement
the interpolating function pass through your data.
Using exponentials creates other problems. It transforms the
problem into one that will not have explicit algebraic solutions
for the coefficients. That means you will need iterative methods
which will only converge well if you provide sufficiently close
initial values as a starting point. Also, these methods tend not
to be very stable numerically when more than a few exponential
terms are required.
You can achieve smooth interpolating functions with cubic
polynomials using a smoothed natural cubic spline. If you are
interested in this approach, I've written a package for my use
which has the needed functions that I would be willing to send
you offline. Note, while this will yield a smooth function
(continuous first and second derivatives) it will not pass
through your data points when you require all positive values
for the range of your data.
DrMajorBob
data = Sort@RandomReal[{0, 1}, {10, 2}];
{min, max} = data[[{1, -1}, 1]];
Clear[f, g]
g = Interpolation[data];
f[x_] := Max[0, g@x]
Plot[f@x, {x, min, max}, PlotRange -> All]
Bobby
On Mon, 11 Jan 2010 17:52:33 -0600, dantimatter <goo...@dantimatter.com>
wrote:
> Hi All,
>
> I would like to construct an interpolation that can only take on
> positive values. The data that I'm trying to interpolate is all
> positive, but for some reason the Interpolation[] function 'wants' to
> connect the dots by dipping below zero. Is there a way to force
> positivity? Or should I simply adjust the InterpolationOrder until it
> works?
>
> Thanks,
> Dan
>
dantimatter
dan
Noqsi
One way is to transform the data. Choose a positive real function f,
apply its inverse to the data. Tranform the interpolation back using
f.
Exp[] often is a good choice for f. Using Valeri's example data:
data = {{0, 2}, {1, 1}, {2, .2}, {3, .01}, {4, 1}, {5, 0.9}, {6,
0.05}, {7, 0.1}, {8, 1}}
logdata = data /. {x_, y_} -> {x, Log[y]}
interp[x_] = Exp[Interpolation[logdata][x]]
If you plot interp[x], you'll see that it avoids zero.
dh
Sorry there is a type. Should read InterpolationOrder->1.
Daniel
DrMajorBob
peaks, reaching far above any of the data:
data = Sort@RandomReal[{0.1, 1}, {20, 2}];
{min, max} = data[[{1, -1}, 1]]
f = Interpolation@data;
logdata = data /. {x_, y_} -> {x, Log[y]};
interp = Exp[Interpolation[logdata]@#] &;
Show[Plot[Evaluate@Through[{f, interp}@x], {x, min, max},
PlotRange -> All]]
(run it several times)
The same thing happens here, but not as dramatically:
data = Sort@RandomReal[{0.1, 1}, {20, 2}];
{min, max} = data[[{1, -1}, 1]]
f = Interpolation@data;
logdata = data /. {x_, y_} -> {x, Sqrt[y]};
interp = (Interpolation[logdata]@#)^2 &;
Show[Plot[Evaluate@Through[{f, interp}@x], {x, min, max},
PlotRange -> All]]
Bobby
On Wed, 13 Jan 2010 04:58:01 -0600, Noqsi <j...@noqsi.com> wrote:
> On Jan 11, 4:52 pm, dantimatter <goo...@dantimatter.com> wrote:
> One way is to transform the data. Choose a positive real function f,
> apply its inverse to the data. Tranform the interpolation back using
> f.
>
> Exp[] often is a good choice for f. Using Valeri's example data:
>
> data = {{0, 2}, {1, 1}, {2, .2}, {3, .01}, {4, 1}, {5, 0.9}, {6,
> 0.05}, {7, 0.1}, {8, 1}}
> logdata = data /. {x_, y_} -> {x, Log[y]}
> interp[x_] = Exp[Interpolation[logdata][x]]
>
> If you plot interp[x], you'll see that it avoids zero.
>
schochet123
> For some data, that works pretty well; for other samples it has HUGE
> peaks, reaching far above any of the data:
This problem can be avoided by using a transformation that is close to
the identity for positive values:
trans[x_]=(Sqrt[1 + x^2] + x)/2
inv[y_]=(4 y^2 - 1)/(4 y)
Simplify[{trans[inv[y]], inv[trans[x]]}, y > 0]
To compare the effects of this new transformation with the old one,
modify the test code slightly:
tryinterp := (data = Sort@RandomReal[{0.1, 1}, {20, 2}];
{min, max} = data[[{1, -1}, 1]];
f = Interpolation@data;
logdata = data /. {x_, y_} -> {x, Log[y]};
invdata = data /. {x_, y_} -> {x, inv[y]};
expinterp = Exp[Interpolation[logdata]@#] &;
transinterp = trans[Interpolation[data]@#] &;
plota =
Plot[Evaluate@Through[{f, expinterp}@x], {x, min, max},
PlotRange -> All];
plotb = Plot[Evaluate@Through[{f, transinterp}@x], {x, min, max},
PlotRange -> All]; GraphicsRow[{plota, plotb}])
Now evaluate tryinterp several times
Steve
DrMajorBob
Bobby
On Fri, 15 Jan 2010 02:18:54 -0600, schochet123 <schoc...@gmail.com>
wrote:
Ray Koopman
> On Jan 14, 12:46 pm, DrMajorBob <btre...@austin.rr.com> wrote:
>
>> For some data, that works pretty well; for other samples
>> it has HUGE peaks, reaching far above any of the data:
>
> This problem can be avoided by using a transformation that
> is close to the identity for positive values:
>
>
> inv[y_] = (4 y^2 - 1)/(4 y)
There's a whole family of transformations here:
trans[c_,x_] = (Sqrt[4c + x^2) + x)/2
inv[c_,y_] = y - c/y
Extremely small values of c give results that approach those obtained
by chopping the interpolation at 0.
Noqsi
These transformations have the property that the interpolation behaves
differently for small and large numbers, with the transition at a
scale given by Sqrt[c]. The Log/Exp pair I suggested is scaleless.
{Sqrt[x],y^2} is another scaleless transformation.
All such transformations will exhibit artifacts from some point of
view given some data set: there is no flawless universal interpolation
procedure.
The {Sqrt[x],y^2} pair is useful for data that's proportional to
counts of independent events (like electric current over a potential
barrier), as it approximately gives each point its proper "weight" in
the analysis, given the correlation of the shot noise variance with
the expectation. However, the interpolation may "bounce" off the x-
axis...