s-baumann / multivariatefunctions.jl

A framework for multivariate functions together with constructors for schumaker splines, OLS, Chebyshev, MARS splines for approximation.

MultivariateFunctions.jl

Build	Coverage	Documentation

This implements single algebra and evaluation on Multivariate functions. There are a few ways in which it can be used.

• This can be used for approximation functions. It can currently implement OLS functions, Chebyshev polynomials, the Schumaker shape preserving spline and basic interpolation schemes. It can also do Recursive Partitioning and create Multivariate Adaptive Regression (MARS) Splines. It could be extended to implement other approximation schemes.
• All basic algebra and calculus on a MultivariateFunction can be done analytically.
• Newton's method is implemented so that roots and optima can be found using analytical Jacobians and Hessians.

This is unlikely to be updated again. I didn't really understand Julia well when I wrote it. I think it is all correct but it is pretty inefficient. So use with care.

multivariatefunctions.jl's Issues

Broken examples of MARS

Issue

MARS example seems broken.

julia> using MultivariateFunctions


julia> using Random


julia> using DataFrames


julia> using Distributions


julia> using DataStructures


julia> Random.seed!(1992)

MersenneTwister(1992)

julia> nObs = 1000

1000

julia> dd = DataFrame()

0×0 DataFrame


julia> dd[:x] = rand( Normal(),nObs) + 0.1 .* rand( Normal(),nObs)

1000-element Vector{Float64}:
 -1.5194924275498334
  0.5844430822023597
 -0.5497375840875416
  1.2817822410251947
 -0.19060992619057496
 -0.47048055970036007
  0.6818583626457018
 -0.8554162326299839
 -1.3695617445991064
  1.0953038238849109
 -1.699099302747622
 -0.6780117585650899
  1.4192925709333952
  0.42807807496395617
  1.0352222109740405
  ⋮
 -0.10688015795519358
  0.9430607035612268
 -0.20183783555815388
  0.6460559374899443
  0.026095413811565132
 -0.6833593785111719
  0.07431322076156123
  0.2014010440942761
 -0.7793609005754738
  0.39473020003581993
 -0.5449445196399865
 -1.0193446965787274
 -1.2568550834023509
 -1.178925006150981
 -1.7098961465470663

julia> dd[:z] = rand( Normal(),nObs) + 0.1 .* rand( Normal(),nObs)

1000-element Vector{Float64}:
 -1.1750705521302567
  0.9069556021870959
  0.20960467727966883
 -0.016935093159266515
  0.32734357936941416
 -0.4583427779582831
 -1.2525094667318868
 -0.8855902909295208
 -0.7514670213109499
  0.37216278394183117
 -0.28668599160737784
  0.685074075486414
  0.3184272652769276
 -0.5276955662875091
  0.38991220787275427
  ⋮
  1.3861777686416614
  1.0765345340164898
 -0.1652608047507846
  0.719043828393897
  0.8701332884169406
  0.767772517601333
  1.0669198445676553
  0.11335691897963737
  1.5910107187781224
 -0.7128766324283243
 -0.7931006582401271
  0.4976775640182734
 -0.6536206642622754
 -0.31287850605786915
 -0.7850092399544946

julia> dd[:w] = (0.5 .* rand( Normal(),nObs)) .+ 0.7.*(dd[:z] .- dd[:x]) + 0.1 .* rand( Normal(),nObs)

1000-element Vector{Float64}:
 -0.0019189980122132633
  0.4725494429492194
 -0.0675946792898773
 -1.157130519654445
  0.14522767354377503
 -0.005382762866542205
 -1.2005096892042109
  0.13085627620806112
  0.05612174029581857
 -0.6172526046137324
  1.4778484610424196
  1.054709698408579
 -1.2271058467494558
 -0.5353888688122925
 -0.7294356191124297
  ⋮
  0.615069475148343
  0.2181079071913068
  0.09486662765477802
  0.4697963855574916
  0.15935631987862492
  0.5755374838341898
  0.7882977393176024
 -0.2396016872530803
  1.5271293107037496
 -0.8947695470014021
 -1.072537042317404
  0.04893462310292454
  0.026454764346978646
  0.5885500404818271
  1.3937893199788063

julia> dd[:y] = (dd[:x] .*dd[:w] ) .* (dd[:z] .- dd[:w]) .+ dd[:x] + rand( Normal(),nObs)

1000-element Vector{Float64}:
 -2.3487840899451675
  1.0087441497559786
 -2.280470624644659
 -0.9242658326625207
  0.815181257773693
 -0.5031071711999816
  0.30192536890203986
 -2.035262792736697
  0.03194223795679596
  1.13595311468184
  3.7442843912932613
  0.6009073972935214
 -0.44822283045481404
  0.6405293235510074
 -0.16269142011610555
  ⋮
 -2.4459349487643323
  1.7869773765861698
 -2.046579053983459
  0.9006917341756171
  0.41441684404908374
 -0.010334525049836674
  1.0941056908369529
  1.4286676478281712
 -0.26487929690458845
  0.5752509993328735
 -0.5676424153457388
 -1.9430060207426862
 -2.6246440845770023
 -1.5385613151541713
  4.127923473876518

julia> dd[7,:y] = 1.0

1.0

julia> y = :y

:y

julia> x_variables = Set{Symbol}([:w, :x, :z])
Set{Symbol} with 3 elements:
  :w
  :z
  :x

julia> number_of_divisions = 7
7

julia> rp_4, rp_reg_4 = create_recursive_partitioning(dd, y, x_variables, number_of_divisions; rel_tol = 1e-3)

ERROR: MethodError: Cannot `convert` an object of type String to an object of type Symbol
Closest candidates are:
  convert(::Type{S}, ::T) where {S, T<:CategoricalValue} at /Users/jinrae/.julia/packages/CategoricalArrays/0ZAbp/src/value.jl:68
  convert(::Type{T}, ::T) where T at essentials.jl:205
  Symbol(::String) at boot.jl:478
  ...
Stacktrace:
  [1] setindex!(h::Dict{Symbol, Float64}, v0::Float64, key0::String)
    @ Base ./dict.jl:374
  [2] evaluate(f::Piecewise_Function, coordinates::DataFrame)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/0_structs_and_generic_reversals.jl:717
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [5] getindex
    @ ./broadcast.jl:575 [inlined]
  [6] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [7] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [8] copyto!(dest::Vector{Vector{Union{Float64, MultivariateFunction}}}, bc::Base.Broadcast.Broadcasted{Nothing, Tuple{Base.OneTo{Int64}}, typeof(evaluate), Tuple{Vector{Piecewise_Function}, Base.RefValue{DataFrame}}})
    @ Base.Broadcast ./broadcast.jl:983
  [9] copyto!
    @ ./broadcast.jl:936 [inlined]
 [10] copy
    @ ./broadcast.jl:908 [inlined]
 [11] materialize
    @ ./broadcast.jl:883 [inlined]
 [12] create_ols_approximation(dd::DataFrame, y::Symbol, model::Vector{Piecewise_Function}; allowrankdeficient::Bool)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/5_ols_regression.jl:43
 [13] create_ols_approximation
    @ ~/.julia/packages/MultivariateFunctions/a2IXB/src/5_ols_regression.jl:43 [inlined]
 [14] optimise_given_specific_split(dd::DataFrame, y::Symbol, array_of_funcs::Vector{Piecewise_Function}, ind::Int64, split_variable::Symbol, split_point::Float64, SplitFunction::typeof(MultivariateFunctions.add_split_with_step_function), removeSplitFunction::Bool)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:28
 [15] #56
    @ ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:36 [inlined]
 [16] optimize(f::MultivariateFunctions.var"#56#57"{DataFrame, Symbol, Vector{Piecewise_Function}, Int64, Symbol, typeof(MultivariateFunctions.add_split_with_step_function), Bool}, x_lower::Float64, x_upper::Float64, mo::Optim.Brent; rel_tol::Float64, abs_tol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, callback::Nothing, show_every::Int64, extended_trace::Bool)
    @ Optim ~/.julia/packages/Optim/TNmSw/src/univariate/solvers/brent.jl:69
 [17] #optimize#80
    @ ~/.julia/packages/Optim/TNmSw/src/univariate/optimize/interface.jl:21 [inlined]
 [18] optimise_split(dd::DataFrame, y::Symbol, array_of_funcs::Vector{Piecewise_Function}, ind::Int64, split_variable::Symbol, rel_tol::Float64, SplitFunction::Function, removeSplitFunction::Bool)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:36
 [19] create_recursive_partitioning(dd::DataFrame, y::Symbol, x_variables::Set{Symbol}, MaxM::Int64; rel_tol::Float64)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:65
 [20] top-level scope
    @ REPL[17]:1

julia>

julia> rp_1, rp_reg_1 = create_mars_spline(dd, y, x_variables, number_of_divisions; rel_tol = 1e-3)
ERROR: MethodError: Cannot `convert` an object of type String to an object of type Symbol
Closest candidates are:
  convert(::Type{S}, ::T) where {S, T<:CategoricalValue} at /Users/jinrae/.julia/packages/CategoricalArrays/0ZAbp/src/value.jl:68
  convert(::Type{T}, ::T) where T at essentials.jl:205
  Symbol(::String) at boot.jl:478
  ...
Stacktrace:
  [1] setindex!(h::Dict{Symbol, Float64}, v0::Float64, key0::String)
    @ Base ./dict.jl:374
  [2] evaluate(f::Piecewise_Function, coordinates::DataFrame)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/0_structs_and_generic_reversals.jl:717
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [5] getindex
    @ ./broadcast.jl:575 [inlined]
  [6] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [7] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [8] copyto!(dest::Vector{Vector{Union{Float64, MultivariateFunction}}}, bc::Base.Broadcast.Broadcasted{Nothing, Tuple{Base.OneTo{Int64}}, typeof(evaluate), Tuple{Vector{Piecewise_Function}, Base.RefValue{DataFrame}}})
    @ Base.Broadcast ./broadcast.jl:983
  [9] copyto!
    @ ./broadcast.jl:936 [inlined]
 [10] copy
    @ ./broadcast.jl:908 [inlined]
 [11] materialize
    @ ./broadcast.jl:883 [inlined]
 [12] create_ols_approximation(dd::DataFrame, y::Symbol, model::Vector{Piecewise_Function}; allowrankdeficient::Bool)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/5_ols_regression.jl:43
 [13] create_ols_approximation
    @ ~/.julia/packages/MultivariateFunctions/a2IXB/src/5_ols_regression.jl:43 [inlined]
 [14] optimise_given_specific_split(dd::DataFrame, y::Symbol, array_of_funcs::Vector{Piecewise_Function}, ind::Int64, split_variable::Symbol, split_point::Float64, SplitFunction::typeof(MultivariateFunctions.add_split_with_max_function), removeSplitFunction::Bool)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:28
 [15] #56
    @ ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:36 [inlined]
 [16] optimize(f::MultivariateFunctions.var"#56#57"{DataFrame, Symbol, Vector{Piecewise_Function}, Int64, Symbol, typeof(MultivariateFunctions.add_split_with_max_function), Bool}, x_lower::Float64, x_upper::Float64, mo::Optim.Brent; rel_tol::Float64, abs_tol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, callback::Nothing, show_every::Int64, extended_trace::Bool)
    @ Optim ~/.julia/packages/Optim/TNmSw/src/univariate/solvers/brent.jl:69
 [17] #optimize#80
    @ ~/.julia/packages/Optim/TNmSw/src/univariate/optimize/interface.jl:21 [inlined]
 [18] optimise_split(dd::DataFrame, y::Symbol, array_of_funcs::Vector{Piecewise_Function}, ind::Int64, split_variable::Symbol, rel_tol::Float64, SplitFunction::Function, removeSplitFunction::Bool)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:36
 [19] create_mars_spline(dd::DataFrame, y::Symbol, x_variables::Set{Symbol}, MaxM::Int64; rel_tol::Float64)
    @ MultivariateFunctions ~/.julia/packages/MultivariateFunctions/a2IXB/src/6_HighDimensionalApproximation.jl:104
 [20] top-level scope
    @ REPL[18]:1

Notes

• versioninfo()

Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, westmere)

• MultivariateFunctions v0.1.7

Overhaul optimization to reduce computation in MARS

See start of page 6 of https://statistics.stanford.edu/sites/g/files/sbiybj6031/f/LCS%20110.pdf where this is said to reduce the amount of OLS computations.

Implement Fast MARS

Implement Fast MARS parent priority queue. See page 6 of:
https://statistics.stanford.edu/sites/g/files/sbiybj6031/f/LCS%20110.pdf