HMMBase.jl
HMMBase.jl copied to clipboard
I can not recover the original parameters
Hi, First of all, thank you for this nice package. I would like to use this package in a research project but I can not recover the original parameters when I use fit_mle. I create synthetic data using the function rand() and then I try to fit it with different initial conditions. For a HMM with bernulli distribution. The Transition matrix is correctly fit it but the fitted emission probabilities aare close to 1. Probably I am doing something wrong but I can not find my error. Here my code:
P=[0.7 0.3 0.1 0.9]
T=[0.7 0.3 0.1 0.9]
NDataSets=100 Nconditions=20 Nstates=2 Nout=2 Ntrials=1000 LlOriginal=zeros(NDataSets) LlBest=zeros(NDataSets)
TBest=zeros(NDataSets,Nstates,Nstates) PiBest=zeros(NDataSets,Nstates) PBest=zeros(NDataSets,Nstates,Nout) TFinal=zeros(Nconditions,Nstates,Nstates) PiFinal=zeros(Nconditions,Nstates) PFinal=zeros(Nconditions,Nstates,Nout) Ll=zeros(Nconditions)
for idata in 1:NDataSets println("iDataSet: ",idata)
#create a new hmm object
hmm = HMM(T[:,:], [Categorical(P[1,:]), Categorical(P[2,:])])
println(hmm.A)
println(hmm.B[1].p)
println(hmm.B[2].p)
choice,state=rand(hmm, Ntrials, seq = true) #create synthetic data
LlOriginal[idata]=-loglikelihood(hmm,choice) #LL of the original parameters
for icondition in 1:Nconditions
#randomize initial condition
aux=rand(4)
hmm.A[1,1]=aux[1]
hmm.A[1,2]=1-aux[1]
hmm.A[2,1]=aux[2]
hmm.A[2,2]=1-aux[2]
hmm.B[1].p[1]=aux[3]
hmm.B[1].p[2]=1-aux[3]
hmm.B[2].p[1]=aux[4]
hmm.B[2].p[2]=1-aux[4]
hmm2,history=fit_mle(hmm, choice) #fit
#println(hmm2.A)
TFinal[icondition,:,:]=hmm2.A
PFinal[icondition,1,:]=hmm2.B[1].p
PFinal[icondition,2,:]=hmm2.B[2].p
PiFinal[icondition,:]=hmm2.a
Ll[icondition]=-loglikelihood(hmm2,choice)
end
#best parameters fit of the Nconditions initial conditions
imin=findall(x->x==minimum(Ll),Ll)[1]
TBest[idata,:,:]=TFinal[imin,:,:]
PBest[idata,:,:]=PFinal[imin,:,:]
PiBest[idata,:,:]=PiFinal[imin,:,:]
LlBest[idata]=Ll[imin]
end
plot([T[1,1],T[1,1]],[1,3],"r-") plot([T[2,2],T[2,2]],[1,3],"r-")
plot([P[1,1],P[1,1]],[3,5],"b--") plot([P[2,2],P[2,2]],[3,5],"b--")
I think that in
choice,state=rand(hmm, Ntrials, seq = true) #create synthetic data
it should be
state,choice=rand(hmm, Ntrials, seq = true) #create synthetic data
where choice
are the observations and state
the unobserved data. Next line you compute the likelihood with states instead of observations.
See
https://github.com/maxmouchet/HMMBase.jl/blob/8ad8254d95f30c56dbff874fb6e8040969849e05/examples/basic_usage.jl#L24
I did not try to verify if it solves the problem.