Suppose I have a time series plot, such as:
using Plots
t = 0:.1:10
n = length(t)
m = sin.(t)
y = m + randn(n) * .3
plot(t, m)
scatter!(t, y)
How can I add a gradient to the alpha value of the ribbon based on the density at each value along the x-axis? An example of what I am trying to accomplish can be found here.
This is an attempt to reproduce your example using Plots.jl and to learn Julia:
using Distributions, SmoothingSplines, Plots; gr()
f(x) = sin(x) # baseline function for this example
# GENERATE INPUT DATA:
σy = 1.0 # stdev relative to curve f(x)
npoints = 100
nsamples = 50
X = vcat([rand(Uniform(0, 4π), npoints) for _ in 1:nsamples]...)
Y = @. f.(X) + rand(Normal(0, σy*abs(cos(X))+0.1)) # variable stdev
# FIT NORMAL DISTRIBUTIONS WITH PARAMETERS: μ(x) and σ(x)
nx = 100; # number of divisions of x-axis to estimate distribution
x = LinRange(0, 4π, nx)
xᵢ = [(x[i], x[i+1]) for i in 1:length(x)-1]
x₀ = mean.(xᵢ)
μ = similar(x₀); σ = similar(x₀)
for (i,xᵢ) in pairs(xᵢ)
ix = xᵢ[1] .<= X .<= xᵢ[2]
h = fit(Normal, Y[ix])
μ[i] = h.μ; σ[i] = h.σ
end
# FIT SMOOTHING SPLINES TO ABOVE
using SmoothingSplines
splμ = fit(SmoothingSpline,x₀, μ, 0.05) # λ=0.05
splσ = fit(SmoothingSpline,x₀, σ, 0.02) # λ=0.02
μp = SmoothingSplines.predict(splμ, x₀)
σp = SmoothingSplines.predict(splσ, x₀)
# PLOT DATA:
p1 = scatter(X,Y, ms=1, msw=-0, mc=:blue, label=false)
plot!(f, lc=:cyan, lw=1.5, label="sin(x) baseline")
plot!(x₀, μ, lc=:red, lw=1, ls=:dash, label="mean μ(x)")
plot!(x₀, σ, lc=:lime, lw=1, ls=:dash, label="stdev σ(x)")
plot!(x₀, μp, lc=:red, lw=1, label="spline μp(x)")
plot!(x₀, σp, lc=:lime, lw=1, label="spline σp(x)")
# PLOT RIBBONS WITH ALPHA TRANSPARENCY SCALED:
N = 300 # number of ribbon slices
qq = @. quantile(Normal(μp, σp), LinRange(0.01, 0.99,N))
α = [LinRange(0.01,0.4,N÷2); LinRange(0.4,0.01,N÷2)]
p2 = plot(legend=false)
for i in 2:N-1
yᵢ = getindex.(qq, i)
dy = yᵢ - getindex.(qq, i-1)
plot!(x₀, yᵢ - dy/2, lw=0, color=:blue, fillalpha=α[i], ribbon=dy)
end
scatter!(X,Y, ms=1, msw=-0, mc=:darkblue, label=false)
plot(p1, p2, size=(1200,600), dpi=300, widen=false)
Looks good. Thanks!