Note the pattern.
The Lagrange Polynomials can be efficiently evaluated using Neville's Algorithm, or Divided Differences. See Hoffman, page 201, 204.
The linear system is given as follows.
using Plots
using PyFormattedStrings
xdata = [ 0.34482759, 0.68965517, 1.03448276, 1.37931034,
1.72413793, 2.06896552, 2.4137931 , 2.75862069, 3.10344828,
3.44827586, 3.79310345, 4.13793103, 4.48275862, 4.82758621,
5.17241379, 5.51724138, 5.86206897, 6.20689655, 6.55172414,
6.89655172, 7.24137931, 7.5862069 , 7.93103448, 8.27586207,
8.62068966, 8.96551724, 9.31034483, 9.65517241, 10. ]
ydata = [ 4.55459724, 5.2125257 , 5.85608367, 6.28256403,
5.3510909 , 7.33644195, 8.74308925, 8.80680797, 7.91569673,
8.80173903, 8.09422165, 8.75590574, 10.68847556, 10.08827538,
9.67774597, 9.75957529, 10.68388337, 10.39037968, 10.83938351,
11.8776604 , 11.97929326, 13.7156162 , 12.63125149, 13.24891817,
12.23343525, 12.55352271, 13.23694 , 14.91190807, 14.93770334]
#-----------------------
A = [length(xdata) sum(log.(xdata));
sum(log.(xdata)) sum(log.(xdata).^2)]
b = [ sum(log.(ydata)),
sum(log.(xdata).*log.(ydata)) ]
p = A\b
α = exp(p[1])
β = p[2]
println(f"α = {α:.3f}, β = {β:.3f}")
α = 5.713, β = 0.366
xx = LinRange(0,10,1000)
yy = α*xx.^β
a = 5
b = 0.45
xb = LinRange(0,10,30)
yb = (5*xb.^0.45) # .+ 3*(rand(length(x)) .- 0.5)
scatter(xdata,ydata, label="data", lw=2)
plot!(xx,yy, label="model")
plot!(xb,yb, linestyle=:dash, color="gray", label="real", lw=2)
plot!(legend_foreground_color=nothing)
plot!(xlabel="x", ylabel="y")