Use of Measurements.jl values works in some cases and fails in others.

[Boxylmer]

Hi!

When using Measurement types, Polynomials plays nicely with Measurements.jl until the polynomials get large enough / complicated enough in fitting that Vandermonde is used. Is there any way to skip this and use a nieve fitting method when I'm needing to use these types to propagate error?

Working minimal example: You will find that the following works

julia> Polynomials.fit([1, 2, 3] .± 0.3, [4, 5, 6.3].±0.2, 2)
> Polynomials.Polynomial(3.3 ± 1.6 + 0.55 ± 1.9*x + 0.15 ± 0.49*x^2)

and this will not

julia> Polynomials.fit([1, 2, 3] .± 0.3, [4, 5, 6.3].±0.2, 4)
ERROR: MethodError: no method matching Float64(::Measurement{Float64})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
 [1] convert(#unused#::Type{Float64}, x::Measurement{Float64})
   @ Base .\number.jl:7
 [2] setindex!(A::Vector{Float64}, x::Measurement{Float64}, i1::Int64)
   @ Base .\array.jl:966
 [3] vander(P::Type{Polynomials.Polynomial}, x::Vector{Measurement{Float64}}, degs::UnitRange{Int64})

[Boxylmer]

After additional crawling through the source code, I have realized this is not feasible, I apologize for asking a bit too early!

For those coming across this: My strategy for anyone wanting to combine Measurements.jl values when fitting in Polynomials.jl is to strip the uncertainties from the measurement types, weight them accordingly ($1/\sigma^2$), then use an appropriate error propagation method through the resulting fit polynomial. (i.e., jackknifing, bootstrapping, or acquiring the covariance matrix through the scaled inverse hessian of the maximum likelihood estimator.)

[jverzani]

Yes, that is a plan. We could do this in an extension with v1.9. If you were interested, a PR would be quite welcome

[Boxylmer]

@jverzani I actually am extremely confused now. Since my MWE actually runs correctly inside my package. The only difference I can find is that the package is using 3.2.10.

MWE in question:

using Polynomials
using Measurements

x = [0, 27.3, 66.2, 111.3, 134, 202.6, 256.8, 296.8, 358.2, 407.4, 453.4, 501.4] .* 0.00689476
y = [0 ± 0, 0.010053518 ± 0.000577937, 0.01908272 ± 0.001100238, 0.02783634 ± 0.001609498, 0.030801961 ± 0.00178267, 0.042389159 ± 0.002462366, 0.050273318 ± 0.002927624, 0.05573762 ± 0.003251388, 0.065036662 ± 0.003804798, 0.070552007 ± 0.004134467, 0.076878536 ± 0.004513929, 0.082775442 ± 0.004868871, ]
fit(x, y, 4)

It looks like something in 3.2.11 or 3.2.12 broke compatibility. My understanding of the issue is limited, but it looks like the strict typing on the matrix used in solving Vandermonde.

In the most recent version this function is where bad things are happening

Polynomials.jl/src/polynomials/standard-basis.jl

Lines 514 to 529 in fd31502

	function vander(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, degs) where {T <: Number}
	A = Matrix{T}(undef, length(x), length(degs))
	Aᵢ = one.(x)
	i′ = 1
	for i ∈ 0:maximum(degs)
	if i ∈ degs
	A[:, i′] = Aᵢ
	i′ += 1
	end
	for (i, xᵢ) ∈ enumerate(x)
	Aᵢ[i] *= xᵢ
	end
	end
	A
	end

Looks like A_i is being initialized as a matrix of floats and thus can't handle Measurements. Though the block in 3.2.10 is the exact same. I'm not sure what changed that's not allowing the input to vander to be of type measurement. Any ideas?

[Boxylmer]

I apologize, the versioning necessary is actually [email protected]

[Boxylmer]

I found it!

JuliaPhysics/Measurements.jl#134

using one() as of Measurements.jl 1.9 will no longer return a measurement type. This is why it starts to fail! To fix this, we only need to determine the type in another way. Is this possible?

[jverzani]

Yes, thanks for that detective work. I just put in a fix replacing one.(x) with ones(T, length(x)).