2010-09-03 3 views
2

Я пишу модель Проверку, которая опирается на вычисление коэффициента, который используется интенсивно алгоритмов, заключается в следующем:Как оптимизировать это вычисление

[альтернативный текст] [1]

где q is double, t a double and k a int. e обозначает экспоненциальную функцию. Этот коэффициент используется на этапах, в которых q и t не меняются, а k всегда начинается с 0 до суммы всех предыдущих коэффициентов (этого шага) достигает 1.

Моя первая реализация была буквальной:

let rec fact k = 
    match k with 
    0 | 1 -> 1 
    | n -> n * (fact (k - 1)) 

let coeff q t k = exp(-. q *. t) *. ((q *. t) ** (float k)) /. float (fact k) 

Конечно, это не так долго, поскольку вычисление всего факториала было просто неосуществимым, когда k преодолел небольшой порог (15-20): очевидно, результаты начали сходить с ума. Поэтому я перестроил все это делая дополнительных подразделений:

let rec div_by_fact v d = 
    match d with 
    1. | 0. -> v 
    | d -> div_by_fact (v /. d) (d -. 1.) 

let coeff q t k = div_by_fact (exp(-. q *. t) *. ((q *. t) ** (float k))) (float k) 

Эта версия работает достаточно хорошо, когда q и t достаточно «нормальный», но когда что-то становится странно, например, q = 50.0 и t = 100.0, и я начинаю его вычислить из k = 0 to 100 я получаю серию из 0, за которой следуют NaN с определенного номера до конца.

Конечно, это вызвано операциями с числами, которые начинают слишком близко к 0 или подобным проблемам.

У вас есть идея, как я могу оптимизировать формулу, чтобы иметь возможность давать достаточно точные результаты по широкому распространению входов?

Все должно быть уже 64 бит (поскольку я использую OCaml, который по умолчанию использует удвоение). Возможно, есть способ использовать 128-битные двойники, но я не знаю, как это сделать.

Я использую OCaml, но вы можете предлагать идеи на любом языке, который вы хотите: C, C++, Java и т. Д. Я использовал их все.

+0

Что касается принятия 'log' + http://en.wikipedia.org/wiki/Stirling%27s_approximation – Anycorn

+0

sum (k = 0, ...) x^k/k!== exp (x), поэтому он почти выглядит так, как будто вы делаете exp (-qt) * exp (qt) = 1. Или я что-то упускаю? Заметим, что если qt велико, вы можете использовать exp (qt) = exp (qt/2)^2, т. Е. Вы можете разделить qt на 2 достаточно раза, чтобы сделать серию короткими, а затем квадратизировать это количество раз, чтобы получить желаемый ответ. Не уверен, что это полезно для того, что вы делаете. –

ответ

3
qt^k/k! = e^[log[qt^k/k!]] 
log[qt^k/k!] = log[qt^k] - log[k!] // log[k!] ~ klnk - k by stirling 
      ~ k ln(qt) - (k lnk - k) 
      ~ k ln(qt/k) - k 

для небольших значений k, приближение Стирлинга неточно. Однако, поскольку вы, кажется, используете конечный известный диапазон, вы можете вычислить log[k!] и поместить его в массив, избегая каких-либо ошибок.

, конечно, есть несколько вариантов, которые вы можете сделать дальше.

+0

Я пробовал аналитический apporach, как вы предлагали, но до сих пор не удалось получить аналогичные результаты. Я буду продолжать пытаться! В то же время я разложил операцию, чтобы держать в границах поплавки. – Jack

1

Это не ответ (я считаю), но дает только разъяснение. Если я что-то не понял, я удалю его после вашего комментария.

Как я понимаю, вы пытаетесь вычислить п, например, в виде следующей суммы равен 1.

alt text

Как вы можете видеть, что приближается к 1 асимптотически, он никогда не будет EQUAL до 1. Пожалуйста, исправьте меня, если я неправильно понял ваш вопрос.

+0

Да, вы правы, я использовал равные, но я имел в виду достаточно близко к 1.0. Фактически я хочу выбрать ошибку, например, например, 'e = 1e-6', чтобы' sum> 1 - e' .. – Jack

Смежные вопросы