# Understanding mutating functions correctly

I am taking the next step with Julia, and trying to learn more. Up to now, I have used it mostly just as a calculator, but I recently started to make my own modules and so on.

One thing I am having trouble understanding is how mutating functions work.

In the code below, I send in a mutable struct, GLMFit, and want to set some of its fields/attributes with this function (full code at the end).

If I do not include this part after the loop is finished, the \beta-, \mu-, and \eta-vectors have the same value as before the loop.

        M.Fit.β = β
M.Fit.μ = μ
M.Fit.η = η


But the \beta is still changed inside the function, before the loop, since I set the first element equal to the mean of y, with this code:

if sum(β) == 0
β[1] = sum(y)/length(y)
end


It is just the changes to \beta that happen inside the loop which are not captured, unless I specifically include the M.Fit.β = β after the loop.

Am I doing something wrong, or is the M.Fit.β = β and same for the other vectors strictly necessary in these cases?

The full code for the function:

function FisherScoring!(M::GLMFit, maxiter=25, tol=1e-10)
dist, lp, link = M.Model.d, M.Model.lp, M.Model.g
X = lp.X
β = M.Fit.β
μ = M.Fit.μ
y = M.Fit.y
η = M.Fit.η
p = length(β)

local J, Jinv, w

if sum(β) == 0
β[1] = sum(y)/length(y)
end

for i in 1:maxiter

η = X * β
w = Var(dist, μ)

U = Score(X, y, μ, dist, link)
J = FisherInformation(X, w, dist, link)
Jinv = inv(J)

βtemp = β

β = β + Jinv * U
δ = β - βtemp
if sqrt(δ' *  δ) < tol
println(@sprintf "Stopping after %d iterations" i)
M.Fit.SE = sqrt.(diag(Jinv))
break
else
L = Likelihood(dist, μ, y)
println(@sprintf "Fisher Iteration %d 𝐿(β) = %.12f" i L)
M.Fit.SE = sqrt.(diag(Jinv))
end
end
M.Fit.β = β
M.Fit.μ = μ
M.Fit.η = η
M.Fit.𝐃 = Deviance(dist, μ, y)
M.Fit.DoF = length(y) - length(β)
L = Likelihood(dist, μ, y)
M.Fit.AIC = AIC(length(β), L)
M.IsFitted = true
end


This depends on what β is — if it is an array, then β = M.Fit.β should be sufficient to make sure they are the same, except that you create new one with

β = β + Jinv * U


(In the future, please include a self-contained MWE, which should make it much easier to help you.)

2 Likes

It is an array. I totally forgot to use the .= syntax to fill it instead of overwriting it.

βtemp = β
β = β + Jinv * U
δ = β - βtemp


If you mutate β, then δ will always be a zero vector, a better way would be

δ = Jinv * U
β .+= δ


Another thing to consider is that inv(A) * b is usually less numerically stable than A \ b (the latter means "solve Ax = b for x", i.e. it is mathematically equivalent to the former).

3 Likes