Understanding mutating functions correctly

New to Julia
1

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 * β
            μ = InverseLinkFunction(link, η)
            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
2

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
3

Thank you for your help!

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

4

Beware, though, about the lines:

β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