作者KirinGuess (Kirin)
看板Statistics
標題Re: [問題] R跑probit模型
時間Tue Jan 12 16:55:02 2010
※ 引述《k84526991 (哇哈哈)》之銘言：
: 要如何用r跑probit模型呢??
: 看計量理論的書 必須要用maxliklihood 去找估計系數值
: 請問這要麼用r跑呢??
: 我不要用glm(.....)
: 有沒有人可以教我詳細的語法
以下是我以前練習寫的，
但我不敢確定正不正確，
特別是在的估計MLE的地方。
需要的話您可以參考一下。
如果有錯或可以改得更精簡的地方，
也請大家給我些意見，謝謝。
----------------------------------------------------------------------------
我們先隨機產生一個資料檔(DLFP)，
包括一個依變數和七個自變數。
LFP <-rbinom(1000,1,0.5)
K5 <-rbinom(1000,3,0.5)
K618<-rbinom(1000,8,0.5)
AGE <-round(runif(1000, min=30, max=60))
WC <-rbinom(1000,1,0.5)
HC <-rbinom(1000,1,0.5)
LWG <-rbinom(1000,6,0.5) -3
INC <-round(runif(1000, min=0, max=96))
DLFP<-data.frame(cbind(LFP,K5,K618,AGE,WC,HC,LWG,INC))
XS<-model.matrix(LFP~K5+K618+AGE+WC+HC+LWG+INC,
DLFP) #DATA MATRIX
以下是LOGIT模型的估計：
FN<-function(beta)
-1*(t(LFP)%*%log(plogis(XS%*%beta))+
t(1-LFP)%*%log(1-plogis(XS%*%beta))
) #-1*Likelihood function
MLE<-nlm(FN,p=rep(0,ncol(XS)),
hessian=T,iterlim=1000)
#Newton-type algorithm
NAMES<-as.matrix(labels(XS)[[2]]) #Names of Indep. Variables
COEF <-MLE$estimate #Coef.
SE <-sqrt(diag(solve(MLE$hessian))) #Std. Error
Z <-COEF/SE #Z
P <-pnorm(abs(Z),lower.tail=F)*2 #P-value
LOGIT<-matrix(c(COEF,SE,Z,P),
nrow=length(NAMES),
dim=list(c(NAMES),
c("Coef.","S.E.","Z","P-value"))
) #Table
LOGIT
以下是PROBIT模型的估計：
FN2<-function(beta)
-1*sum(t(LFP)%*%log(pnorm(XS%*%beta))+
t(1-LFP)%*%log(1-pnorm(XS%*%beta))
) #-1*Likelihood function
MLE2<-nlm(FN2,p=rep(0,ncol(XS)),
hessian=T,iterlim=1000)
#Newton-type algorithm
NAMES2<-as.matrix(labels(XS)[[2]]) #Names of Indep. Variables
COEF2 <-MLE2$estimate #Coef.
SE2 <-sqrt(diag(solve(MLE2$hessian))) #S.E.
Z2 <-COEF2/SE2 #Z
P2 <-pnorm(abs(Z2),lower.tail=F)*2 #P-value
PROBIT<-matrix(c(COEF2,SE2,Z2,P2),
nrow=length(NAMES2),
dim=list(c(NAMES2),
c("Coef.","S.E.","Z","P-value"))
) #TABLE
PROBIT
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 114.33.213.179
1F：推 k84526991:感謝 我試試 01/12 17:59
※ 編輯: KirinGuess 來自: 140.123.197.80 (02/03 19:46)