using Distributions, Zygote, Plots, StatsPlots, Random
Censored Data
https://juliastats.org/Distributions.jl/stable/censored/
Random.seed!(123)
= collect(-1:0.01:3)
line
= 1.5
censoring = rand(Normal(1,0.5),500)
actual
= min.(censoring,actual)
observed
histogram(observed,
=:true,
normalize= "Observed, censored data",
label =0.5,
alpha=:topleft,
legend=:png
fmt
)plot!(line,
pdf.(Normal(1,0.5),line),
=3,
lw=:blue,
c= "True Distribution"
label
)vline!([censoring],
=3, c="red",
lw=:dot,
ls="Censoring point"
label )
= mean(observed)
mean_full = std(observed)
std_full
= mean(observed[observed .< censoring])
mean_red = std(observed[observed .< censoring])
std_red
= [0.,0.]
ps
for i in 1:50
= Zygote.gradient(p->-mean(logpdf.(censored.(Normal(p[1],exp(p[2])),-Inf,ones(length(observed)).*censoring),observed)),ps)[1]
grads
.-= 0.1 .* grads
ps
end
= ps[1]
mean_model = exp(ps[2])
std_model
histogram(observed,
=:true,
normalize= "Observed, censored data",
label =0.5,
alpha=:topleft,
legend=:png,
fmt=(900,600)
size
)
vline!([censoring],
=3, c="red",
lw=:dot,
ls="Censoring point"
label
)
plot!(line,
pdf.(Normal(1,0.5),line),
=3,
lw=:blue,
c= "True Distribution"
label
)
plot!(line,
pdf.(Normal(mean_model,std_model),line),
=3,
lw=:green,
c=:dash,
s= "With proper censoring model"
label
)
plot!(line,
pdf.(Normal(mean_full,std_full),line),
=1,
lw=:orange,
c= "Directly from data"
label
)
plot!(line,
pdf.(Normal(mean_red,std_red),line),
=1,
lw=:purple,
c= "Without censored observations"
label )