作者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)