作者f496328mm (為什麼會流淚)
看板R_Language
標題[問題] MLE 做 聯立方程式 估計 (已解決)
時間Wed Sep 7 15:56:29 2016
聯立方程式 or system equation
用MLE去做 我目前只會做一條equation的MLE
我目前上網找的只能做到一條equation的MLE估計
目前code如下:
#==========================================
library(stats4)
set.seed(1001)
N <- 100
x <- runif(N)
y <- 5 * x + 3 + rnorm(N)
LL <- function(beta0, beta1, mu, sigma) {
# Find residuals
#
R = y - x * beta1 - beta0
#
# Calculate the likelihood for the residuals (with mu and sigma as parameters)
#
R = suppressWarnings(dnorm(R, mu, sigma))
#
# Sum the log likelihoods for all of the data points
#
-sum(log(R))
}
fit <- mle(LL, start = list(beta0 = 3, beta1 = 1, mu = 0, sigma=1))
summary(fit)
#==========================================
R = y - x * beta1 - beta0
這是我的一條equation
我目前只會做到一條
如果要做一條以上，要如何做?
#=========================================
另外我也有找到這個package "systemfit"
它可以一次做兩個以上的迴歸
不過 我的理解是 他是分開做的
因為我分別用lm去做 估計出來的參數也一樣
所以沒什麼意義 可能是我理解錯誤
不知道這個函數的用意是什麼
謝謝
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 114.42.98.232
※ 文章網址: https://webptt.com/m.aspx?n=bbs/R_Language/M.1473234992.A.2CE.html
1F：推 ksherry: Google Full information maximum likelihood 09/07 17:23
2F：→ ksherry: 聯立方程組 你還是要自己先求出均衡解 09/07 17:23
3F：推 qhair: 我猜你要的是BBsolve 09/07 17:39
不是BBsolve喔 剛剛看了一下 不過還是謝謝
因為是估計值 所以至少要有Std. Error p-value
我正在試試看分開用nls估計 看差異大不大
只要還是想一起估計 因為互相為限制式 可是我不會多條equation估計><
※ 編輯: f496328mm (114.42.98.232), 09/07/2016 19:05:26
自己的問題自己回
研究出來啦 煩了我好幾個禮拜
http://imgur.com/xUatDl5
先做模擬預設要估計的beta 兩條equation
最後估計出近似原始設的beta
code如下:
library(stats4)
set.seed(1001)
N <- 100
x1 <- rnorm(N,10,1)
x2 <- rnorm(N,5,1)
beta0=3
beta1=5
beta2=50
y1 <- beta1 * x1 + beta0 + rnorm(N)
y2 <- beta2 * x2 + beta1*x1 + beta0 + rnorm(N)
LL <- function(beta0, beta1,beta2, mu, sigma) {
R1 = y1 - x1 * beta1 - beta0
R2 = y2 - x2 * beta2 - x1*beta1 - beta0
#
R1 = suppressWarnings(dnorm(R1, mu, sigma, log = TRUE))
R2 = suppressWarnings(dnorm(R2, mu, sigma, log = TRUE))
#
-sum(R1,R2)
}
fit <- mle2(LL, start = list( beta0 = 3, beta1 = 1,
beta2 = 1, mu = 0, sigma = 1))
summary(fit)
※ 編輯: f496328mm (114.42.98.232), 09/07/2016 19:42:41
※ 編輯: f496328mm (114.42.98.232), 09/07/2016 19:42:54
4F：→ celestialgod: 發文格式... 09/07 20:19