[問題] optim()

看板R_Language作者 (廢物)時間11年前 (2013/12/06 10:47), 編輯推噓0(005)
留言5則, 2人參與, 最新討論串1/1
[軟體熟悉度]: 入門(寫過其他程式,只是對語法不熟悉) [問題敘述]: optim()偵測到的變數維度與預料中不同 各位板上先進好~ 我現在是經濟系的研究助理~ 要幫老闆用MitISEM這個package去做Baysian Estimation~ 其中我的目標函數的程式碼如下: SChang <- function(theta, y, Z, w, t, dim, log=TRUE){ if (is.vector(theta)) theta <- matrix(theta, nrow = 1) n <- length(y) bx <- theta[,1] bz <- matrix(theta[,2:(1+dim[1])],nrow=dim[1],ncol=1) mu_x <- theta[,2+dim[1]]*matrix(rep(1,n),nrow=n, ncol=1) del <- matrix(theta[,(3+dim[1]):(2+dim[1]+dim[2])],nrow=dim[2],ncol=1) sig_v <- abs(theta[,3+dim[1]+dim[2]]) ## sig_v^2 sig_u <- abs(theta[,4+dim[1]+dim[2]]) ## sig_u^2 sig_epo <- abs(theta[,5+dim[1]+dim[2]]) ## sig_epo^2 sig_x <- abs(theta[,6+dim[1]+dim[2]]) ## sig_x^2 u <- matrix(theta[,(7+dim[1]+dim[2]):(6+dim[1]+dim[2]+n)],nrow=n, ncol=1) x <- matrix(theta[,(7+dim[1]+dim[2]+n):(6+dim[1]+dim[2]+n+n)], nrow=n, ncol=1) if (prod(u >0) > 0) p <- -(n*log(sig_v^0.5)+n*log(sig_u^0.5)+n*log(sig_x^0.5)+n*log(sig_epo^0.5)+ (y+u-Z %*% bz-x*bx) %*% t(y+u-Z %*% bz-x*bx)/(2*sig_v)+ (w-x) %*% t(w-x)/(2*sig_epo)+ (x-mu_x) %*% t(x-mu_x)/(2*sig_x)+ (u-t %*% del) %*% t(u-t %*% del)/(2*sig_u)+ sum(log(pnorm(u/sig_u^0.5)))) else p <- -Inf if (!log) p <- exp(p) as.vector(p) } 其中y, Z, w, t是我的data~ theta在我的case中應該是個2014維的向量~ 也是我要估計的函數的變數~ 接下來我就把它丟到MitISEM中做估計: app.Schang <- MitISEM(KERNEL=SChang, mu0=mu0, y=y, Z=Z, w=w, t=t, dim=dim) 結果顯示: 錯誤在optim(par = mu0, fn = KERNEL, method = method, control = control, : optim 內的目的函數長度被評估為 998001,而不是 1 是我目標函數寫錯嗎? 998001與我預期的2014相去甚遠~ 我到底是哪裡做錯了? [關鍵字]: optim() -- -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.112.4.183 ※ 編輯: Dboy 來自: 140.112.4.183 (12/06 10:48)

12/06 11:28, , 1F
p 回傳的是 logL value, 不是vector of length(y).
12/06 11:28, 1F

12/06 11:31, , 2F
這是MitISEM這個package的要求~它要求要有取log這個動作~所以
12/06 11:31, 2F

12/06 11:31, , 3F
我才這樣寫~還是說我誤會1樓的意思了?
12/06 11:31, 3F

12/06 11:33, , 4F
而且optim本來就要求目標函數傳回scalar而不是向量不是嗎?
12/06 11:33, 4F

12/06 12:27, , 5F
我發現我transpose寫錯地方......感謝各位~
12/06 12:27, 5F
文章代碼(AID): #1IeJh2wt (R_Language)
文章代碼(AID): #1IeJh2wt (R_Language)