Pontius

test_pontius.jl

using LinearRegressionKit, StatsModels, DataFrames, CSV

y = [0.11019, 0.21956, 0.32949, 0.43899, 0.54803, 0.65694, 0.76562, 0.87487, 0.98292, 1.09146, 1.20001, 1.30822, 1.41599, 1.52399, 1.63194, 1.73947, 1.84646, 1.95392, 2.06128, 2.16844, 0.11052, 0.22018, 0.32939, 0.43886, 0.54798, 0.65739, 0.76596, 0.87474, 0.983, 1.0915, 1.20004, 1.30818, 1.41613, 1.52408, 1.63159, 1.73965, 1.84696, 1.95445, 2.06177, 2.16829]
x = [150000, 300000, 450000, 600000, 750000, 900000, 1050000, 1200000, 1350000, 1500000, 1650000, 1800000, 1950000, 2100000, 2250000, 2400000, 2550000, 2700000, 2850000, 3000000, 150000, 300000, 450000, 600000, 750000, 900000, 1050000, 1200000, 1350000, 1500000, 1650000, 1800000, 1950000, 2100000, 2250000, 2400000, 2550000, 2700000, 2850000, 3000000]

df = DataFrame(x= x, y= y)

f = @formula(y ~ x + x^2)
lrk= regress(f, df, req_stats=["default"])
lrk

delta(a,b) = abs(a-b) ≈ 0. ?  0. : abs(a-b)

println("Pontius - LRK")

println(string("delta B0 ", delta(0.673565789473684E-03, lrk.coefs[1])))
println(string("delta B1 ", delta(0.732059160401003E-06, lrk.coefs[2])))
println(string("delta B2 ", delta(-0.316081871345029E-14, lrk.coefs[3])))
println(string("delta std err B0 ", delta(0.107938612033077E-03, lrk.stderrors[1])))
println(string("delta std err B1 ", delta(0.157817399981659E-09, lrk.stderrors[2])))
println(string("delta std err B2 ", delta(0.486652849992036E-16, lrk.stderrors[3])))
println(string("delta resid std dev ", delta(0.205177424076185E-03, lrk.RMSE)))
println(string("delta R2 ", delta(0.999999900178537, lrk.R2)))

using GLM
println("Pontius - GLM")
lr = GLM.lm(f, df)

println(string("delta B0 ", delta(0.673565789473684E-03, coef(lr)[1])))
println(string("delta B1 ", delta(0.732059160401003E-06, coef(lr)[2])))
println(string("delta B2 ", delta(-0.316081871345029E-14, coef(lr)[3])))
println(string("delta std err B0 ", delta(0.107938612033077E-03, stderror(lr)[1])))
println(string("delta std err B1 ", delta(0.157817399981659E-09, stderror(lr)[2])))
println(string("delta std err B2 ", delta(0.486652849992036E-16, stderror(lr)[3])))
println(string("delta resid std dev ", delta(0.205177424076185E-03, deviance(lr)/dof_residual(lr))))
println(string("delta R2 ", delta(0.999999900178537, r2(lr))))