{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Illustration of robust polynomial fitting \n",
"2019-02-28 Jeff Fessler, University of Michigan "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra: Diagonal, svd\n",
"using Random: seed!\n",
"using Plots\n",
"include(ENV[\"HOME\"] * \"/l/g/teach/w19-598-opt/hw/auto/test/ncg_inv.jl\")\n",
"using Distributions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"s = (t) -> atan(4*(t-0.5)); # nonlinear function"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"seed!(1) # seed rng\n",
"M = 15\n",
"tm = sort(rand(M)) # M random sample locations\n",
"y = s.(tm) + 0.1 * randn(M); # noisy samples"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t0 = LinRange(0, 1, 101) # fine sampling for showing curve\n",
"scatter(tm, y, color=:blue,\n",
" label=\"y\", xlabel=\"t\", ylabel=\"y\", ylim=[-1.3, 1.3])\n",
"plot!(t0, s.(t0), color=:blue, label=\"s(t)\", legend=:topleft)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15×4 Array{Float64,2}:\n",
" 1.0 0.00790928 6.25568e-5 4.94779e-7\n",
" 1.0 0.210968 0.0445076 0.00938968\n",
" 1.0 0.236033 0.0557117 0.0131498 \n",
" 1.0 0.251662 0.0633339 0.0159387 \n",
" 1.0 0.28119 0.0790679 0.0222331 \n",
" 1.0 0.312707 0.0977856 0.0305783 \n",
" 1.0 0.346517 0.120074 0.0416077 \n",
" 1.0 0.424718 0.180385 0.0766128 \n",
" 1.0 0.437108 0.191063 0.0835153 \n",
" 1.0 0.488613 0.238742 0.116653 \n",
" 1.0 0.555751 0.308859 0.171649 \n",
" 1.0 0.773223 0.597874 0.46229 \n",
" 1.0 0.951916 0.906145 0.862574 \n",
" 1.0 0.986666 0.973511 0.96053 \n",
" 1.0 0.999905 0.999809 0.999714 "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"deg = 3 # polynomial degree\n",
"Afun = (tt) -> [t.^i for t in tt, i in 0:deg] # matrix of monomials\n",
"A = Afun(tm) # M × 4 matrix"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Float64,1}:\n",
" -1.1217772852297472\n",
" -0.7987637026605704\n",
" 9.106395885988642 \n",
" -6.099750221152859 "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xh = A \\ y # ordinary LS solution"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!(t0, Afun(t0)*xh, line=:red, label=\"ordinary LS cubic fit\")"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Float64,1}:\n",
" -1.3165917881067561\n",
" 3.2310676853223086\n",
" -3.94749649713934 \n",
" 3.1398714413842437"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"yc = copy(y)\n",
"yc[12] = -0.5\n",
"#seed!(8)\n",
"#yc = s.(tm) + 0.2 * rand(Laplace(0,1), M)\n",
"xc = A \\ yc # ordinary LS solution"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter(tm, yc, color=:blue,\n",
" label=\"y\", xlabel=\"t\", ylabel=\"y\", ylim=[-1.3, 1.3])\n",
"plot!(t0, s.(t0), color=:blue, label=\"s(t)\", legend=:topleft)\n",
"plot!(t0, Afun(t0)*xc, line=:red, label=\"ordinary LS cubic fit\")"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([-1.14924, -0.517269, 8.83346, -6.09636], Any[0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"delta = 0.01\n",
"dpot = (z,delta) -> z / (1 + abs(z)/delta) # Fair potential\n",
"xr,_ = ncg_inv([A], [v -> dpot.(v - yc,delta)], [1.], xc, niter=100)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!(t0, Afun(t0)*xr, line=:magenta, label=\"robust cubic fit\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"#savefig(\"tmp.pdf\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.1.0",
"language": "julia",
"name": "julia-1.1"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.1.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}