R: maximum entropy regression linear and non linear
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
No Responses Yet to “R: maximum entropy regression linear and non linear”