dat<- read.table ("Contracts_Red.txt",header=TRUE) Anzahl<- c(4,3,3,3,3,2,2) Zielstates<- read.table("Zielstates.txt", header=FALSE) dist <- function (x) { mean<- max(0.00001,x[2]) sd<- max(0.00001,x[3]) w<- plnorm(x[1]+0.5, meanlog = mean, sdlog = sd, lower.tail = FALSE, log.p = FALSE) w } dist1 <- function (x) { mean<- max(0.00001,x[2]) sd<- max(0.00001,x[3]) if (x[1]<=0.5) w<- 1 else w<- plnorm(x[1]-0.5, meanlog = mean, sdlog = sd, lower.tail = FALSE, log.p = FALSE) w } dens <- function (x) { w<- dist1(x)-dist(x) w } L<- function(z) { #z[1]=state, z[i]=parameter for transitions to other states LL<- 0 for (i in 1:5199) if (dat[i,25]==z[1]){ for (m in 1:Anzahl[z[1]+1]) { x<- c(dat[i,20],z[2*m],z[2*m+1]) if (Zielstates[z[1]+1,m]==dat[i,24]) LL<- LL-log(dens(x)) else LL<- LL-log(dist(x))} if (dat[i,22]>0) for (m in 1:Anzahl[z[1]+1]){ x<- c(dat[i,22],z[2*m],z[2*m+1]) LL<- LL-log(dist(x))} } return(LL) } B <- array(0,c(7,9)) C <- array(0,c(7,8)) for (n in 1:7){ #estimation State<- n-1; Min<- 100000; Lmin<- function(y) { z<- c(State,y) r<- L(z) return(r) } dim<- 2*Anzahl[n] a<- array(0,c(dim)) anf<- array(0,c(dim)) u<- array(1,c(dim)) if (n==1) anf<- c(4.857793, 1.280511, 4.515230, 1.264155, 4.643286, 1.256131, 3.85951, 1.398584) if (n==2) anf<- c(3.986773, 1.208599, 3.803702, 1.179337, 2.982379, 1.953849) if (n==3) anf<- c(4.644558, 1.488709, 3.812662, 1.188853, 2.747024, 1.917147) if (n==4) anf<- c(4.482829, 1.488381, 3.695483, 1.160581, 4.004806, 0.930976) if (n==5) anf<- c(4.468907, 1.472684, 4.421664, 1.447728, 3.323349, 1.943854) if (n==6) anf<- c(4.177719, 1.335833, 2.868238, 1.269539) if (n==7) anf<- c(4.576287, 1.589189, 2.460096,1.162619) for (iter in 1:2) { resL <- optim(par=anf,fn=Lmin,gr=NULL,method=("BFGS"),control=list(trace=1,fnscale=1,parscale=u,alpha=1,beta=0.5,gamma=2,abstol=0.00001,reltol=0.00001), hessian=FALSE) anf<- resL$par for (k in 1:(2*Anzahl[n])) if (anf[k]<0) anf[k]<- 0.01 res <- optim(par=anf,fn=Lmin,gr=NULL,method=("Nelder-Mead"),control=list(trace=1,fnscale=1,parscale=u,alpha=1,beta=0.5,gamma=2,abstol=0.00001,reltol=0.00001), hessian=FALSE) anf<- res$par for (k in 1:(2*Anzahl[n])) if (anf[k]<0) anf[k]<- 0.1 } resL <- optim(par=anf,fn=Lmin,gr=NULL,method=("BFGS"),control=list(trace=1,fnscale=1,parscale=u,alpha=1,beta=0.5,gamma=2,abstol=0.00001,reltol=0.00001), hessian=TRUE) d<-det(100*resL$hessian) if (abs(d)>10^-20) a<-diag(solve(resL$hessian)) z<- resL$par for (m in 1:(2*Anzahl[n])){ B[n,m]<- z[m] C[n,m]<- if (a[m]<0) -1 else sqrt(a[m])*10 } B[n,9]<- resL$val } print(B) print(C) # source("LogNorm.R")