В июне 2014 года, как это обычно бывает по чётным месяцам, был проведён конкурс по функциональному программированию, который проводится под эгидой Фонда Поддержки Функционального Программирования ФП(ФП). Традиционно я хотел бы подвести итоги конкурса и рассказать о решении конкурсной задачи при помощи языка программирования Haskell. Так что всех заинтересованных я приглашаю ознакомиться с этой небольшой заметкой.
В качестве задачи на конкурс была предложена задача по поиску закономерностей в ряду проявлений казалось бы «случайного» события. Но как и всё в этом мире чисто случайными являются, видимо, результаты измерения квантовых состояний, так что во всём другом можно найти какие-то закономерности. Так и здесь. Был дан список дат, когда произошло некоторое событие, и предлагалось дать ответы на два вопроса:
- Каков минимальный период, в котором частотная вероятность проявления события хотя бы в один день периода равна или более 50 %?
- Необходимо было дать прогноз проявления события с даты конкурса до конца текущего года.
Только два конкурсанта смогли предоставить решения. Впрочем, оба они были неправильными, поскольку правильным ответом на первый вопрос является число 24. А вот вторая задача будет обработана в конце года, когда будет явлена статистика по проявлениям событий. Так что приз за первый вопрос остался неразыгранным, а приз за второй вопрос будет предоставлен тому конкурсанту, прогноз которого наберёт больше очков, в следующем году.
Ну а здесь остаётся рассмотреть решение этих задач на языке программирования Haskell.
Поиск периода
Прежде всего был определён список проявлений события:
dates :: [String]
dates = ["27.09.2013", "06.10.2013", "23.10.2013", "06.11.2013", "26.11.2013",
"27.11.2013", "21.12.2013", "30.12.2013", "06.01.2014", "16.01.2014",
"17.01.2014", "21.01.2014", "25.01.2014", "26.01.2014", "05.02.2014",
"11.02.2014", "21.02.2014", "02.03.2014", "07.03.2014", "30.03.2014",
"08.04.2014", "18.04.2014", "23.04.2014", "27.04.2014", "02.05.2014",
"15.05.2014", "17.05.2014", "18.05.2014", "19.05.2014", "20.05.2014",
"25.05.2014", "26.05.2014", "28.05.2014"]
Даже беглый взгляд покажет, что в этом списке нет явных закономерностей. Но вот медитирование над ним несколько дней кряду может выявить что-то. Чтобы не предаваться бессмысленной трате времени на медитации, как раз и решено было написать программу для поиска закономерностей в периодах различной длины. Начнём писать эту программу, как это обычно бывает в разработке на функциональных языках программирования, сверху вниз. Вот главная функция программы:
main :: IO ()
main = do putStrLn ("В последовательности дат обнаружены " ++
"закономерности в минимальном периоде длиной " ++
show (length findMinimalPeriod) ++ " дней.")
revealSequences findMinimalPeriod
Здесь есть вызов двух функций, которые мы рассмотрим далее. Первая, findMinimalPeriod, возвращает период минимальной длины для заданного порога вероятности. Внимательный читатель заметит, что ей не передаётся никаких аргументов, поэтому порог вероятности определён где-то в виде константы. В общем-то, это порочная практика для парадигмы функционального программирования, но иногда в исследовательских целях это помогает быстро достичь результатов. Так что вот определение:
significance :: Int
significance = 5
lowProbability :: Float
lowProbability = 0.5
findMinimalPeriod :: [(Int, Float)]
findMinimalPeriod = head $
filter (l -> maximum (map snd l) >= lowProbability) $
map process [1 .. interval `div` significance]
Константа significance определяет минимальную «высоту» цилиндра, на который наматывается шкала времени (ведь для того, чтобы найти периоды, можно отметить даты проявления события на длинной ленте, которую потом намотать спиралью на цилиндр с заданной длиной окружности, которая определяет период; соответственно закономерности будут выглядеть как вертикальные линии). Ну а константа lowProbability задаёт минимальный порог вероятности проявления события. Сама же функция findMinimalPeriod берёт голову списка, полученного после фильтрации списка на наличие вероятности не менее заданного порога, который (список) получен при помощи обработки (функция process) чисел от 1 до некоторой верхней границы.
Верхняя граница определяется при помощи функции interval, определение которой следующее:
interval :: Int
interval = finish - start
where
start = stringToDayOfYear $ head dates
finish = 365 + stringToDayOfYear (last dates)
Как видно, здесь мы считаем длину интервала, в котором заданы даты проявлений события. Вычитаем из последней даты первую (возможно, что надо было бы прибавить 1). Эта функция не очень хороша, поскольку в ней есть число 365, а это значит, что она не универсальна. Ну да ладно. К тому же, предыдущая функция (findMinimalPeriod) вообще написана из рук вон плохо и может выбросить ошибку времени исполнения из-за отсутствия проверки на пустоту списка, передаваемого в функцию head.
Теперь перейдём к определению функции process:
process :: Fractional a => Int -> [(Int, a)]
process p = map (l -> (fst $ head l, ((/) `on` fromInteger) (sum $ map snd l) (toEnum $ length l))) $
groupBy ((==) `on` fst) $
sortBy (comparing fst) $
map (first (`mod` p) . i -> if i `elem` ds3
then (i, 1)
else (i, 0)) [head ds3 .. last ds3]
where
ds1 = map stringToDayOfYear dates
ds2 = uncurry (++) $ second (map (+ 365)) $ span (>= head ds1) ds1
ds3 = map (subtract (head ds2)) ds2
Функция получает на вход длину периода, а возвращает список пар (гистограмму), в которых первым элементом является номер дня в периоде, а вторым — частотная вероятность проявления события в этот день. При помощи локальных определений ds1, ds2 и ds3 строится список последовательных номеров дней проявления события, начиная от первого дня в списке dates. Далее этот список подвергается такой процедуре. Для всех номеров от 1 до номера последней даты проявления события ищется остаток от деления на длину периода. Для всех таких остатков ставится флаг 0, если в соответствующий день проявления события не было, и флаг 1 — если было. Затем список остатков с флагами группируется по остаткам, после чего группы схлопываются в пары вида (номер дня в периоде, вероятность проявления события). Всё.
Тут надо рассмотреть ещё две сервисные функции:
double :: a -> (a, a)
double x = (x, x)
stringToDayOfYear :: String -> Int
stringToDayOfYear = uncurry (monthAndDayToDayOfYear False) .
(read . take 2 . drop 3 *** read . take 2) .
double
Про первую и говорить нечего (странно только то, что её определения нет в стандартном модуле Prelude; хотя это и понятно, поскольку оно настолько тривиально). Вторая функция переводит дату из строкового представления "DD.MM.YYYY" в числовое, принятое в модуле Data.Time.Calendar.MonthDay, при помощи которого обрабатываются даты.
Наконец, определим функцию revealSequences:
revealSequences :: [(Int, Float)] -> IO ()
revealSequences ps = do let l = length ps
(d1, p1) = maximumBy (comparing snd) ps
(d2, p2) = maximumBy (comparing snd) $ delete (d1, p1) ps
putStrLn ("Максимальное проявление события (вероятность: " ++
show p1 ++ ") происходит на " ++ show (d1 + 1) ++
"-й день " ++ show l ++ "-дневного периода.")
putStrLn ("Второй максимум (вероятность: " ++ show p2 ++
") происходит на " ++ show (d2 + 1) ++ "-й день.")
Как видно, она просто находит и выводит на экран два пика проявления события в заданном списке (гистограмме). Эта функция реализована исключительно для удобства исследования и уже не имеет отношения к решению поставленной в конкурсе задачи.
Прогнозирование
Теперь перейдём ко второму вопросу конкурса — прогнозированию дат проявления события до конца года. Надо отметить, что вообще задача прогнозирования является делом неблагодарным. Один вопрос, когда у нас есть чёткая закономерность (хотя бы даже и вероятностная), и можно более или менее применить детерминированную формулу. Совсем другое дело, когда есть некоторое (совсем небольшое) количество значений, на основании которых надо построить возможное будущее. Точек бифуркаций мириады, так что даже допуски и интервальная арифметика не помогут. Тем не менее, можно придумать какой-нибудь метод и попробовать оценить его применимость уже задним умом, когда пройдёт определённый период времени, и мы сможем сравнить навангованные значения с фактически произошедшими.
Поступим здесь именно таким образом. В качестве метода можно предложить такой. Рассмотрим все возможные периоды, начиная от минимального найденного на предыдущем этапе (меньшие него, честно говоря, не вызывают доверия с точки зрения статистической правдоподобности) и заканчивая длиной известного периода наблюдений, делённой на норму допуска («толщину», как мы её назвали ранее). Для каждого из таких периодов получим вероятности проявления событий в каждый день периода, как мы уже научились делать это при помощи функции process. Далее просто найдём среднюю вероятность для каждого дня по всем вероятностям, найденным для каждого рассматриваемого периода.
Именно этот метод реализует следующая функция:
forecast :: FilePath -> String -> String -> IO ()
forecast fp sd fd = do let (b, e) = (length findMinimalPeriod, interval `div` significance)
writeFile fp $
unlines $
map (((n, q) -> let (m, d) = dayOfYearToMonthAndDay False (n - 365)
in prettyShowInt d ++ "." ++
prettyShowInt m ++ ".2014: " ++
prettyShowFloat q) . second (/ toEnum (e - b + 1))) $
foldr1 (zipWith ((d, q1) (_, q2) -> (d, q1 + q2))) $
map getProbabilities [b..e]
where
getProbabilities p = let ds = stringToDayOfYear $ head dates
fs = 365 + stringToDayOfYear (last dates)
d1 = 365 + stringToDayOfYear sd
d2 = 365 + stringToDayOfYear fd
in drop (interval + (d1 - fs)) $
zipWith (x (_, q) -> (x, q)) [ds..d2] $
cycle $
process p
leadingZero :: String -> String
leadingZero [c] = '0' : [c]
leadingZero c = c
prettyShowInt i = leadingZero $ show i
prettyShowFloat f = let (d, r) = span (/= '.') $ show (f * 100)
in leadingZero d ++ take 5 (r ++ cycle "0")
Её определение выглядит несколько монструозно, но ядром вычислений является локальное определение getProbabilities (из названия должно быть понятно, какому шагу метода оно соответствует). Остальное — всего лишь обвязка для вывода полученных значений в файл в заданном условиями конкурса формате.
В общем-то, всё. Теперь осталось только дождаться конца года и сравнить прогноз с фактом.
Заключение
Само собой разумеется, что описанный метод прогнозирования может быть проверен только на практике. В принципе, прогнозирование в таких задачах является делом неблагодарным, поэтому здесь либо не пытаются заниматься прогнозами, либо используют танцы с бубнами и прочие эвристики. Собственно, как это было сделано в представленном решении. Ну и будущее покажет.
Как обычно исходный код представленного решения всякий желающий может получить здесь.
Мои предыдущие статьи о конкурсах по ФП на Хаброхабре:
- Конкурсы в 2011 году: Альманах конкурсов по ФП за 2011 год
- Конкурсы в 2012 году: Альманах конкурсов по ФП за 2012 год
- Конкурсы в 2013 году: Альманах конкурсов по ФП за 2013 год
- Задача о выкидывании негра за борт
- Апрельский конкурс по ФП: поэтический вечер
Автор: Darkus
