Как оптимизировать это вычисление - PullRequest
Новая
       12

Как оптимизировать это вычисление

2 голосов
/ 03 сентября 2010

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

! [Alt text] [1]

где q является двойным, t тоже двойным и k 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 и т. Д. Iвполне использовал их все.

Ответы [ 2 ]

3 голосов
/ 03 сентября 2010 
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!] и поместить его в массив, избегая каких-либо ошибок.

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

1 голос
/ 04 сентября 2010

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

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

alt text

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

...