SCM

SCM Repository

[depmix] Annotation of /trunk/R/depmix.fitted.R
ViewVC logotype

Annotation of /trunk/R/depmix.fitted.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (view) (download)

1 : ingmarviss 10
2 :     #
3 :     # Ingmar Visser, 29-2-2008
4 :     #
5 :    
6 :     #
7 :     # DEPMIX.FITTED CLASS
8 :     #
9 :    
10 :     setClass("depmix.fitted",
11 :     representation(message="character", # convergence information
12 :     conMat="matrix", # constraint matrix on the parameters for general linear constraints
13 :     posterior="matrix" # posterior probabilities for the states
14 :     ),
15 :     contains="depmix"
16 :     )
17 :    
18 :     setMethod("fit",signature(object="depmix"),
19 : ingmarviss 64 function(object,w=NULL,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,...) {
20 : ingmarviss 10
21 :     # decide on method to be used
22 :     # em has to be added here under certain conditions
23 :     # in particular when there are linear constraints donlp should be used
24 : ingmarviss 64 # otherwise EM is good
25 : ingmarviss 10
26 : ingmarviss 64 # can/does EM deal with fixed constraints??? it should be possible for sure
27 :     if(is.null(method)) {
28 :     if(object@stationary&is.null(equal)&is.null(conrows)) {
29 :     method="EM"
30 :     } else {
31 :     method="donlp"
32 :     }
33 :     }
34 : ingmarviss 10
35 :     # determine which parameters are fixed
36 :     if(!is.null(fixed)) {
37 :     if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
38 :     } else {
39 :     if(!is.null(equal)) {
40 :     if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
41 :     fixed <- !pa2conr(equal)$free
42 :     } else {
43 :     fixed <- getpars(object,"fixed")
44 :     }
45 :     }
46 :     # set those fixed parameters in the appropriate submodels
47 :     object <- setpars(object,fixed,which="fixed")
48 :    
49 : ingmarviss 64 if(method=="EM") {
50 :     object <- em(object,verbose=TRUE)
51 : ingmarviss 10 }
52 :    
53 :     if(method=="donlp") {
54 : ingmarviss 64 # get the full set of parameters
55 :     allpars <- getpars(object)
56 :     # get the reduced set of parameters, ie the ones that will be optimized
57 :     pars <- allpars[!fixed]
58 : ingmarviss 10
59 : ingmarviss 64 # set bounds, if any
60 :     par.u <- rep(+Inf, length(pars))
61 :     par.l <- rep(-Inf, length(pars))
62 :    
63 :     # make loglike function that only depends on pars
64 :     logl <- function(pars) {
65 :     allpars[!fixed] <- pars
66 :     object <- setpars(object,allpars)
67 :     -logLik(object)
68 :     }
69 :    
70 : ingmarviss 10 if(!require(Rdonlp2)) stop("donlp optimization requires the 'Rdonlp2' package")
71 :    
72 :     # make constraint matrix and its upper and lower bounds
73 :     lincon <- matrix(0,nr=0,nc=npar(object))
74 :     lin.u <- numeric(0)
75 :     lin.l <- numeric(0)
76 :    
77 :     # incorporate equality constraints, if any
78 :     if(!is.null(equal)) {
79 :     if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
80 :     equal <- pa2conr(equal)$conr
81 :     lincon <- rbind(lincon,equal)
82 :     lin.u <- c(lin.u,rep(0,nrow(equal)))
83 :     lin.l <- c(lin.l,rep(0,nrow(equal)))
84 :     }
85 :    
86 :     # incorporate general linear constraints, if any
87 :     if(!is.null(conrows)) {
88 :     if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
89 :     lincon <- rbind(lincon,conrows)
90 :     if(conrows.upper==0) {
91 :     lin.u <- c(lin.u,rep(0,nrow(conrows)))
92 :     } else {
93 :     if(length(conrows.upper)!=nrow(conrows)) stop("'conrows.upper does not have correct length")
94 :     lin.u <- c(lin.u,conrows.upper)
95 :     }
96 :     if(conrows.lower==0) {
97 :     lin.l <- c(lin.l,rep(0,nrow(conrows)))
98 :     } else {
99 :     if(length(conrows.lower)!=nrow(conrows)) stop("'conrows.lower does not have correct length")
100 :     lin.l <- c(lin.l,conrows.lower)
101 :     }
102 :     }
103 :    
104 :     # select only those columns of the constraint matrix that correspond to non-fixed parameters
105 :     linconFull <- lincon
106 :     lincon <- lincon[,!fixed]
107 :    
108 :     # set donlp2 control parameters
109 :     cntrl <- donlp2.control(hessian=FALSE,difftype=2,report=TRUE)
110 :    
111 :     mycontrol <- function(info) {
112 :     return(TRUE)
113 :     }
114 :    
115 :     # optimize the parameters
116 :     result <- donlp2(pars,logl,
117 :     par.upper=par.u,
118 :     par.lower=par.l,
119 :     A=lincon,
120 :     lin.upper=lin.u,
121 :     lin.lower=lin.l,
122 :     control=cntrl,
123 :     control.fun=mycontrol
124 :     )
125 :    
126 :     class(object) <- "depmix.fitted"
127 :    
128 :     object@conMat <- linconFull
129 :     object@message <- result$message
130 :    
131 :     # put the result back into the model
132 :     allpars[!fixed] <- result$par
133 :     object <- setpars(object,allpars)
134 :    
135 :     }
136 :    
137 : ingmarviss 64 # FIX ME
138 :     object@posterior <- matrix(0,1,1)
139 :    
140 : ingmarviss 10 return(object)
141 :     }
142 :     )
143 :    
144 :     # convert conpat vector to rows of constraint matrix
145 :     pa2conr <- function(x) {
146 :     fix=as.logical(x)
147 :     x=replace(x,list=which(x==1),0)
148 :     un=setdiff(unique(x),0)
149 :     y=matrix(0,0,length(x))
150 :     for(i in un) {
151 :     z=which(x==i)
152 :     for(j in 2:length(z)) {
153 :     k=rep(0,length(x))
154 :     k[z[1]]=1
155 :     k[z[j]]=-1
156 :     y=rbind(y,k)
157 :     }
158 :     }
159 :     pa = list(free=fix,conr=y)
160 :     return(pa)
161 :     }
162 :    
163 :    
164 :     setMethod("show","depmix.fitted",
165 :     function(object) {
166 :     cat("Convergence info:",object@message,"\n")
167 : ingmarviss 64 print(logLik(object))
168 :     cat("AIC: ", AIC(object),"\n")
169 :     cat("BIC: ", AIC(object),"\n")
170 :     }
171 :     )
172 :    
173 :     setMethod("summary","depmix.fitted",
174 :     function(object) {
175 : ingmarviss 10 cat("Initial state probabilties model \n")
176 :     print(object@prior)
177 :     cat("\n")
178 :     for(i in 1:object@nstates) {
179 :     cat("Transition model for state (component)", i,"\n")
180 :     print(object@transition[[i]])
181 :     cat("\n")
182 :     }
183 :     cat("\n")
184 :     for(i in 1:object@nstates) {
185 :     cat("Response model(s) for state", i,"\n\n")
186 :     for(j in 1:object@nresp) {
187 :     cat("Response model for response",j,"\n")
188 :     print(object@response[[i]][[j]])
189 :     cat("\n")
190 :     }
191 :     cat("\n")
192 :     }
193 :     }
194 :     )

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge