Подбор кривых: Найдите самую гладкую функцию, которая удовлетворяет список ограничений

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

Если единственные критерии, которые должны выполняться:

1. CDF должен быть "достаточно гладким" (см. Ниже)
2. CDF должен быть неубывающим
3. CDF должен пройти через "полосу ошибок" y-интервалы
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)? Опять же, возможных решений бесконечно много. Может потребоваться добавить некоторые дополнительные критерии, чтобы сделать уникальный выбор. Дополнительные критерии также, вероятно, усложнили бы проблему, чем она есть сейчас.

Комментарий 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 в некоторой степени удовлетворяется, хотя, безусловно, могут быть более плавные решения.

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

Моя реализация 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]

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

_{(источник: yootles.com )}

Нормальный и логистический почти накладываются друг на друга, а логнормальный - это синяя кривая.

Они не совсем идеальны. В частности, они не совсем монотонны. Вот график производных:

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

_{(источник: yootles.com ])}

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

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