I can not recover the original parameters

[genisprat]

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--")

[dmetivie]

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

HMMBase.jl/examples/basic_usage.jl

Line 24 in 8ad8254

z, y = rand(hmm, 500, seq = true)

I did not try to verify if it solves the problem.