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 https://github.com/maxmouchet/HMMBase.jl/blob/8ad8254d95f30c56dbff874fb6e8040969849e05/examples/basic_usage.jl#L24 I did not try to verify if it solves the problem.