# I made a EM Algorithm for Gaussian Mixture Models, but that is not working

**URL:** <https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205>\
**Category:** General Usage\
**Tags:** question\
**Created:** [February 13, 2021, 11:37am UTC](https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205 "2021-02-13T11:37:00Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![ryuma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryuma/32/21949_2.png) [@ryuma](https://discourse.julialang.org/u/ryuma)\
**Post date:** [February 13, 2021, 11:37am UTC](https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205/1 "2021-02-13T11:37:00Z")

</div>

I made a EM Algorithm for Gaussian Mixture Models, but that is not working.

error massages → matrix is not Hermitian; Cholesky factorization failed.

why error happen, please teach me.

```julia
using LinearAlgebra
using Distributions

struct GMMsResults
    post_prob
    mean
    sigma
    mixing_coef
end

function GMMs!( data, num_class ;seed=1, max_iter=100 )
    
    num_data, dims = size( data )
    rng = MersenneTwister( seed )
    
    # 初期化
    post_prob = fill( 1/num_class, num_data, num_class ) 
    mean = rand( rng, dims, num_class )
    sigma = [Matrix{Float64}( I, dims, dims ) for _ in 1:num_class] # 配列にΣがクラス数分入っている
    mixing_coef = fill( 1/num_class, num_class )
    
    for _ in 1:max_iter

        update_post_prob!( data, num_class, post_prob, mean, sigma, mixing_coef )
        sum_class_assignment = sum( post_prob', dims=2 )
        update_mean!( data, num_class, post_prob, mean, sum_class_assignment )
        update_sigma!( data, num_class, post_prob, mean, sigma, sum_class_assignment )
        update_mixing_coef!( data, num_class, mixing_coef, sum_class_assignment )
        
    end
    
    return GMMsResults( post_prob, mean, sigma, mixing_coef )
end

function update_post_prob!( data, num_class, post_prob, mean, sigma, mixing_coef )
    
    num_data = size( data )[1]
    
    for n in 1:num_data
        sum_tmp = 0
        for k in 1:num_class
            sum_tmp += mixing_coef[k] * pdf( MvNormal( mean[:,k], sigma[k] ), data[n,:] )
        end
        
        for k in 1:num_class
            post_prob[n,k] = ( mixing_coef[k] * pdf( MvNormal( mean[:,k], sigma[k] ), data[n,:] ) ) / sum_tmp
        end
    end
    
end

function update_mean!( data, num_class, post_prob, mean, sum_class_assignment )
    
    num_data = size( data )[1]
    
    for k in 1:num_class
        sum_tmp = zeros( num_class )
        for n in 1:num_data
            sum_tmp += post_prob[n,k] .* data[n,:]
        end    

        mean[:,k] = sum_tmp ./ sum_class_assignment[k]
    end
    
end

function update_sigma!( data, num_class, post_prob, mean, sigma, sum_class_assignment )
    
    num_data, dims = size( data )

    for k in 1:num_class
        sum_tmp = zeros( dims, dims )
        for n in 1:num_data
            sum_tmp += post_prob[n,k] .* ( data[n,:] - mean[:, k] ) * ( data[n,:] - mean[:, k] )'
        end    

        sigma[k] = sum_tmp ./ sum_class_assignment[k]
    end
    
end

function update_mixing_coef!( data, num_class, mixing_coef, sum_class_assignment )
    
    num_data = size( data )[1]
    
    for k in 1:num_class
        mixing_coef[k] = sum_class_assignment[k] / num_data
    end
    
end

```

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [February 13, 2021, 9:36pm UTC](https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205/2 "2021-02-13T21:36:39Z")

</div>

Does it happens with any data or just in specific cases ?

You can compare your algo with this one: [Clustering · BetaML.jl Documentation](https://sylvaticus.github.io/BetaML.jl/dev/Clustering.html#BetaML.Clustering.gmm-Tuple%7BAny,Any%7D)

If I remember correctly I had the same problem when I didn’t implemented a minimum variance/covariance in the Gaussian…

---

<div class="post-metadata">

**Author:** ![ryuma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryuma/32/21949_2.png) [@ryuma](https://discourse.julialang.org/u/ryuma)\
**Post date:** [February 14, 2021, 2:32am UTC](https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205/3 "2021-02-14T02:32:13Z")

</div>

Thanks for reply.  
I used any data. data is generated by multivariate normal distribution.

e.g.

```julia
μ = ones( 2 )
Σ = I
x = MvNormal( μ, Σ );
y = MvNormal( μ .+ 5 , Σ + ones( 2,2 ) );

X = rand( x, 50 ) 
Y = rand( y, 50 ) 
data = hcat( X, Y )';

plot( aspect_ratio = 1 )
scatter!( data[:,1], data[:,2], markersize=8 )
covellipse!( x.μ, x.Σ )
covellipse!( y.μ, y.Σ )

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/5/75f7f205e30fdba1fa745a147b81a69c35dc471d.png)

I`ll compare my alogo with that one.

---

<div class="post-metadata">

**Author:** ![ryuma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryuma/32/21949_2.png) [@ryuma](https://discourse.julialang.org/u/ryuma)\
**Post date:** [February 14, 2021, 8:03am UTC](https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205/4 "2021-02-14T08:03:31Z")

</div>

I solved this problem. A cause of the error is MvNormal in Distribution.jl . When I gave MvNormal positive-definite matrix, MvNormal returned error. So, I replaced MvNormal with gaussian I implemented. My algo is correctly working.

---

<div class="post-metadata">

**Author:** ![ryuma](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryuma/32/21949_2.png) [@ryuma](https://discourse.julialang.org/u/ryuma)\
**Post date:** [February 17, 2021, 11:18am UTC](https://discourse.julialang.org/t/i-made-a-em-algorithm-for-gaussian-mixture-models-but-that-is-not-working/55205/5 "2021-02-17T11:18:47Z")

</div>

Another solution is using a type::Symmetric. When my algo calculate a variance/covariance, it caused a numeric error on Float. That’s why MvNormal() returned error.