Some handy utilities for calculating and drawing linear regressions in Julia, including residual analysis.

linreg.jl

import LsqFit    
import Plots    
    
function fwhm(sigma::Real)::Float64    
    # Calculate FWHM from sigma.    
    2sqrt(2log(2))*sigma    
end    
    
const M64 = Measurements.Measurement{Float64}      
                                  
                    
mfl = (v,e) -> Measurements.measurement(v, e)      
mean = a -> sum(a)/length(a)      
         
function nv(x)      
    Measurements.value.(x)            
end                                   
         
function sd(x)      
    Measurements.uncertainty.(x)      
end                                  
                                     
struct LinReg                
    m::Measurements.Measurement      
    b::Measurements.Measurement      
    chisq::Float64      
    samplesize::Float64                                      
    chisq_dof::Float64                                                 
end                                                                    
                                                                 
function linreg(xs::Vector{<:Real}, ys::Vector{<:Real})          
    m = (mean(xs.*ys)-mean(xs)*mean(ys))/(mean(xs.^2)-mean(xs)^2)      
    b = mean(ys)-m*mean(xs)       
    chisq = sum( (nv(ys) .- (nv(m) .*nv(xs) .+ nv(b))).^2 /      
                (sd(ys).^2 .+ (nv(m) .* sd(xs)).^2) )          
    samplesize = length(xs)      
    chisq_dof = samplesize-2                                  
    return LinReg(m, b, nv(chisq), samplesize, chisq_dof)      
end                                     
                                        
function linreg_residual_plot(xs, ys; figsize=(600,350))          
    lr = linreg(xs, ys)                                           
    fitys = nv(lr.m .* xs .+ lr.b)        
    resids = nv(ys .- fitys)                               
    residerrs = sqrt.((sd(ys)).^2 + (nv(lr.m) .* sd(xs)).^2)         
    resids = mfl.(resids, residerrs)                                                
                                                                                                                                       
    lay = Plots.grid(2, 1, heights=[0.7, 0.3])                                                                                         
                                                                                                                                       
    Plots.plot(xs, ys, t=:scatter,                                                   
               label="Daten, χ²/ndof = $(round(lr.chisq, digits=1))/$(lr.chisq_dof) = $(round(lr.chisq/lr.chisq_dof, digits=1))",      
               legend=:top, xlabel="t / µs", ylabel="s / mm")                                                                          
    dataplot = Plots.plot!(xs, fitys, label="mx + b = ($(lr.m)) x + ($(lr.b))")      
                                                                                                          
    Plots.plot(nv(xs), resids, legend=false, t=:scatter,                    
               ylabel="Daten - (mx + b)")                                                                   
    residplot = Plots.plot!(nv(xs), zeros(size(xs)), color="gray")      
    Plots.plot(dataplot, residplot, layout=lay, size=figsize, link=:x)                                                                                     
end