Свернуть работу с интерполированными функциями в Mathematica - PullRequest
Новая
       40

Свернуть работу с интерполированными функциями в Mathematica

3 голосов
/ 25 сентября 2010

Я использую Mathematica 7.

У меня есть интерполированная функция, вот пример: 

pressures = 
  WeatherData["Chicago", "Pressure", {2010, 8}] // 
     DeleteCases[#, {_, _Missing}] & // 
    Map[{AbsoluteTime[#[[1]]], #[[2]]} &, #] & // Interpolation;

Я бы хотел вычислить его производную, что прямо вперед: 

dpressures = D[pressures[x], x]

Теперь, если вы строите эту функцию 

Plot[3600*dpressures, {x, AbsoluteTime[{2010, 8, 2}], AbsoluteTime[{2010, 8, 30}]}]

(извините, не знаете, как разместить изображение из Mathematica, и у вас нет времени, чтобы выяснить это.) Вы обнаружите, что очень шумно.Итак, я хотел бы сгладить это.Моей первой мыслью было использовать Convolve и интегрировать его с ядром Гаусса, что-то вроде следующего: 

a = Convolve[PDF[NormalDistribution[0, 5], x], 3600*dpressures, x, y]

Возвращает 

360 Sqrt[2/\[Pi]] Convolve[E^(-(x^2/50)), InterpolatingFunction[{{3.48961266 10^9, 3.49228746 10^9}},<>], ][x], x, y]

Что мне кажется разумнымК сожалению, я полагаю, что где-то допустил ошибку, потому что результат, который я получаю, кажется не поддающимся оценке.То есть: 

a /. y -> AbsoluteTime[{2010, 8, 2}]

Возвращает 

360 Sqrt[2/\[Pi]] Convolve[E^(-(x^2/50)), InterpolatingFunction[{{3.48961266 10^9, 3.49228746 10^9}},<>][x], x, 3489696000]]

Что просто не то, что я искал, я ожидаю число от -1 до 1.

1 Ответ

5 голосов
/ 26 сентября 2010

Convolve ищет закрытую форму для свертки.Вы можете попробовать числовую свертку, начиная с чего-то вроде 

NConvolve[f_, g_, x_, y_?NumericQ] := 
 NIntegrate[f (g /. x -> y - x), {x, -Infinity, Infinity}]

Однако для этой шумной негладкой функции численное интегрирование будет бороться.(Он не будет работать с настройками по умолчанию и будет медленным даже при тщательно выбранных настройках.)

Я предлагаю вам работать непосредственно с базовыми данными вместо интерполяции зашумленных данных.

границы вашего временного диапазона: 

In[89]:= {lower = Min[First[pressures]], upper = Max[First[pressures]]}    
Out[89]= {3.48961*10^9, 3.49229*10^9}

Используйте интерполяцию для получения образцов каждый час *: 

data = Table[pressures[x], {x, lower, upper, 3600}];

Теперь сравните 

ListLinePlot[Differences[data]]

со сглаженной версией более 5часовое окно: 

ListLinePlot[GaussianFilter[Differences[data], 5]]
  • Вы можете использовать InterpolationOrder -> 1 для шумных данных.
...