MCMC и байесова статистика в BASIC

Original author: Rasmus Bååth
  • Translation
BASIC был одним из самых распространенных языков программирования. В 80-х он шел в стандартном наборе программ на компьютере (например, Commodore 64 и Apple II), а в 90х и DOS и Windows 95 включали в себя QBasic IDE.

QBasic был также моим первым языком программирования. Я не программировал на Бейсике уже почти 20 лет и решил вспомнить этот действительно странный язык. Поскольку я провел много времени за байесовскими алгоритмами, я подумал, что будет интересно увидеть как байесовская аналитика будет выглядеть в утилите 20-летней давности.

image


Здесь мы реализуем алгоритм Метрополиса-Гастингса, стандартный метод Монте-Карло по схеме марковских цепей (MCMC), который используется в байесовских моделях, в Бейсике.

Я применил распределение Лапласа к самому очаровательному набору данных: количеству волчат на одну стаю. Выборка состояла из 16 стай. В итоге я подсчитал и вывел результат, по-прежнему используя Бейсик (источник).

image

Где достать BASIC


Есть много различных версий BASIC, но я взял тот, на котором вырос: Microsoft QBasic 1.1. QBasic имеет много расширенных функций таких как типы, определяемые пользователем и (внимание!) functions. Но я не использовал никаких дополнительных преференций, я ведь собирался писать олдскульный BASIC, используя строчную нумерацию и GOTO, что означает, что код можно легко импортировать в, скажем, Commodore 64 BASIC.

Скачать QBasic можно здесь. К тому же это бесплатно. Если конечно вы не используете DOS до сих пор, то следующим шагом будет установка эмулятора DOSBox. После того, как вы запустите QBASIC.EXE вашему взору предстанет очень дружественный, ярко-синий интерфейс. Его можно потестить с помощью кастомного скрипта. Он очистит экран и напечатает “HELLO WORLD”.

image

Обратите внимание, что в QBasic нумерация строк не так важна, но в более старом BASIC (например в том, что стоял на Commodore 64) она необходима, поскольку программа выполняется по порядку.

Что мы будем реализовывать


Мы будем реализовывать алгоритм Метрополиса-Гастингса на классической MCMC, с которым вы столкнетесь в первую очередь, если будете изучать вычислительные методы для байесовских моделей.

Байесовская модель, которую мы будем представлять, имеет простое одномерное распределение Лапласа (просто чтобы дать гауссовскому распределению отдохнуть). Распределение Лапласа, как и Гауссовское, конечно и симметрично, но имеет более острый пик и пологие хвосты.

Оно обладает двумя параметрами: параметр сдвига μ, который определяет также его среднее значение и медиану и размах b, который определяет ширину распределения.

image

В Лаплассовском распределении медиана является максимальной вероятностной оценкой параметра μ. Чтобы включить это в простейшую байесовскую модель мы должны иметь два параметрa. Здесь я буду упрощенно применять Uniform(−∞,∞) к μ и log(b), это, P(μ),P(log(b))∝1. Полная модель:

x∼Laplace(μ,b)
μ∼Uniform(−∞,∞)
log(b)∼Uniform(−∞,∞)

Повторюсь, что данные, которые мы будем использовать, возможно самые милые, с которыми я работал. Они содержат количество волчат в выборке из 16 волчьих стай(source):

image

Реализация в R


Прежде чем углубиться в BASIC, реализуем все в R:

# Набор данных о волчатах
x <- c(5, 8, 7, 5, 3, 4, 3, 9, 5, 8, 5, 6, 5, 6, 4, 7)

# The log posterior density of the Laplace distribution model, when assuming 
# uniform/flat priors. The Laplace distribution is not part of base R but is
# available in the VGAM package.
model <- function(pars) {
  sum(VGAM::dlaplace(x, pars[1], exp(pars[2]), log = TRUE))
}

# The Metropolis-Hastings algorithm using a Uniform(-0.5, 0.5) proposal distribution 
metrop <- function(n_samples, model, inits) {
  samples <- matrix(NA, nrow = n_samples, ncol = length(inits))
  samples[1,] <- inits
  for(i in 2:n_samples) {
    curr_log_dens <- model(samples[i - 1, ])
    proposal <- samples[i - 1, ] + runif(length(inits), -0.5, 0.5)
    proposal_log_dens <- model(proposal)
    if(runif(1) < exp(proposal_log_dens - curr_log_dens)) {
      samples[i, ] <- proposal
    } else {
      samples[i, ] <- samples[i - 1, ]
    }
  }
  samples
}

samples <- metrop(n_samples = 1000, model, inits = c(0,0))
# Plotting a traceplot
plot(samples[,1], type = "l", ylab = expression(Location ~ mu), col = "blue")
# Calculating median posterior and 95% CI discarding the first 250 draws as "burnin".
quantile(samples[250:1000,1], c(0.025, 0.5, 0.975))


image

##  2.5%   50% 97.5% 
## 4.489 5.184 6.144

(Этот скрипт просто показывает, что должно получиться на Бейсике. Если вы хотите попробовать метода Метрополиса-Гастингса на MCMC в R, воспользуйтесь MCMCmetrop1R в MCMCpack package, или этим Metropolis-Hastings script.

Модель в BASIC


Начнем с очистки экрана (CLS
) и определения переменных. DIM
определяет массивы и матрицы.

image

Даже в этом коротком куске кода можно обнаружить целую кучу странностей:
  • Поскольку мы пытаемся придерживаться олдскула, то везде стоит нумерация строк. Она не добавляется автоматически, а прописывается мной и иногда я пропускаю строчки. Отсюда хороший лайф-хак — каждый раз добавлять к строчке по 10, чтобы можно было вставить номер в пропущенную строку и не переписывать остальные (например 75
    на скриншоте). Поскольку в старом Бейсике нет функций, нумерация строк будет использоваться командами GOSUB
    и GOTO.
  • Все операторы прописываются в верхнем регистре и если написать оператор в нижнем регистре, он исправится на верхний.
  • Почему SAMPLES!
    В QBasic все переменные являются числами и чтобы изменить их тип нужно добавить символ к имени переменной. Так что THIS$
    это строка, а THAT!
    — это число с плавающей точкой. Поскольку параметры являются непрерывными, то большинство переменных будут заканчиваться на "!".
  • Команда DATA
    — это неплохой способ поместить данные в программу, позже их можно извлечь с помощью функции READ.
    Если у вас более одного набора данных, желаю удачи.

Продолжим определение модели. Мы представим ее в виде подпрограммы, которая заканчивается выражением RETURN
. Она вызывается выражением GOSUB
,
которое переходит на , выполняет код до выражения RETURN
и затем переходит на строку, которая идет после GOSUB.
Это похоже на недо-функцию, которая не может передать аргументы или вернуть значения. Но все хорошо, пока все переменные глобальны. Подпрограмма может обращаться и устанавливать любые переменные, пока вы не используете одинаковые имена переменных для двух разных вещей. Подпрограмма ниже предполагает, что вектор PARAMS!
перезаписывает переменную LOGDENS!.
Подпрограмма также использует READ X!,
чтобы считывать значения из DATA
и RESTORE
чтобы сбросить READ.


image

А вот и алгоритм Метрополиса-Гастингса. Ничего нового, кроме GOSUB 520,
который используется чтобы перейти в подпрограмму модели на строку 520. Далее мы используем RND,
чтобы равномерно сгенерировать рандомные числа в промежутке от 0 до 1.

image

Подсчет и вывод на экран


Теперь можно собрать все вместе с помощью следующего кода:

image

Код ниже берет матрицу SAMPLES!
и рассчитывает среднее значение, стандартное отклонение и 95 интервал доверия.

image

Еще немного странностей от Бейсика: PRINT "HELLO", "YOU!" выведет на экран HELLO YOU!, а PRINT "HELLO"; "YOU!" выведет HELLOYOU! Разделитель "," или ";" определяет будет ли пробел между выводимыми объектами.

И наконец мы хотим вывести график, чтобы удостоверится, что цепь Маркова выглядит нормально. И хоть в Бейсике нет даже встроенной функции для сортировки, зато есть колоссальная поддержка графики. Кусок кода ниже просто масштабирует координаты выборки, чтобы они подходили под разрешение 320 x 200. Обратите внимание на забавный синтаксис для отрисовки линий: LINE (, )-(, ).


image

Мы наверное захотим записать реализацию на диск:

image

Выполнение


Теперь наконец мы готовы собрать все воедино, сформировать и вывести.

image

Самый нижний GOTO
просто переходит на последнюю строку программы, чтобы выйти из нее. Если мы установим эмуляцию скорости 33 MHz, как в 386 компе (согласно этой странице), то генерация выборки, подсчет и вывод займут около минуты. Вот результат:

image
image

Так что там с волчатами


Этот пример был не совсем об анализе данных, но глядя на результаты можно предположить, что среднее количество щенков в стае колеблется от 4-5 до 6. Я должен также отметить, что обе реализации на BASIC и на R похожи:

image

Здесь полная программа в текстовом формате (MCMCDEMO.BAS),

Конспект


  1. Мы можем работать с довольно продвинутой байесовой статистикой с помощью утилиты, которую использовали 20 лет назад.
  2. Программирование с помощью GOTO
    и GOSUB
    очень-очень олдскульное.
  3. Очень приятно не писать множество скобок (Но с другой строны приходится писать буквенные END IF
    и NEXT I.
    )
  4. Запуск программы на 386 компьютере прошел довольно быстро.
  5. Очень тяжело писать код, когда все переменные глобальны.
.io
Company
AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 20

    –1
    Когда мне пришлось в последний раз столкнуться с написанием чего-то c помощью QBASIC, я немного заморочился и при команде «Build» из Sublime Text запускался Dosbox с бэйсиком, в котором и запускалась программа. Лучше бы в качестве первого языка в некоторых учебных заведениях преподавали Python.
      +3
      Первый язык должен быть строго типизированным
        +4
        Еще один неотличатор строгой типизации от слабой и динамической от статической?)
          +3
          ОК. Первый язык должен быть строго и статически типизирован. Чем за большее число ошибок даёт по рукам с внятным объяснением компилятор, тем лучше.
            0
            С таким подходом до звания программист будут доходить одни только упоротые зануды.
            /me упоротый зануда
              +1
              Лучше, когда компилятор подробно и обстоятельно объясняет человеку, почему тот дурак, чем вылезает непонятная ошибка в рантайме.
                0
                Только не первый язык. Задача первого — дать понимание в общих чертах, здесь битье по рукам не нужно. А вот первый серьезный язык — другое дело.
                Не работает этот принцип «с азов к сложному поступательно» при обучении программированию. Сначала надо влезть на «среднюю» ступеньку, что-то сделать реально работающее, а уже потом нырять вниз, к основам, и переходить к серьезным языкам.
                +1
                >С таким подходом до звания программист будут доходить одни только упоротые зануды.
                Мир станет лучше.
                  0
                  Мир не станет лучше, если все программисты будут упоротыми занудами.
            0
            Из‐за стандартных типов и общего принципа «explicit is better then implicit» Python считается строго типизированным. Впрочем, свои классы вы можете писать так, что ваш код строго типизированным не будет.
              0
              И первый учитель программирования не должен знать что вообще происходит в программе и читать урок с книжки…
                0
                Даже если первый язык должен быть строго типизирован, то всё равно лучше Python, чем BASIC в 2015 году.
                  0
                  Не, в первом языке невидимые символы не должны влиять на семантику.
              0
              Зачем нумеровать строки? QBASIC поддерживает метки.
              Похоже, Вы из любви к искусству выбрали стиль совместимый с более старыми BASIC?
                0
                В статье так и написано, что автор хотел писать на старом BASIC.
                  0
                  А зачем тогда насиловать QBASIC? Взять какой-нибудь GWBASIC 1.0. Его ничуть не сложнее найти, чем QBASIC…
                +4
                IF...ENDIF — не олдскульные блоки. Олдскульный бейсик поддерживал только однострочный IF… THEN 150 ELSE 250 (GOTO можно было опустить). Или, например IF I>NBURNIN THEN COL=1.
                А блоки — это как раз синтаксис Q/Quick/Visual basic.
                  0
                  яростно плюсую, зачем лезть куда-то не разобравшись, да еще и делать громкие неверные заявления.
                  Коммент ниже про суффиксы тоже верный, все числа по умолчанию вещественные.

                  Что забавного в синтаксисе line я так и не понял. Странно, что автор не упомянул как раз забавную особенность, что последняя точка сохраняется и первый аргумент (пару x1,y1) можно не указывать.

                  Как раз на qbasic можно писать через sub и без глобальных переменных, вызов меню управления sub'ами по F2, если не ошибаюсь.

                  А главное — зачем писать на qb в стиле старого basic? возьмите gwbasic, он вроде прекрасно работает под виндой (пробовал на xp вроде)
                  –6
                  По прочтению только один вопрос — зачем?
                  Времени вагон и девать его некуда?
                  Ведь можно изучить что-то ещё более полезное и более продвинутое.
                    +3
                    Да вроде же там переменные вещественные по умолчанию. Суффиксы надо указывать только для строк ($) или для целых (%). Хотя, может это не тот диалект.

                    Для полной олдскульности попробуйте обойтись именами переменных из одной буквы или одной буквы и цифры (вида M или L1).

                    Only users with full accounts can post comments. Log in, please.