Подгонка кривой: найдите самую гладкую функцию, которая удовлетворяет списку ограничений - PullRequest
7 голосов
/ 24 апреля 2010

Рассмотрим множество неубывающих сюръективных (на) функций от (-inf, inf) до [0,1]. (Типичные CDF s удовлетворяют этому свойству.) Другими словами, для любого действительного числа x 0 логистическая функция , пожалуй, самый известный пример.

Теперь нам даны некоторые ограничения в виде списка значений x, а для каждого значения x - пара значений y, между которыми должна находиться функция. Мы можем представить это в виде списка {x, ymin, ymax} троек, таких как

constraints = {{0, 0, 0}, {1, 0.00311936, 0.00416369}, {2, 0.0847077, 0.109064}, 
 {3, 0.272142, 0.354692}, {4, 0.53198, 0.646113}, {5, 0.623413, 0.743102}, 
 {6, 0.744714, 0.905966}}

Графически это выглядит так:

ограничения на cdf http://yootles.com/outbox/cdffit1.png

Теперь мы ищем кривую, которая учитывает эти ограничения. Например:

встроенный cdf http://yootles.com/outbox/cdffit2.png

Давайте сначала попробуем простую интерполяцию через середины ограничений:

mids = ({#1, Mean[{#2,#3}]}&) @@@ constraints
f = Interpolation[mids, InterpolationOrder->0]

График, f выглядит так:

интерполированный cdf http://yootles.com/outbox/cdffit3.png

Эта функция не сюръективна. Кроме того, мы хотели бы, чтобы это было гладко. Мы можем увеличить порядок интерполяции, но теперь он нарушает ограничение, что его диапазон равен [0,1]:

интерполированный cdf с более высоким порядком интерполяции http://yootles.com/outbox/cdffit4.png

Таким образом, цель состоит в том, чтобы найти наиболее гладкую функцию , которая удовлетворяет ограничениям:

  1. неубывающая.
  2. Стремится к 0, когда x приближается к отрицательной бесконечности, и стремится к 1, когда x приближается к бесконечности.
  3. Проходит через данный список y-ошибок-баров.

Первый пример, который я привел выше, кажется хорошим кандидатом, но я сделал это с помощью функции FindFit Mathematica, предполагая lognormal CDF . Это хорошо работает в этом конкретном примере, но в общем случае необязательно иметь логарифмический CDF, который удовлетворяет ограничениям.

Ответы [ 3 ]

5 голосов
/ 24 апреля 2010

Не думаю, что вы указали достаточно критериев, чтобы сделать желаемый CDF уникальным.

Если единственные критерии, которые должны соблюдаться, это:

  1. CDF должен быть «довольно гладким» (см. Ниже)
  2. CDF должен быть неубывающим
  3. CDF должен проходить через y-интервалы "bar error"
  4. CDF должен стремиться к 0 при x -> -Infinity
  5. CDF должен стремиться к 1, поскольку x -> Infinity.

тогда, возможно, вы могли бы использовать Монотонная кубическая интерполяция . Это даст вам функцию C ^ 2 (дважды непрерывно дифференцируемую), которая, в отличие от кубических сплайнов, гарантированно будет монотонным, если даны монотонные данные.

Это оставляет открытым вопрос, какие именно данные вы должны использовать для генерации монотонная кубическая интерполяция. Если вы берете центральную точку (среднее значение) каждой ошибки бар, вы гарантированно, что результирующие точки данных являются монотонными повышение? Если нет, то вы можете сделать произвольный выбор, чтобы гарантировать что выбранные вами точки монотонно растут (поскольку критерии не заставляют наше решение быть уникальным).

Что теперь делать с последней точкой данных? Есть ли X, который гарантированно быть больше любого x в наборе данных ограничений? Возможно, вы можете снова сделать произвольный выбор удобства и выберите очень большой X и поместите (X, 1) в качестве конечная точка данных.

Комментарий 1: Ваша проблема может быть разбита на 2 подзадачи:

  1. Учитывая точные точки (x_i, y_i), через которые должен пройти CDF, как вы генерируете CDF? Я подозреваю, что существует бесконечно много возможных решений, даже с ограничением бесконечной гладкости.

  2. Как выбрать y (x_i, y_i) для y-ошибок? Опять же, существует бесконечно много возможных решений. Возможно, потребуется добавить некоторые дополнительные критерии, чтобы сделать уникальный выбор. Дополнительные критерии также могут сделать проблему еще сложнее, чем в настоящее время.

Комментарий 2: Вот способ использовать монотонную кубическую интерполяцию и удовлетворять критериям 4 и 5:

Монотонная кубическая интерполяция (назовем это f) отображает R -> R .

Пусть CDF(x) = exp(-exp(f(x))). Тогда CDF: R --> (0,1). Если бы мы могли найти соответствующее f, то, определив CDF таким образом, мы могли бы удовлетворить критерии 4 и 5.

Чтобы найти f, преобразуйте ограничения CDF (x_0,y_0),...,(x_n,y_n), используя преобразование xhat_i = x_i, yhat_i = log(-log(y_i)). Это обратное преобразование CDF. Если y_i увеличиваются, то yhat_i уменьшаются.

Теперь примените монотонную кубическую интерполяцию к точкам данных (x_hat, y_hat) для генерации f. Затем, наконец, определите CDF(x) = exp(-exp(f(x))). Это будет монотонно возрастающая функция от R -> (0,1), которая проходит через точки (x_i, y_i).

Это, я думаю, удовлетворяет всем критериям 2--5. Критерий 1 в некоторой степени удовлетворен, хотя, безусловно, могут существовать более плавные решения.

4 голосов
/ 25 апреля 2010

Я нашел решение, которое дает разумные результаты для различных входов. Я начинаю с подгонки модели - сначала к нижним пределам ограничений, а затем к верхним. Я буду называть среднее значение этих двух встроенных функций «идеальной функцией». Я использую эту идеальную функцию для экстраполяции слева и справа от того, где заканчиваются ограничения, а также для интерполяции между любыми пропусками в ограничениях. Я вычисляю значения для идеальной функции через равные промежутки времени, включая все ограничения, от того, где функция почти равна нулю слева, до почти одного справа. При ограничениях я обрезаю эти значения по мере необходимости, чтобы удовлетворить ограничениям. Наконец, я создаю интерполяционную функцию, которая проходит через эти значения.

Моя реализация Mathematica следует.
Сначала пара вспомогательных функций:

(* Distance from x to the nearest member of list l. *)
listdist[x_, l_List] := Min[Abs[x - #] & /@ l]

(* Return a value x for the variable var such that expr/.var->x is at least (or
   at most, if dir is -1) t. *)
invertish[expr_, var_, t_, dir_:1] := Module[{x = dir},
  While[dir*(expr /. var -> x) < dir*t, x *= 2];
  x]

А вот и основная функция:

(* Return a non-decreasing interpolating function that maps from the
   reals to [0,1] and that is as close as possible to expr[var] without
   violating the given constraints (a list of {x,ymin,ymax} triples).
   The model, expr, will have free parameters, params, so first do a
   model fit to choose the parameters to satisfy the constraints as well
   as possible. *)
cfit[constraints_, expr_, params_, var_] := 
Block[{xlist,bots,tops,loparams,hiparams,lofit,hifit,xmin,xmax,gap,aug,bests},
  xlist = First /@ constraints;
  bots = Most /@ constraints; (* bottom points of the constraints *)
  tops = constraints /. {x_, _, ymax_} -> {x, ymax};
  (* fit a model to the lower bounds of the constraints, and 
     to the upper bounds *)
  loparams = FindFit[bots, expr, params, var];
  hiparams = FindFit[tops, expr, params, var];
  lofit[z_] = (expr /. loparams /. var -> z);
  hifit[z_] = (expr /. hiparams /. var -> z);
  (* find x-values where the fitted function is very close to 0 and to 1 *)
  {xmin, xmax} = {
    Min@Append[xlist, invertish[expr /. hiparams, var, 10^-6, -1]],
    Max@Append[xlist, invertish[expr /. loparams, var, 1-10^-6]]};
  (* the smallest gap between x-values in constraints *)
  gap = Min[(#2 - #1 &) @@@ Partition[Sort[xlist], 2, 1]];
  (* augment the constraints to fill in any gaps and extrapolate so there are 
     constraints everywhere from where the function is almost 0 to where it's 
     almost 1 *)
  aug = SortBy[Join[constraints, Select[Table[{x, lofit[x], hifit[x]}, 
                                              {x, xmin,xmax, gap}], 
                                        listdist[#[[1]],xlist]>gap&]], First];
  (* pick a y-value from each constraint that is as close as possible to 
     the mean of lofit and hifit *)
  bests = ({#1, Clip[(lofit[#1] + hifit[#1])/2, {#2, #3}]} &) @@@ aug;
  Interpolation[bests, InterpolationOrder -> 3]]

Например, мы можем соответствовать логнормальной, нормальной или логистической функции:

g1 = cfit[constraints, CDF[LogNormalDistribution[mu,sigma], z], {mu,sigma}, z]
g2 = cfit[constraints, CDF[NormalDistribution[mu,sigma], z], {mu,sigma}, z]
g3 = cfit[constraints, 1/(1 + c*Exp[-k*z]), {c,k}, z]

Вот как это выглядит для моего первоначального списка примеров ограничений:

ограниченное соответствие логнормальной, нормальной и логистической функции http://yootles.com/outbox/cdffit5.png

Норма и логистика почти друг на друга, а логнормал - синяя кривая.

Это не совсем идеально. В частности, они не совсем однообразные. Вот сюжет производных:

Plot[{g1'[x], g2'[x], g3'[x]}, {x, 0, 10}]

производные встроенных функций http://yootles.com/outbox/cdffit6.png

Это показывает некоторое отсутствие гладкости, а также небольшую немонотонность вблизи нуля. Я приветствую улучшения этого решения!

0 голосов
/ 24 апреля 2010

Вы можете попытаться провести кривую Безье через средние точки. В частности, я думаю, что вы хотите C2 непрерывной кривой.

...