Nothing
##################################################################
#Functions for the parameterized m-step
#################################################################
get.rates<-function(param.values, covariate.array){
#############################################
#Author: JL, 11/15/2011
#gets the linear predictors corresponding to exp(theta)=exp(beta^TZ)
#INPUTS: covariate.array=array of dim(num.beta, num.theta, num.subjects)
# that for each subject and each transition contains the values of the covariates
# param.values =vector with parameter values
#OUTPUTS: matrix with rates for each subject dim(num.transitions.types, num.subjects)
##############################################################
linear.pred=apply(covariate.array*param.values,c(2,3),"sum")
return(exp(linear.pred))
}
rate_score<-function(rates,durations, deriv.array, transitions, transition.codes){
##############################################
#Author: JL, 11/15/2011
#Function to get the rate matrix score function
#INPUTS: rates=matrix with rates for each transition and each subject (dim(num.transitions, num.subjects))
# transition.codes= matrix with ni, and nj columns encoding the transitions
# corresponding to each rate
# transition= array of dim(num.states,num.states,num.subjects) with the transitions for each subject
# deriv.array = array of dim (num.betas,num.thetas, num.subjects) that contains the deriv. of theta WRT beta
# duration= matrix of dim(num,states, num.subjects) with the duration spent in each state
##############################################
suff.stat.score=array(NA,dim=c(dim(transition.codes)[1],dim(durations)[2]))
score=vector(length=c(dim(deriv.array)[1]))
for(l in 1:dim(transition.codes)[1]){
suff.stat.score[l,]=transitions[transition.codes[l,"ni"],transition.codes[l,"nj"],]-rates[l,]*durations[transition.codes[l,"ni"],]
}
for(m in 1:dim(durations)[2]){
score=score+deriv.array[,,m]%*%suff.stat.score[,m]
}
return(score)
}
##############################################
rate_hessian<-function(rates,durations, deriv.array, transition.codes){
##############################################
#Author: JL, 11/15/2011
#INPUTS: rates=matrix with rates for each transition and each subject (dim(num.transitions, num.subjects))
# deriv.array = array of dim (num.betas,num.thetas, num.subjects) that contains the deriv. of theta WRT beta
# transition.codes= matrix with ni, and nj columns encoding the transitions
# corresponding to each rate
# durations= matrix of dim(num,states, num.subjects) with the duration spent in each state
#OUTPUTS: hessian matrix for rate model parameters
##############################################
hessian=array(0,dim=c(dim(deriv.array)[1],dim(deriv.array)[1]))
if(dim(deriv.array)[1]==0){
return(hessian)
}else{
for(m in 1:dim(durations)[2]){
if(dim(deriv.array)[1]==1){
temp=deriv.array[,,m]
hessian=hessian+t(temp)%*%(temp*((-1)*rates[,m]*durations[transition.codes[,"ni"],m]))
}else{
hessian=hessian+deriv.array[,,m]%*%(t(deriv.array[,,m])*((-1)*rates[,m]*durations[transition.codes[,"ni"],m]))
}
}
return(hessian)
}
}
NR_update<-function(current.params,score, hessian){
##############################################
#Author: JL, 11/15/2011
#INPUTS: score= score vector
# hessian= hessian matrix
# current.params =vector of current values of the parameters
#OUTPUTS: next newton raphson update for parameters (vector of length num.parameters)
##############################################
return(current.params-solve(hessian)%*%score)
}
get.prob.array<-function(covariate.array, param.values){
#############################################
#get array of multinomial probabilities (not including the reference state) for all subects
##############################################
#Author: JL, 11/15/2011
#INPUTS:
# covariate_array: array of covariates for all subjects of dim(num.params,num.states,num.subjects)
# param.values: vector of values for model parameters (betas)
#OUTPUTS: array of dim(num.states,num.subects)
linear.pred=apply(covariate.array*param.values,c(2,3),"sum")
exp.linear.pred=exp(linear.pred)
return(t(t(exp.linear.pred)*1/(1+apply(exp.linear.pred,2,"sum"))))
}
multinom.score.old<-function(prob.matrix,deriv.array,counts,N_vector){
#############################################
#multinomial matrix score
##############################################
#Author: JL, 11/15/2011
#INPUTS: prob.matrix: matrix of dim(num.states, num.subjects) with probabilities calculated
# for each state/subject
# deriv.array = array of dim (num.betas,num.thetas, num.subjects) that contains the deriv. of theta WRT beta
# counts: array of category counts dim(num.states, num.sub) (excluding reference category)
# N_vector: vector (length =num.subjects) of multinomial sample size n for each subject
#OUTPUTS: score vector for multinomial parameters
num.subjects=dim(counts)[2]
suff.stat.score=array(0,dim=c(dim(counts)[1],dim(counts)[2]))
param.score=rep(0,times=dim(deriv.array)[1])
suff.stat.score=(counts-prob.matrix*N_vector)
for(m in 1:num.subjects){
param.score=param.score+matrix(deriv.array[,,m])%*%suff.stat.score[,m]
}
return(param.score)
}
multinom.score<-function(prob.matrix, deriv.array, counts, N_vector)
{
num.subjects = dim(counts)[2]
suff.stat.score = array(0, dim = c(dim(counts)[1], dim(counts)[2]))
param.score = rep(0, times = dim(deriv.array)[1])
suff.stat.score = (counts - prob.matrix * N_vector)
for (m in 1:num.subjects) {
param.score = param.score + matrix(deriv.array[, , m],ncol=dim(deriv.array)[2]) %*%
suff.stat.score[, m]
}
return(param.score)
}
multinom.covariance.array<-function(prob.matrix,N_vector){
##############################################################
#multinomial covariance array
############################################################
#Author: JL, 11/15/2011
#INPUTS: prob.matrix: matrix of dim(num.states, num.subjects) with probabilities calculated
# for each state/subject
# N_vector: vector (length =num.subjects) of multinomial sample size n for each subject
#OUTPUT array (nums.state,num.states,num.subjects) of multinomial variance-covariance matrices
out=array(0,dim=c(dim(prob.matrix)[1],dim(prob.matrix)[1],dim(prob.matrix)[2]))
for(m in 1:dim(prob.matrix)[2]){
if(dim(out)[1]==1){
out[,,m]=prob.matrix[,m]*(1-prob.matrix[,m])
}else{
out[,,m]=(-1)*prob.matrix[,m]%*%t(prob.matrix[,m])
diag(out[,,m])=prob.matrix[,m]*(1-prob.matrix[,m])
}
out[,,m]=out[,,m]*N_vector[m]
}
return(out)
}
multinom.hessian<-function(prob.matrix, deriv.array,N_vector){
#############################################
#multinomial matrix hessian
##############################################
#Author: JL, 11/15/2011
#INPUTS: prob.matrix: matrix of dim(num.states, num.subjects) with probabilities calculated
# for each subject
# deriv.array = array of dim (num.betas,num.thetas, num.subjects)
# that contains the deriv. of theta WRT beta
# N_vector: vector (length =num.subjects) of multinomial sample size n for each subject
#OUTPUTS: array of hessian matrix for multinomial parameters (dim num.param, num.params, num.subjects)
#covariance_array=array(0,dim=c(dim(prob.matrix)[1],dim(prob.matrix)[1],dim(prob.matrix)[2]))
#covariance_array=multinom.covariance.array(prob.matrix,N_vector)
#hessian=array(0,dim=c(dim(deriv.array)[1],dim(deriv.array)[1]))
#for(m in 1:dim(deriv.array)[3]){
# hessian=hessian+deriv.array[,,m]%*%covariance_array[,,m]%*%t(deriv.array[,,m])
#}
#return(-1*hessian)
myDims=dim(deriv.array)
out=.Call("multinom_hessian",prob.matrix,deriv.array,myDims,N_vector)
return(out)
}
emission.score.hessian.old<-function(emission_updates,param.values,emission,deriv.array,no.NR=0){
#######################################################################
#emission score and hessian updates
##############################################
#Author: JL, 11/20/2011
#INPUTS:
# emission_updates: matrix (dim(num.hidden.states,num.obs.states,num.subjects)) with the emission count data
# param.values: vector of values for model parameters
# emission: emission object with covariate array, param.types
# deriv.array = array of dim (num.betas,num.thetas, num.subjects)
# that contains the deriv. of theta WRT beta
# no.NR (indicator =1 if not doing NR update )
#OUTPUTS: a list with two elements: score (vector of length(num.params)) and hessian (matrix of dim(num.params,num.params))
##########################################
if(!no.NR){
N_array=apply(emission_updates,c(1,3),"sum")
}
score=rep(0,times=dim(deriv.array)[1])
hessian=array(0,dim=c(dim(deriv.array)[1],dim(deriv.array)[1]))
hidden_states=unique(emission$emission.states["i"])
all.states=c(0,0)
prob.array=rep(0,length=dim(emission$covariate.array)[3])
for(k in 1:dim(hidden_states)[1]){
states=emission$emission.states[emission$emission.states$i==hidden_states[k,"i"],]
cov.array.states=array(emission$covariate.array[,as.numeric(rownames(states)),],dim=c(dim(emission$covariate.array)[1],dim(states)[1],dim(emission$covariate.array)[3]))
deriv.array.states=array(deriv.array[,as.numeric(rownames(states)),],dim=c(dim(deriv.array)[1],dim(states)[1],dim(deriv.array)[3]))
prob.matrix=get.prob.array(cov.array.states, param.values)
if(!no.NR){
N_vector=N_array[hidden_states[k,"i"],]
counts=matrix(emission_updates[states$i,states$j,],ncol=dim(emission_updates)[3])
score=score+multinom.score(prob.matrix,deriv.array.states,counts,N_vector)
hessian=hessian+multinom.hessian(prob.matrix,deriv.array.states,N_vector)
}
all.states=rbind(all.states,states)
prob.array=rbind(prob.array,prob.matrix)
}
rownames(all.states)=seq(0,(dim(all.states)[1]-1))
return(list(score=score,hessian=hessian,states=all.states[-1,],prob.array=prob.array[-1,]))
}
emission.score.hessian=function (emission_updates, param.values, emission, deriv.array, no.NR = 0)
{
if (!no.NR) {
N_array = apply(emission_updates, c(1, 3), "sum")
}
score = rep(0, times = dim(deriv.array)[1])
hessian = array(0, dim = c(dim(deriv.array)[1], dim(deriv.array)[1]))
hidden_states = unique(emission$emission.states["i"])
all.states = c(0, 0)
prob.array = rep(0, length = dim(emission$covariate.array)[3])
for (k in 1:dim(hidden_states)[1]) {
states = emission$emission.states[emission$emission.states$i ==
hidden_states[k, "i"], ]
cov.array.states = array(emission$covariate.array[, as.numeric(rownames(states)),
], dim = c(dim(emission$covariate.array)[1], dim(states)[1],
dim(emission$covariate.array)[3]))
deriv.array.states = array(deriv.array[, as.numeric(rownames(states)),
], dim = c(dim(deriv.array)[1], dim(states)[1], dim(deriv.array)[3]))
prob.matrix = get.prob.array(cov.array.states, param.values)
if (!no.NR) {
N_vector = N_array[hidden_states[k, "i"], ]
counts=matrix(ncol=dim(emission_updates)[3],nrow=dim(states)[1])
for(l in 1:dim(states)[1]){
counts[l,] = matrix(emission_updates[states$i[l], states$j[l],], ncol = dim(emission_updates)[3])
}
score = score + multinom.score(prob.matrix, deriv.array.states,
counts, N_vector)
hessian = hessian + multinom.hessian(prob.matrix,
deriv.array.states, N_vector)
}
all.states = rbind(all.states, states)
prob.array = rbind(prob.array, prob.matrix)
}
rownames(all.states) = seq(0, (dim(all.states)[1] - 1))
return(list(score = score, hessian = hessian, states = all.states[-1,
], prob.array = prob.array[-1, ]))
}
update.multinom.params<-function(current.params,covariate.array,deriv.array,param.types,counts=NULL,N_vector=NULL,max.it,no.NR=0){
#############################################
#update the multinomial parameters
##############################################
#Author: JL, 11/15/2011
#INPUTS: current.params: initial values of the parameters
# covariate_array: array of covariates for all subjects of dim(num.states,num.states,num.subjects)
# deriv_array: array of covariates for all subjects corresponding to non-fixed parameters
# param.types: vector for each parameter coding NR (0), simple update (1), or fixed (2)
# counts: array of category counts dim(num.states, num.sub) (excluding reference category)
# N_vector: vector (length =num.subjects) of multinomial sample size n for each subject
# max.it=maximum number of NR updates
# no.NR (indicator =1 if not doing NR update )
#OUTPUTS: updated values of the parameters
#######################################################
#number of parameters for NR update
num.NR.params=dim(deriv.array)[1]
#deriv array for NR update
########CHANGE
#deriv.array=array(covariate.array[param.types==0,,],dim=c(sum(param.types==0),dim(covariate.array)[2],dim(covariate.array)[3]))
hessian=NULL
score=array(0,dim=c(num.NR.params,(max.it+1)))
NR.params=array(0,dim=c(num.NR.params,(max.it+1)))
NR.params[,1]=current.params[param.types==0]
#score=array(0,dim=c(length(current.params),(max.it+1)))
params=array(0,dim=c(length(current.params),(max.it+1)))
params[,1]=current.params
for(i in 1:(max.it)){
prob.matrix=get.prob.array(covariate.array,params[,i])
if(!no.NR){
score[,i]=multinom.score(prob.matrix,deriv.array,counts,N_vector)
hessian=multinom.hessian(prob.matrix,deriv.array,N_vector)
NR.params[,i+1]=NR_update(current.params=NR.params[,i],score[,i], hessian)
params[,i+1]=current.params
params[param.types==0,i+1]<-NR.params[,i+1]
}
}
return(list(params=params,prob.matrix=prob.matrix,hessian=hessian))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.