This code is for implementation of maximum entropy non linear and linear estimation of regression in R. You will need Alabama package.
## PARA REGRESIONES LINEALES.
## Funcion GENERALIZED MAXIMUM ENTROPY
gme.v1=function(Y=Y,X=X,M=M,C=C,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
## Support Matrix of dim k x M.
if (opt2==”noprior”) {
z1=rep(0,M)
nz1=floor(length(z1)/2)
dum<-1/nz1
dum<-seq(0,1,by=(dum))*C
z1<-c(-rev(dum[2:(nz1+1)]),dum)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## inicial po.
p1o=rep(1,length(z1))
p1o=p1o/sum(p1o)
po<-rep(p1o,dim(X)[2])
}
if (opt2==”prior”) {
z1=rep(0,M)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## rownames
dum1<-rep(“beta”,dim(X)[2])
dum2<-seq(1:dim(X)[2])
dum1<-paste(dum1,dum2,sep=”")
rownames(Z)<-dum1
## colnames
dum1<-rep(“p”,M)
dum2<-seq(1:M)
dum1<-paste(dum1,dum2,sep=”")
dum2<-rep(dum1,dim(X)[2])
colnames(Z)<-dum2
##
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
dum=dum+(M)
}
## inicial po.
p1o=rep(0,length(z1))
po<-rep(p1o,dim(X)[2])
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
po[dum1:(dum1+(M-1))]<-matprior2[i,]
dum=dum+M
}
}
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
V=kronecker(diag(1,length(Y)), t(v1))
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))
## Valor inicial total.
pwo<-c(po,wo)
## k
k=length(po)
## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
## tama?o de restricciones.
res.1<-dim(V)[2]
## restricciones de probabilidades.
res.2<-2
res<-res.1+res.2
h<-rep(NA,res)
## Restricciones lineales.
h[1:res.1]<-X%*%(Z%*%p)+(V%*%w)-Y
## Restricciones de probabilidades.
h[(res.1+1)]<-sum(p,na.rm=TRUE)-1
h[(res.1+2)]<-sum(w,na.rm=TRUE)-1
h
}
## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
h<-rep(NA,(2*len))
## Restricciones lineales.
h[1:len]<-pw
h[(len+1):(2*len)]<- -pw+1
}
H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
-Objective
}
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Grad.p<-log(p)+1
Grad.w<-log(w)+1
Grad<-c(Grad.p,Grad.w)
Grad
}
## El default de auglag y optim es minimizar.
################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}
if (opt1==”gr”) {
## con gradiente
results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]
e<-V%*%results$par[(k+1):length(results$par)]
list(beta=beta,e=e)
}
## Generalized Maximum Entropy no lineal Especializada para
gme.nl=function(Y=Y,X=X,M=M,C=C,fun=fun,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
## inicializacion de I y L
I=X[,"I"]
L=X[,"L"]
## Support Matrix of dim k x M.
if (opt2==”noprior”) {
z1=rep(0,M)
nz1=floor(length(z1)/2)
dum<-1/nz1
dum<-seq(0,1,by=(dum))*C
z1<-c(-rev(dum[2:(nz1+1)]),dum)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## inicial po.
p1o=rep(1,length(z1))
p1o=p1o/sum(p1o)
po<-rep(p1o,dim(X)[2])
}
if (opt2==”prior”) {
z1=rep(0,M)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## rownames
dum1<-rep(“beta”,dim(X)[2])
dum2<-seq(1:dim(X)[2])
dum1<-paste(dum1,dum2,sep=”")
rownames(Z)<-dum1
## colnames
dum1<-rep(“p”,M)
dum2<-seq(1:M)
dum1<-paste(dum1,dum2,sep=”")
dum2<-rep(dum1,dim(X)[2])
colnames(Z)<-dum2
##
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
dum=dum+(M)
}
## inicial po.
p1o=rep(0,length(z1))
po<-rep(p1o,dim(X)[2])
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
po[dum1:(dum1+(M-1))]<-matprior2[i,]
dum=dum+M
}
}
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
V=kronecker(diag(1,length(Y)), t(v1))
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))
## Valor inicial total.
pwo<-c(po,wo)
## k
k=length(po)
## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
## tama?o de restricciones.
res.1<-dim(V)[2]
## restricciones de probabilidades.
res.2<-2
res<-res.1+res.2
h<-rep(NA,res)
## Restricciones lineales.
theta<-(Z%*%p)
h[1:res.1]<-fun(theta=theta,X=X)+(V%*%w)-Y
## Restricciones de probabilidades.
h[(res.1+1)]<-sum(p,na.rm=TRUE)-1
h[(res.1+2)]<-sum(w,na.rm=TRUE)-1
h
}
## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
h<-rep(NA,(2*len))
## Restricciones lineales.
h[1:len]<-pw
h[(len+1):(2*len)]<- -pw+1
}
H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
-Objective
}
## El default de auglag y optim es minimizar.
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Grad.p<-log(p)+1
Grad.w<-log(w)+1
Grad<-c(Grad.p,Grad.w)
Grad
}
## El default de auglag y optim es minimizar.
################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
if (opt1==”gr”) {
## con gradiente
results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]
e<-V%*%results$par[(k+1):length(results$par)]
list(beta=beta,e=e)
}
#gme.nl=cmpfun(gme.nl)
##
##
## Generalized Maximum Entropy no lineal v2 especializado para stock de capital.
gme.nl2=function(Y=Y,X=X,M=M,C=C,fun=fun,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
## inicializacion de I y L
I=X[,"I"]
L=X[,"L"]
## Support Matrix of dim k x M.
if (opt2==”noprior”) {
z1=rep(0,M)
nz1=floor(length(z1)/2)
dum<-1/nz1
dum<-seq(0,1,by=(dum))*C
z1<-c(-rev(dum[2:(nz1+1)]),dum)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## inicial po.
p1o=rep(1,length(z1))
p1o=p1o/sum(p1o)
po<-rep(p1o,dim(X)[2])
}
if (opt2==”prior”) {
z1=rep(0,M)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## rownames
dum1<-rep(“beta”,dim(X)[2])
dum2<-seq(1:dim(X)[2])
dum1<-paste(dum1,dum2,sep=”")
rownames(Z)<-dum1
## colnames
dum1<-rep(“p”,M)
dum2<-seq(1:M)
dum1<-paste(dum1,dum2,sep=”")
dum2<-rep(dum1,dim(X)[2])
colnames(Z)<-dum2
##
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
dum=dum+(M)
}
## inicial po.
p1o=rep(0,length(z1))
po<-rep(p1o,dim(X)[2])
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
po[dum1:(dum1+(M-1))]<-matprior2[i,]
dum=dum+M
}
}
## Matriz para obtener restricciones de probabilidades
z.id=rep(1,M)
Idp=kronecker(diag(1,dim(X)[2]), t(z.id))
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
V=kronecker(diag(1,length(Y)), t(v1))
## Matriz para obtener restricciones de probabilidades
v1.id=rep(1,length(v1))
Idw=kronecker(diag(1,length(Y)), t(v1.id))
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))
## Valor inicial total.
pwo<-c(po,wo)
## k
k=length(po)
## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
## tama?o de restricciones.
res.1<-length(Y)
## restricciones de probabilidades.
res.2<-dim(Z)[1]
res<-res.1+res.2+res.1
h<-rep(NA,res)
## Restricciones lineales.
theta<-(Z%*%p)
h[1:res.1]<-fun(theta=theta,X=X)+(V%*%w)-Y
## Restricciones de probabilidades.
h[(res.1+1):(res.1+(res.2))]<-as.numeric(Idp%*%p-rep(1,times=(res.2/2)))
h[(res.1+(res.2)+1):length(h)]<-as.numeric(Idw%*%w-rep(1,times=res.1))
h
}
## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
h<-rep(NA,(2*len))
## Restricciones lineales.
h[1:len]<-pw
h[(len+1):(2*len)]<- -pw+1
}
H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
-Objective
}
## El default de auglag y optim es minimizar.
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Grad.p<-log(p)+1
Grad.w<-log(w)+1
Grad<-c(Grad.p,Grad.w)
Grad
}
## El default de auglag y optim es minimizar.
################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
if (opt1==”gr”) {
## con gradiente
results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]
e<-V%*%results$par[(k+1):length(results$par)]
Hn=-results$value/log(length(results$par))
list(beta=beta,e=e,Hn=Hn)
}
#gme.nl=cmpfun(gme.nl)
gme.v2=function(Y=Y,X=X,M=M,C=C,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
## Support Matrix of dim k x M.
if (opt2==”noprior”) {
z1=rep(0,M)
nz1=floor(length(z1)/2)
dum<-1/nz1
dum<-seq(0,1,by=(dum))*C
z1<-c(-rev(dum[2:(nz1+1)]),dum)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## inicial po.
p1o=rep(1,length(z1))
p1o=p1o/sum(p1o)
po<-rep(p1o,dim(X)[2])
}
if (opt2==”prior”) {
z1=rep(0,M)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## rownames
dum1<-rep(“beta”,dim(X)[2])
dum2<-seq(1:dim(X)[2])
dum1<-paste(dum1,dum2,sep=”")
rownames(Z)<-dum1
## colnames
dum1<-rep(“p”,M)
dum2<-seq(1:M)
dum1<-paste(dum1,dum2,sep=”")
dum2<-rep(dum1,dim(X)[2])
colnames(Z)<-dum2
##
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
dum=dum+(M)
}
## inicial po.
p1o=rep(0,length(z1))
po<-rep(p1o,dim(X)[2])
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
po[dum1:(dum1+(M-1))]<-matprior2[i,]
dum=dum+M
}
}
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
V=kronecker(diag(1,length(Y)), t(v1))
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))
## Valor inicial total.
pwo<-c(po,wo)
## k
k=length(po)
## matrices de restricciones
## Matriz para obtener restricciones de probabilidades
v1.id=rep(1,length(v1))
Idw=kronecker(diag(1,length(Y)), t(v1.id))
## Matriz para obtener restricciones de probabilidades
z.id=rep(1,M)
Idp=kronecker(diag(1,dim(X)[2]), t(z.id))
## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
## tama?o de restricciones.
res.1<-length(Y)
## restricciones de probabilidades.
res.2<-dim(Z)[1]
res<-res.1+res.2+res.1
h<-rep(NA,res)
## Restricciones lineales.
h[1:res.1]<-X%*%(Z%*%p)+(V%*%w)-Y
## Restricciones de probabilidades.
## Restricciones de probabilidades.
h[(res.1+1):(res.1+(res.2))]<-as.numeric(Idp%*%p-rep(1,times=(res.2/2)))
h[(res.1+(res.2)+1):length(h)]<-as.numeric(Idw%*%w-rep(1,times=res.1))
h
}
## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
h<-rep(NA,(2*len))
## Restricciones lineales.
h[1:len]<-pw
h[(len+1):(2*len)]<- -pw+1
}
H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
-Objective
}
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Grad.p<-log(p)+1
Grad.w<-log(w)+1
Grad<-c(Grad.p,Grad.w)
Grad
}
## El default de auglag y optim es minimizar.
################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}
if (opt1==”gr”) {
## con gradiente
results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]
e<-V%*%results$par[(k+1):length(results$par)]
Hn<–results$value
list(beta=beta,e=e,Hn=Hn)
}
#———————————————————————-
# entropia normalizada, recordad que la distribucion uniforme maximiza la entropia 1/n, entonces
# se normaliza la entropia dividiendo para ln n . Hn = H/ln n
# -results$fn/log(dim(X)[1])
## Generalized Maximum Entropy no lineal v2 especializado para stock de capital.
gme.nl3=function(Y=Y,X=X,M=M,C=C,fun=fun,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2,fix.t=fix.t) {
## inicializacion de I y L
I=X[,"I"]
L=X[,"L"]
## Support Matrix of dim k x M.
if (opt2==”noprior”) {
z1=rep(0,M)
nz1=floor(length(z1)/2)
dum<-1/nz1
dum<-seq(0,1,by=(dum))*C
z1<-c(-rev(dum[2:(nz1+1)]),dum)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## inicial po.
p1o=rep(1,length(z1))
p1o=p1o/sum(p1o)
po<-rep(p1o,dim(X)[2])
}
if (opt2==”prior”) {
z1=rep(0,M)
Z=kronecker(diag(1,dim(X)[2]), t(z1))
## rownames
dum1<-rep(“beta”,dim(X)[2])
dum2<-seq(1:dim(X)[2])
dum1<-paste(dum1,dum2,sep=”")
rownames(Z)<-dum1
## colnames
dum1<-rep(“p”,M)
dum2<-seq(1:M)
dum1<-paste(dum1,dum2,sep=”")
dum2<-rep(dum1,dim(X)[2])
colnames(Z)<-dum2
##
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
dum=dum+(M)
}
## inicial po.
p1o=rep(0,length(z1))
po<-rep(p1o,dim(X)[2])
dum=0
for (i in 1:dim(X)[2]) {
dum1=1+dum
po[dum1:(dum1+(M-1))]<-matprior2[i,]
dum=dum+M
}
}
## Matriz para obtener restricciones de probabilidades
z.id=rep(1,M)
Idp=kronecker(diag(1,dim(X)[2]), t(z.id))
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
V=kronecker(diag(1,length(Y)), t(v1))
## Matriz para obtener restricciones de probabilidades
v1.id=rep(1,length(v1))
Idw=kronecker(diag(1,length(Y)), t(v1.id))
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))
## Valor inicial total.
pwo<-c(po,wo)
## k
k=length(po)
## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
## tama?o de restricciones.
res.1<-length(Y)
## restricciones de probabilidades.
res.2<-dim(Z)[1]
res<-res.1+res.2+res.1
h<-rep(NA,res)
## Restricciones lineales.
theta<-(Z%*%p)
h[1:res.1]<-fun(theta=theta,X=X)+(V%*%w)-Y
## Restricciones de probabilidades.
h[(res.1+1):(res.1+(res.2))]<-as.numeric(Idp%*%p-rep(1,times=(res.2/2)))
h[(res.1+(res.2)+1):length(h)]<-as.numeric(Idw%*%w-rep(1,times=res.1))
# fix.t
hf<-as.numeric(V%*%w)[fix.t]
h<-c(h,hf)
h
}
## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
h<-rep(NA,(2*len))
## Restricciones lineales.
h[1:len]<-pw
h[(len+1):(2*len)]<- -pw+1
}
H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
-Objective
}
## El default de auglag y optim es minimizar.
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
len=length(pw)
p=pw[1:k]
w=pw[(k+1):len]
Grad.p<-log(p)+1
Grad.w<-log(w)+1
Grad<-c(Grad.p,Grad.w)
Grad
}
## El default de auglag y optim es minimizar.
################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
if (opt1==”gr”) {
## con gradiente
results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]
e<-V%*%results$par[(k+1):length(results$par)]
Hn=-results$value/log(length(results$par))
list(beta=beta,e=e,Hn=Hn)
}
Filed under: R | Leave a Comment
Ubuntu: corrupted Synaptic.
If your Synaptic is corrupted you can follow the instructions in this forum
Filed under: linux | Leave a Comment
This post draws heavily from andrewdyck blog .
You need unix installation files for Stata 11 of course, then open a terminal and:
sudo mkdir /usr/local/stata10
cd /usr/local/stata10
sh /media/cdrom/install
Dont forget to give permissions to read and write and create and delete files to this folder. sudo chmod -R 755 /usr/local/stata10
Alternatively, use gksudo nautilus and create the required folder /usr/local/stata10. Make a right click over stata10 folder and choose properties, and give all the permissions. No other folder is allowed.
After you choose the required options, in Ubuntu first 3. Linux then the version of Stata you have, either Linux 32 or Linux 64.
You’re not finished yet, you need to make the symbolic link necessary to run the Stata gui:
sudo ln -s /usr/lib64/x86_64-linux-gnu/libtiff.so.4.3.3 /usr/lib/libtiff.so.3
If, you’re running a different version of ubuntu you have to locate libtiff.so and make a symbolic link with the required version of Stata (if it already exists libtiff.so.4.3.3 or whatever version there is -i have not tested others but I know v 4 works). If the file doesn’t exist install it. (Synaptic)
Finally, you need the libgtksourceviewmm-1.0-dev through Synaptic or sudo apt-get it!
In the /usr/local/stata10 folder run
./xstata you have the gui, and then run update all
If you want the console
./stata
Finally, you can install ess emacs speaks statistics from the repositories and emacs will help you edit your scripts.
Enjoy, Stata under ubuntu.
Filed under: linux, Stata, ubuntu | Leave a Comment
python: install openopt
1. Download MinGW32
2. Add to the PATH “C:\Mingw32\bin”
3. Download openopt package from the author’s web
4. open a prompt and cd the directory where the files had been uncompressed
5. execute in the openopt folder
setup.py config –compiler=mingw32 build –compiler=mingw32
Filed under: Python | Leave a Comment
This assumes you have installed Rtools and that R can compile packages from source.
1. Download the source files nloptr-src,
from the author’s web. Unzip this in any file for example C:\nloptr.
2. Download the binaries from NLopt
that correspond to the lastest version. This link refers to nlopt-2.2.3-dll.zip.
3. Copy the whole content of nlopt-2.2.3-dll.zip. in C:\nloptr\src\nlopt-2.2
4. In R: run
install.packages(“C:/nloptr”,repos=NULL,type=”source”)
You will receive a warning that the file libnlopt-0.dll is not installed in your system, don’t worry and do the following step.
5. A new folder will be created with the name src-i386
Copy the file C:\nloptr\src\nlopt-2.2\libnlopt-0.dll into C:\nloptr\src-i386\nlopt-2.2\libnlopt-0.dll.
6. In R: run
install.packages(“C:/nloptr”,repos=NULL,type=”source”)
7. It’s a messy intallations but it works fine. If you have better options please let me know.
Filed under: R | Leave a Comment
Excel to LaTex, xl2Latex can be downloaded from the web, it´s a .xla that enables a excel complement that converts an excel table that is selected with the mouse to a LaTex code, that can be copied in the clipboard or saved in a file.
Filed under: LaTex | Leave a Comment
This post is based completely (except for the bold characters) in the thread in http://stackoverflow.com/questions/5223113/using-mysql-in-r-for-windows,
about how to install RMySQL in Windows 7 with MySQL 5.5
The answer is due to:
Yuriy Petrovskiy
“Found solution with help of ran2, who gave me link to common question. The basic process is described here, but there are several hints, So I will describe the whole solution (please change the R version and paths if needed):
1. Install latest RTools from here
2. install MySQL or header and library files of mysql
3. create or edit file C:\Program Files\R\R-2.12.1\etc\Renviron.site and add line like MYSQL_HOME=C:/mysql (path to your mysql files)
4. copy libmysql.lib from mysql/lib to mysql/lib/opt to meet dependencies.
5. copy libmysql.dll to C:\Program Files\R\R-2.13.0\bin\i386 or to windows/system32 directory.
6. run install.packages(‘RMySQL’,type=’source’) and wait while compilation will end.
“
Hope it helps.
Filed under: R | Leave a Comment
The task, recode or create a new variable with a new codification from an original one. In this example we have ciiu industrial code to a user more aggregated groups in string.
Remember that the dictionary must have a one to one correspondence, that is
newcode originalcode
“sector1″ 1010
“sector1″ 1020
// The first step is to import the dictionary file:
insheet using “dict.txt”, tab
// Next, make the stata file into a matrix, if you need string codes use names, for example rownames.
mkmat originalcode, matrix(dict) rownames(stringvariable_newcode)
// Then use the database to be recoded. Do not use clear all because it will delete the matrix.
use “database.dta”, clear
// Identify the original variable to be recoded.
// Then use this code.
// drop if exists the variable where the new codes are going to be written.
capture drop recodmip
gen recodmip=”"
// calculates the row number of the matrix or the number of new codes.
local nrow=rowsof(dict)
forval i=1(1)`nrow’ {
mat b=dict[`i',"originalcode"]
local nam : rownames b
di “`nam’” // the new code in string
local codigo=b[1,1] // makes the matrix into a scalar.
capture replace recodmip=”`nam’” if originalcode==`codigo’
}
Filed under: Stata | Leave a Comment
R: gsub replacing a dot “.”
The gsub command in R v. 2.12.1 when applies the gsub(pattern=”.”,replacement=”",text) recognizes “.” as a wild card for replace all,
to solve this you must use gsub(pattern=”\\.”,replacement=”",text), to replace only the dots.
Filed under: Uncategorized | Leave a Comment
Rpy2 is a great Python package, nonetheless, it has problems running with newer versions of R (for instance 2.12.1), because R comes with a new directory structure. Specifically, it includes both 32 bits and 64 bits versions in separated folders, and the former lib folder was renamed to library. To adjust for this changes, the \rinterface\__init__.py script located under the folder of rpy2\rinterface has to be modified.
The lines that were added have a #ADDED comment next to them.
import os, sys
try:
R_HOME = (os.environ["R_HOME"], )
except KeyError:
R_HOME = os.popen("R RHOME").readlines()
if len(R_HOME) == 0:
if sys.platform == 'win32':
try:
import win32api
import win32con
hkey = win32api.RegOpenKeyEx(win32con.HKEY_LOCAL_MACHINE,
"Software\\R-core\\R",
0, win32con.KEY_QUERY_VALUE )
R_HOME = win32api.RegQueryValueEx(hkey, "InstallPath")[0]
win32api.RegCloseKey( hkey )
except:
raise RuntimeError(
"Unable to determine R version from the registery." +\
"Calling the command 'R RHOME' does not return anything.\n" +\
"This might be because R.exe is nowhere in your Path.")
else:
raise RuntimeError(
"R_HOME not defined, and no R command in the PATH."
)
else:
#Twist if 'R RHOME' spits out a warning
if R_HOME[0].startswith("WARNING"):
R_HOME = R_HOME[1]
else:
R_HOME = R_HOME[0]
R_HOME = R_HOME.strip()
os.environ['R_HOME'] = R_HOME
# Win32-specific code copied from RPy-1.x
if sys.platform == 'win32':
import win32api
#os.environ['PATH'] += ';' + os.path.join(R_HOME, 'bin')
#os.environ['PATH'] += ';' + os.path.join(R_HOME, 'modules')
#os.environ['PATH'] += ';' + os.path.join(R_HOME, 'lib')
os.environ['PATH'] += ';' + os.path.join(R_HOME, 'bin', 'i386') #ADD
os.environ['PATH'] += ';' + os.path.join(R_HOME, 'modules', 'i386') #ADD
os.environ['PATH'] += ';' + os.path.join(R_HOME, 'library') #ADD
# Load the R dll using the explicit path
# First try the bin dir:
Rlib = os.path.join(R_HOME, 'bin', 'R.dll')
# Then the lib dir:
#if not os.path.exists(Rlib):
# Rlib = os.path.join(R_HOME, 'lib', 'R.dll')
# Try bin/i386 subdirectory seen in R 2.12.0 ## ADDED ##
if not os.path.exists(Rlib): ## ADDED ##
Rlib = os.path.join(R_HOME, 'bin', 'i386', 'R.dll') ## ADDED ##
# Otherwise fail out!
if not os.path.exists(Rlib):
raise RuntimeError("Unable to locate R.dll within %s" % R_HOME)
win32api.LoadLibrary( Rlib )
# cleanup the namespace
del(os)
try:
del(win32api)
del(win32con)
except:
pass
from rpy2.rinterface.rinterface import *
class StrSexpVector(SexpVector):
"""
Vector of strings.
"""
def __init__(self, v):
super(StrSexpVector, self).__init__(v, STRSXP)
class IntSexpVector(SexpVector):
"""
Vector of integers.
"""
def __init__(self, v):
super(IntSexpVector, self).__init__(v, INTSXP)
class FloatSexpVector(SexpVector):
"""
Vector of floats.
"""
def __init__(self, v):
super(FloatSexpVector, self).__init__(v, REALSXP)
class BoolSexpVector(SexpVector):
"""
Vector of booleans (logical in R terminology).
"""
def __init__(self, v):
super(BoolSexpVector, self).__init__(v, LGLSXP)
# wrapper because print is strangely not a function
# Python prior to version 3.0
def consolePrint(x):
"""This is the default callback for R's console. It simply writes to stdout."""
sys.stdout.write(x)
setWriteConsole(consolePrint)
def consoleFlush():
sys.stdout.flush()
setFlushConsole(consoleFlush)
def consoleRead(prompt):
input = raw_input(prompt)
input += "\n"
return input
setReadConsole(consoleRead)
Filed under: Uncategorized | 6 Comments
Search
-
Blogroll
Recent Entries
- R: maximum entropy regression linear and non linear
- Ubuntu: corrupted Synaptic.
- Stata: Install Stata 10 MP dynamically linked (gui) in Ubuntu Natty
- python: install openopt
- R: Installing nloptr in R 2.13 under Windows 7.
- LaTex and Lyx: Excel to LaTex xl2Latex
- R: How to install RMySQL in Windows 7 with MySQL 5.5
- Stata: An easy way to recode a variable using a dictionary.
- R: gsub replacing a dot “.”
- Python: Rpy2 with Python 2.6 and R 2.12.1 R_HOME path problem.
- Python: Cubic Splines executable project