Deepseek is a new open-source LLM that is good at Fortran programming, comparable to ChatGPT o1 I think.
Using the free chat interface to give the task
Write a Fortran subroutine thats a finite mixture of normals with a specified # of components to univariate data, using the EM algorithm, and write a main program that tests it for data simulated from a known mixture distribution.
and feeding compiler error messages back to Deepseek, I soon got the working program below. For production use you would probably tell the LLM to assume that common algorithms such as generating normal deviates or sorting are already available, for example in stdlib, so that the LLM does not reinvent the wheel.
module mixture_of_normals
implicit none
integer, parameter :: dp = kind(1.0d0) ! Double precision kind
real(dp), parameter :: pi = 3.141592653589793238462643383279502884197_dp
contains
! Subroutine to fit a finite mixture of normals using the EM algorithm
subroutine fit_mixture_of_normals(data, num_components, weights, means, variances, max_iter, tol, debug)
real(dp), intent(in) :: data(:) ! Input data
integer, intent(in) :: num_components ! Number of components
real(dp), intent(out) :: weights(:), means(:), variances(:) ! Output parameters
integer, intent(in) :: max_iter ! Maximum number of iterations
real(dp), intent(in) :: tol ! Convergence tolerance
logical, intent(in), optional :: debug ! Debugging flag
logical :: debug_ ! Local debugging flag
real(dp), allocatable :: responsibilities(:,:), log_likelihoods(:)
real(dp) :: log_likelihood, prev_log_likelihood
integer :: n, iter
! Initialize local debugging flag
debug_ = .false.
if (present(debug)) debug_ = debug
! Number of data points
n = size(data)
! Allocate memory for responsibilities and log-likelihoods
allocate(responsibilities(n, num_components), log_likelihoods(n))
! Initialize parameters
call initialize_parameters(data, num_components, weights, means, variances)
! EM algorithm
prev_log_likelihood = -huge(1.0_dp)
do iter = 1, max_iter
! E-step: Compute responsibilities
call compute_responsibilities(data, weights, means, variances, responsibilities, log_likelihoods)
! Compute log-likelihood
log_likelihood = sum(log_likelihoods)
! Check for convergence
if (abs(log_likelihood - prev_log_likelihood) < tol) exit
! M-step: Update parameters
call update_parameters(data, responsibilities, weights, means, variances)
! Debugging: Print current parameters
if (debug_) then
print '(a20, i0)', "Iteration:", iter
print '(a20, 100(f12.6, 1x))', "Weights:", weights
print '(a20, 100(f12.6, 1x))', "Means:", means
print '(a20, 100(f12.6, 1x))', "Variances:", variances
print '(a20, f12.6)', "Log-likelihood:", log_likelihood
print *, "----------------------------------------"
end if
prev_log_likelihood = log_likelihood
end do
deallocate(responsibilities, log_likelihoods)
end subroutine fit_mixture_of_normals
! Initialize parameters
subroutine initialize_parameters(data, num_components, weights, means, variances)
real(dp), intent(in) :: data(:)
integer, intent(in) :: num_components
real(dp), intent(out) :: weights(:), means(:), variances(:)
integer :: i, n ! Declare 'n' here
real(dp) :: data_min, data_max
! Number of data points
n = size(data)
data_min = minval(data)
data_max = maxval(data)
! Initialize weights uniformly
weights = 1.0_dp / num_components
! Initialize means and variances
do i = 1, num_components
means(i) = data_min + (data_max - data_min) * (i - 1) / (num_components - 1)
variances(i) = 1.0_dp
end do
end subroutine initialize_parameters
! Compute responsibilities (E-step)
subroutine compute_responsibilities(data, weights, means, variances, responsibilities, log_likelihoods)
real(dp), intent(in) :: data(:), weights(:), means(:), variances(:)
real(dp), intent(out) :: responsibilities(:,:), log_likelihoods(:)
integer :: n, num_components, i, j
real(dp), allocatable :: likelihood(:,:), total_likelihood(:)
n = size(data)
num_components = size(weights)
! Allocate memory for likelihood and total_likelihood
allocate(likelihood(n, num_components), total_likelihood(n))
! Compute likelihood for each component
do j = 1, num_components
do i = 1, n
likelihood(i, j) = weights(j) * normal_pdf(data(i), means(j), variances(j))
end do
end do
! Compute total likelihood for each data point
total_likelihood = sum(likelihood, dim=2)
! Compute responsibilities
do j = 1, num_components
responsibilities(:, j) = likelihood(:, j) / total_likelihood
end do
! Compute log-likelihood
log_likelihoods = log(total_likelihood)
deallocate(likelihood, total_likelihood)
end subroutine compute_responsibilities
! Update parameters (M-step)
subroutine update_parameters(data, responsibilities, weights, means, variances)
real(dp), intent(in) :: data(:), responsibilities(:,:)
real(dp), intent(inout) :: weights(:), means(:), variances(:)
integer :: n, num_components, j
real(dp), allocatable :: sum_responsibilities(:)
n = size(data)
num_components = size(weights)
! Allocate memory for sum_responsibilities
allocate(sum_responsibilities(num_components))
! Update weights
sum_responsibilities = sum(responsibilities, dim=1)
weights = sum_responsibilities / n
! Update means
do j = 1, num_components
means(j) = sum(data * responsibilities(:, j)) / sum_responsibilities(j)
end do
! Update variances
do j = 1, num_components
variances(j) = sum((data - means(j))**2 * responsibilities(:, j)) / sum_responsibilities(j)
end do
deallocate(sum_responsibilities)
end subroutine update_parameters
! Probability density function of a normal distribution
function normal_pdf(x, mean, variance) result(pdf)
real(dp), intent(in) :: x, mean, variance
real(dp) :: pdf
pdf = exp(-0.5_dp * (x - mean)**2 / variance) / sqrt(2.0_dp * pi * variance)
end function normal_pdf
end module mixture_of_normals
program test_mixture_of_normals
use mixture_of_normals
implicit none
integer, parameter :: n = 10**5 ! Number of data points
integer, parameter :: num_components = 2 ! Number of components
integer, parameter :: num_simulations = 5 ! Number of simulations
real(dp) :: data(n)
real(dp) :: weights(num_components), means(num_components), variances(num_components)
real(dp) :: true_weights(num_components), true_means(num_components), true_variances(num_components)
integer :: i, j, idx(num_components)
! True parameters of the mixture distribution
true_weights = [0.99_dp, 0.01_dp]
true_means = [0.0_dp, 0.0_dp]
true_variances = [1.0_dp, 100.0_dp]
print*,"#obs:", n ! line added by me
! Print true parameters once
print '(a20, a)', "True Parameters:", ""
print '(a20, 100(f12.6, 1x))', "Weights:", true_weights
print '(a20, 100(f12.6, 1x))', "Means:", true_means
print '(a20, 100(f12.6, 1x))', "Variances:", true_variances
print *, "========================================"
! Loop for multiple simulations and fittings
do i = 1, num_simulations
print '(a, i0)', "Simulation:", i
print *, "----------------------------------------"
! Simulate data from a known mixture distribution
call simulate_mixture_data(data, true_weights, true_means, true_variances)
! Fit mixture of normals to the data
call fit_mixture_of_normals(data, num_components, weights, means, variances, max_iter=100, tol=1.0e-6_dp, &
debug=.false.)
! Sort components by weight in descending order
do j = 1, num_components
idx(j) = j
end do
call sort_by_weight(weights, means, variances, idx, num_components)
! Print estimated parameters
print '(a20, a)', "Estimated Parameters:", ""
print '(a20, 100(f12.6, 1x))', "Weights:", weights
print '(a20, 100(f12.6, 1x))', "Means:", means
print '(a20, 100(f12.6, 1x))', "Variances:", variances
print *, "========================================"
end do
contains
! Simulate data from a known mixture distribution
subroutine simulate_mixture_data(data, true_weights, true_means, true_variances)
real(dp), intent(out) :: data(:)
real(dp), intent(in) :: true_weights(:), true_means(:), true_variances(:)
integer :: i
! Simulate data
do i = 1, n
call random_number(data(i))
if (data(i) < true_weights(1)) then
data(i) = true_means(1) + sqrt(true_variances(1)) * randn()
else
data(i) = true_means(2) + sqrt(true_variances(2)) * randn()
end if
end do
end subroutine simulate_mixture_data
! Random number generator for standard normal distribution
function randn() result(r)
real(dp) :: r
real(dp) :: u1, u2
call random_number(u1)
call random_number(u2)
r = sqrt(-2.0_dp * log(u1)) * cos(2.0_dp * pi * u2)
end function randn
! Sort components by weight in descending order
subroutine sort_by_weight(weights, means, variances, idx, num_components)
real(dp), intent(inout) :: weights(:), means(:), variances(:)
integer, intent(inout) :: idx(:)
integer, intent(in) :: num_components
integer :: i, j
real(dp) :: temp_weight, temp_mean, temp_variance
integer :: temp_idx
! Bubble sort (simple sorting algorithm)
do i = 1, num_components - 1
do j = i + 1, num_components
if (weights(i) < weights(j)) then
! Swap weights
temp_weight = weights(i)
weights(i) = weights(j)
weights(j) = temp_weight
! Swap means
temp_mean = means(i)
means(i) = means(j)
means(j) = temp_mean
! Swap variances
temp_variance = variances(i)
variances(i) = variances(j)
variances(j) = temp_variance
! Swap indices
temp_idx = idx(i)
idx(i) = idx(j)
idx(j) = temp_idx
end if
end do
end do
end subroutine sort_by_weight
end program test_mixture_of_normals
with sample output
#obs: 100000
True Parameters:
Weights: 0.990000 0.010000
Means: 0.000000 0.000000
Variances: 1.000000 100.000000
========================================
Simulation:1
----------------------------------------
Estimated Parameters
Weights: 0.989977 0.010023
Means: 0.002358 0.044858
Variances: 0.999245 98.262609
========================================
Simulation:2
----------------------------------------
Estimated Parameters
Weights: 0.989784 0.010216
Means: 0.004232 0.504152
Variances: 1.001481 101.687353
========================================
Simulation:3
----------------------------------------
Estimated Parameters
Weights: 0.989921 0.010079
Means: -0.001570 0.154525
Variances: 1.003607 107.735091
========================================
Simulation:4
----------------------------------------
Estimated Parameters
Weights: 0.989776 0.010224
Means: -0.001740 -0.078073
Variances: 0.993489 101.462980
========================================
Simulation:5
----------------------------------------
Estimated Parameters
Weights: 0.989377 0.010623
Means: 0.001496 0.394484
Variances: 0.995740 96.620077
========================================
I think we will soon see agents that use Deepseek to solve coding tasks, so that the user never sees code that does not compile.
ā I could have initially asked it to declare procedures as pure or elemental where possible. Asking it to do so for the code it generated, it declared