Hello,
I have two curves with the values:
col_y = [8.02511788606e7, 9.98481296252e7, 2.328332483329e8, 6.148181496526e8, 9.453171952624e8, 1.2331292887706e9, 1.4394296796998e9, 1.460066157955e9, 1.3570403824485e9, 1.4269447174002e9]
vin_y = [1.434628836798e8, 1.318375623746e8, 1.108524843855e8, 6.91176196241e7, 4.19137382458e7, 1.48676444645e7, 1.00659726053e7, 7.8734512643e6]
Is it possible to find the intersection x between these two curves? what is the approach in julia?
Thank you
Do you have a model for what kind of curves these are (they look sigmoid, do you know what process generated them)? And is the x-scale the same or do they for instance span between the same limits (which is not the same thing since the number of elements is different)? Linear x-scale?
The reason you need to know is the answer will vary wildly depending on which way you interpret the question. A vector of values isn’t a function, in fact there is an uncountable infinity of functions which fit exactly.
The way I would do this is
f_col, f_vin that fits each of the sample vectors (for instance using LsqFit.jl)g(x) = f_col(x) - f_vin(x)
g(x) (for instance using Roots.jl)You have points, not curves. The names suggest they are y-values, but it’s not clear what the corresponding x-values are for each.
Plotting these by assuming x=1:10 for the first and x=1:8 for the second gives:
So there is an intersection just after x=2, but this is based on a linear interpolation between points which Plots does automatically. As Gustaphe says you might have some other function in mind to interpolate between the points, which will result in a different point of intersection.
For this problem, Dierckx is just fine:
using Dierckx
n = length(vin_y)
spl = Spline1D(1:n, col_y[1:n] - vin_y)
x0 = roots(spl) # 2.4181
y0 = Spline1D(1:n, vin_y)(x0) # 1.2583e8
Thanks for the answers. I made a model, but it does not exactly fit the data:
Thus, I wanted to check the time difference of intersection between the actual data and the model…
I’ll check the suggestions.
Here’s my way. I have an experimental version of LsqFit installed, so there will be a slight difference in syntax, but I just realized I’m in a bit of a hurry so I won’t fix that.
using Plots, LsqFit, Roots
default(;fontfamily="Computer Modern", label=nothing, seriestype=:scatter) # Better plot defaults
x_1 = 1:10
x_2 = 1:8
y_1 = [8.02511788606e7, 9.98481296252e7, 2.328332483329e8, 6.148181496526e8, 9.453171952624e8, 1.2331292887706e9, 1.4394296796998e9, 1.460066157955e9, 1.3570403824485e9, 1.4269447174002e9]
y_2 = [1.434628836798e8, 1.318375623746e8, 1.108524843855e8, 6.91176196241e7, 4.19137382458e7, 1.48676444645e7, 1.00659726053e7, 7.8734512643e6]
model(x, p) = 1e9*p[1]/(1+exp(-(x-p[2])/p[3])) # logistic function
f_1 = curve_fit(model, x_1, y_1, [1.0, 1.0, 1.0])
f_2 = curve_fit(model, x_2, y_2, [1e-1, 1.0, -1.0])
x_long = range(0.5,10.5;length=200)
pl_1 = plot(x_1, y_1)
plot!(pl_1, x_2, y_2)
plot!(pl_1, x_long, model.(x_long,Ref(f_1.param)); st=:line)
plot!(pl_1, x_long, model.(x_long,Ref(f_2.param)); st=:line)
g(x) = model(x, f_1.param) - model(x, f_2.param)
pl_2 = plot(g, x_long, st=:line)
plot(pl_1, pl_2; layout=(2,1)
find_zero(g, 2.0) # returns ~ 2.27
If you want the intersection between the actual data series, the answer is there is none. None of your sample points measure the same amount of bacteria of the two species.
Thank you!
the Dierckx solution worked out of the box. By using the actual x values (sorry, I forgot to add them originally)
x = [0.0, 6.465599999999999, 12.544800000000002, 24.0, 29.107200000000002, 45.772800000000004, 53.392799999999994, 68.0616]
I got the exact intersection:
The more canonical solution by gustaphe is giving me troubles. The definition of the model is very similar to the model I have set:
model = true_logistic(x, p) = p[1] ./ (1 .+ exp.(-p[2].*(x.-p[3])))
that is a logistic regression but I set 1 instead of 1e9. When I define g and find_zero:
julia> g(x) = model(col_x[1:8], P_col) - model(col_x[1:8], P_vin)
g (generic function with 1 method)
julia> find_zero(g, col_x[2])
ERROR: Roots.ConvergenceFailed("Stopped")
Stacktrace:
[1] find_zero(fs::Function, x0::Float64, M::Secant, N::AlefeldPotraShi; p::Nothing, tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Roots ~/.julia/packages/Roots/H7pXT/src/order0.jl:57
[2] #find_zero#52
@ ~/.julia/packages/Roots/H7pXT/src/order0.jl:37 [inlined]
[3] find_zero
@ ~/.julia/packages/Roots/H7pXT/src/order0.jl:35 [inlined]
[4] #find_zero#16
@ ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704 [inlined]
[5] find_zero(f::Function, x0::Float64)
@ Roots ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704
[6] top-level scope
@ none:1
julia> find_zero(g, 2.0)
ERROR: Roots.ConvergenceFailed("Stopped")
Stacktrace:
[1] find_zero(fs::Function, x0::Float64, M::Secant, N::AlefeldPotraShi; p::Nothing, tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Roots ~/.julia/packages/Roots/H7pXT/src/order0.jl:57
[2] #find_zero#52
@ ~/.julia/packages/Roots/H7pXT/src/order0.jl:37 [inlined]
[3] find_zero
@ ~/.julia/packages/Roots/H7pXT/src/order0.jl:35 [inlined]
[4] #find_zero#16
@ ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704 [inlined]
[5] find_zero(f::Function, x0::Float64)
@ Roots ~/.julia/packages/Roots/H7pXT/src/find_zero.jl:704
[6] top-level scope
@ none:1
What am I getting wrong?
Likely you have a bad starting point. With only 1 intersection you should be able to use a bracketing method: find_zero(g, 1, 50)
The reason I multiplied by 1e9 is that LsqFit (and any fitting method) is better at finding similarly sized parameters, so I scaled up the model so we can expect p.≈1.0.
Your problem is in the definition of g. It needs to depend on x, not col_x.
Ah, so simple. It worked!
julia> g(x) = model(x, P_col) - model(x, P_vin)
g (generic function with 1 method)
julia> find_zero(g, col_x[2])
10.018996507513226
Thank you
Glad it worked out!
Can I ask what function you’re fitting to in the figure? Because the logistic function seems like a better fit, and from the little I know of microbiology it’s not an unreasonable model, so it seems a little strange to only use the better fit to verify the worse.
Hi, it is the true_logistic; I see that does not fit the data completely but it is the same model that I used for another group of bacteria (where it was more fitting). For consistency, I am sticking to this…
Ah okay, I missed the log axis. The fit is decent 
