Ошибка при использовании оптимальной функции для уравнения правдоподобия на R

Question

Я новичок в R и столкнулся с ошибкой при попытке использовать функцию optim.

У меня есть уравнение вероятности, которое я хотел бы максимизировать, поэтому я реализовал следующий код:

>datafile=read.delim("clipboard")

> log.lik=function(theta, x, mu){
+ b=theta[1]
+ beta=theta[2]
+ tau=theta[3]
+ result=function(b,beta, tau){(-sum(x)/b)-sum(log(-beta*(x-tau)))-sum(log(integrate(exp(-x/b)/(1+exp(-beta(x-tau)))), lower=1500, upper=Inf))}
+ return(result)
+ }
> estimate=c(1,1,1)
> model=optim(par=estimate, fn=log.lik, control=list(fnscale=-1), x=datafile, mu=1500)

Все работает, пока функция optim, при которой я получаю следующее сообщение об ошибке: Ошибка в optim(par = эстимейт, fn = log.lik, control = список (fnscale = -1),: не может привести тип "замыкание" к вектору типа "двойной"

Кто-нибудь может понять, в чем здесь проблема? Любая помощь будет оценена!

Файл данных - это всего лишь один столбец имитируемых финансовых потерь в формате CSV. Когда я вывожу переменную datafile, вот пример того, что я получаю:

       X1946774
1      34949037
2     734018898
3     393502463
4     388573133
5      93213300
6      74982868
7      55322550
8      10828207
9       4530577
10      3786748
11      2041762
12    342745985
13    292313639
14    259569928
15    143871771
16     53691635
17     24489644
18     20506718
19     14281945

Отредактированный код, включающий изменения из комментариев:

 > log.lik=function(theta,x,mu){
    + b=theta[1]
    + beta=theta[2]
    + tau=theta[3]
    + integrand<-function(x,b,beta,tau){exp(-x/b)/(1+exp(-beta*(x-tau)))}
    + result<-(-sum(x)/b)-sum(log(-beta*(x-tau)))-sum(log(integrate(integrand, lower=mu, upper=Inf)))
    + return(result)
    + }
    > model=optim(par=estimate, fn=log.lik, control=list(fnscale=-1), x=datafile, mu=1500)

Answer 1

Слишком долго для комментария.

Во-первых, вы все еще используете integrate(...) неправильно. Эта функция возвращает список (читайте документацию!!). $value Элементом этого списка является интеграл. Таким образом, чтобы найти интеграл f от a в b, используйте:

integrate(f,a,b,...)$value

К сожалению, это наименьшая из ваших проблем. По сути, вы берете слишком много ярлыков. Вы не можете просто собрать некоторый код - вам нужно обратить внимание на математику.

Например, вы составили график значений вашего integrand(...) функция для начальных значений theta, В диапазоне (mu,Inf)?? Если бы вы имели, вы бы увидели, что подынтегральная функция 0 в этом диапазоне, потому что mu=1500 а также b=1, а также exp(-1500/1) численно 0; следовательно, интеграл равен 0, а лог интеграла не определен. Кроме того, ваша целевая функция включает в себя термин log(-beta*(x-tau)), но для beta=tau=1, -beta*(x-tau) < 0 для всех x в вашем наборе данных, и снова журнал не определен.

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