As always, thanks for the suggestion. For some reason the legend includes labels from the data series of the first plot. Here is an example based on your code:

##
Summary

```
using Distributions, SmoothingSplines, Plots;
pyplot()
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
p2
temp = [-10 -20; -30 -40 ]
plot!(p2, temp, ylims = (-3,3), leg=true, label = ["a" "b"])
```

Is there any way to suppress the legend from the first data series?