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 18 - (view) (download)

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

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