From 65f76520554d892c87495a229c7a1eb97733ee3f Mon Sep 17 00:00:00 2001 From: math4mad Date: Fri, 20 Oct 2023 10:07:47 +0800 Subject: [PATCH] add note add understanding bayes note --- docs/index.html | 6 + .../1-classfication-comparison.html | 6 + docs/materials.html | 6 + docs/math/1-double-integral.html | 6 + docs/math/2-function.html | 6 + docs/math/3-oval-spere.html | 6 + docs/math/4-one-variable-integral.html | 6 + docs/math/5-himmelblau-function.html | 6 + docs/schedule.html | 6 + docs/search.json | 128 ++-- docs/statistics/1-normaldist.html | 6 + docs/statistics/10-dirichlet-dist.html | 6 + docs/statistics/12-likelihood-coin-toss.html | 65 +- docs/statistics/13-understand-bayes.html | 665 ++++++++++++++++++ .../figure-html/cell-4-output-1.png | Bin 0 -> 18956 bytes .../figure-html/cell-5-output-1.png | Bin 0 -> 20737 bytes .../figure-html/cell-6-output-1.png | Bin 0 -> 18476 bytes docs/statistics/2-2d-descison-boundary.html | 6 + docs/statistics/3-bio-norm-dist.html | 6 + docs/statistics/4-pair-plots.html | 6 + docs/statistics/5-2d-decision-boundary.html | 6 + .../6-mutli-variable-distributions.html | 6 + docs/statistics/7-mvnormal-sampling.html | 6 + docs/statistics/8-dotplot.html | 6 + docs/statistics/9-2d-decision-boundary.html | 6 + statistics/12-likelihood-coin-toss.qmd | 12 + statistics/13-understand-bayes.qmd | 64 ++ 27 files changed, 993 insertions(+), 55 deletions(-) create mode 100644 docs/statistics/13-understand-bayes.html create mode 100644 docs/statistics/13-understand-bayes_files/figure-html/cell-4-output-1.png create mode 100644 docs/statistics/13-understand-bayes_files/figure-html/cell-5-output-1.png create mode 100644 docs/statistics/13-understand-bayes_files/figure-html/cell-6-output-1.png create mode 100644 statistics/13-understand-bayes.qmd diff --git a/docs/index.html b/docs/index.html index ec11f1c..ecb903d 100644 --- a/docs/index.html +++ b/docs/index.html @@ -172,6 +172,12 @@ 12-likelihood of coin toss + +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/machinelearning/1-classfication-comparison.html b/docs/machinelearning/1-classfication-comparison.html index 39a0b6a..a2b9585 100644 --- a/docs/machinelearning/1-classfication-comparison.html +++ b/docs/machinelearning/1-classfication-comparison.html @@ -206,6 +206,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/materials.html b/docs/materials.html index d24abe8..575ac94 100644 --- a/docs/materials.html +++ b/docs/materials.html @@ -172,6 +172,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/math/1-double-integral.html b/docs/math/1-double-integral.html index 86c1fbd..c4b4b0a 100644 --- a/docs/math/1-double-integral.html +++ b/docs/math/1-double-integral.html @@ -206,6 +206,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/math/2-function.html b/docs/math/2-function.html index af0a082..9149bdc 100644 --- a/docs/math/2-function.html +++ b/docs/math/2-function.html @@ -206,6 +206,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/math/3-oval-spere.html b/docs/math/3-oval-spere.html index 5152237..7076f47 100644 --- a/docs/math/3-oval-spere.html +++ b/docs/math/3-oval-spere.html @@ -206,6 +206,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/math/4-one-variable-integral.html b/docs/math/4-one-variable-integral.html index 97be8d7..626ca41 100644 --- a/docs/math/4-one-variable-integral.html +++ b/docs/math/4-one-variable-integral.html @@ -206,6 +206,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/math/5-himmelblau-function.html b/docs/math/5-himmelblau-function.html index dca56eb..c3f8bf2 100644 --- a/docs/math/5-himmelblau-function.html +++ b/docs/math/5-himmelblau-function.html @@ -206,6 +206,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/schedule.html b/docs/schedule.html index 6b79de4..f68a009 100644 --- a/docs/schedule.html +++ b/docs/schedule.html @@ -172,6 +172,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/search.json b/docs/search.json index 6b8c8ef..7112cd9 100644 --- a/docs/search.json +++ b/docs/search.json @@ -98,46 +98,32 @@ "text": "Code\nusing GLMakie\nusing Distributions\nusing Latexify\n\nGLMakie.activate!()\n\nn=100\nxs=ys=range(-4,4,n)\nμ=[0 , 0]\nΣ=[1 0 ; 0 1]\nd=MvNormal(μ,Σ)\nd2=Normal()\nzs=[pdf(d,[x,y]) for x in xs , y in ys]\n\n# 投影到 x-z 平面数据\nxz=[Point3f(4,y,pdf(d2,y)/2) for y in ys]\n\n# 投影到 y-z 平面数据\nyz=[Point3f(x,4,pdf(d2,x)/2) for x in ys]\n\n\nstr=latexify([0,0])\nfig = Figure(size=(600,600))\nax = Axis3(fig[1, 1],title=L\"\\mu \")\nsurface!(xs,ys,zs;space=:data)\nlines!(ax, xz;linewidth=3)\nlines!(ax, yz;linewidth=3)\n\nfig\n#save(\"./imgs/mutli-normal-dist.png\",fig)" }, { - "objectID": "statistics/5-2d-decision-boundary.html", - "href": "statistics/5-2d-decision-boundary.html", - "title": "5-2d-decision-boundary", + "objectID": "statistics/3-bio-norm-dist.html", + "href": "statistics/3-bio-norm-dist.html", + "title": "3-biovariate-normal-dist", "section": "", - "text": "Code\nimport MLJ:fit!,predict,predict_mode,predict_mean,machine\nusing MLJ,GLMakie,Random,DataFrames\n\niris = load_iris(); \niris =DataFrame(iris);\nnums=100\n\niris[!, :target] = [r.target == \"virginica\" ? 1.0 : 0.0 for r in eachrow(iris)]\niris=coerce(iris, :target=> Multiclass )\ngdf=groupby(iris, :target)\nX,y=iris[:,3:4],iris[:,:target]\n\ncats=levels(y)\n\nfunction boundary_data(df,;n=nums)\n n1=n2=n\n xlow,xhigh=extrema(df[:,1])\n ylow,yhigh=extrema(df[:,2])\n tx = LinRange(xlow,xhigh,n1)\n ty = LinRange(ylow,yhigh,n2)\n x_test = mapreduce(collect, hcat, Iterators.product(tx, ty));\n x_test=MLJ.table(x_test')\n return tx,ty,x_test\nend\ntx,ty,x_test=boundary_data(X)\n\n\nusing MLJ\nLogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels\n#X, y = make_blobs(centers = 2)\nmach = fit!(machine(LogisticClassifier(), X, y))\npredict(mach, X)\nfitted_params(mach)\nprobs=predict(mach, x_test)|>Array #返回分类概率值\nprobs_res=broadcast(pdf, probs, 1.0).|>(d->round(d,digits=2))|>d->reshape(d,nums,nums) #返回概率为1.0(\"virginica\")的概率值\n\n\n[ Info: For silent loading, specify `verbosity=0`. \n[ Info: Training machine(LogisticClassifier(lambda = 2.220446049250313e-16, …), …).\n┌ Info: Solver: MLJLinearModels.LBFGS{Optim.Options{Float64, Nothing}, NamedTuple{(), Tuple{}}}\n│ optim_options: Optim.Options{Float64, Nothing}\n└ lbfgs_options: NamedTuple{(), Tuple{}} NamedTuple()\n\n\nimport MLJLinearModels ✔\n\n\n100×100 Matrix{Float64}:\n 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n ⋮ ⋮ ⋱ ⋮ \n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 … 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.01 0.01 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.01 0.01 0.01 … 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.01 0.01 0.01 0.01 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.01 0.01 0.01 0.01 0.01 0.02 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.01 0.01 0.01 0.02 0.02 0.03 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.01 0.01 0.02 0.02 0.03 0.04 1.0 1.0 1.0 1.0 1.0 1.0 1.0" - }, - { - "objectID": "statistics/5-2d-decision-boundary.html#plot-results", - "href": "statistics/5-2d-decision-boundary.html#plot-results", - "title": "5-2d-decision-boundary", - "section": "plot results", - "text": "plot results\n\n\nCode\n fig=Figure(resolution=(1600,800))\n\n function plot_res_contour()\n\n ax = Axis(fig[1, 1], xlabel=\"petal-length\", ylabel=\"petal-width\", title=\"2d2class-contour\")\n contour!(ax, tx, ty, probs_res; labels=true)\n colors = [:red, :blue]\n for i in 1:2\n scatter!(ax, gdf[i][:, 3], gdf[i][:, 4], color=(colors[i], 0.8), marker=:circle, markersize=10, strokewidth=1, strokecolor=:black, label=gdf[i][1, 5] == 1 ? \"virginica\" : \"non-virginica\")\n end\n axislegend(ax, position=:lt)\n #save(\"./imgs/iris-logreg-2d-2class-contourf.png\",fig)\n fig\n end\n\n function plot_res_contourf()\n\n ax = Axis(fig[1, 2], xlabel=\"petal-length\", ylabel=\"petal-width\", title=\"2d2class-contourf\")\n contourf!(ax, tx, ty, probs_res; levels=6, colormap=(:heat, 0.5))\n #contourf!(ax,tx,ty,yhat;levels=length(cats),colormap=(:heat,0.5))\n colors = [:red, :blue]\n for i in 1:2\n scatter!(ax, gdf[i][:, 3], gdf[i][:, 4], color=(colors[i], 0.8), marker=:circle, markersize=10, strokewidth=1, strokecolor=:black, label=gdf[i][1, 5] == 1.0 ? \"virginica\" : \"non-virginica\")\n end\n axislegend(ax, position=:lt)\n #save(\"./imgs/iris-logreg-2d-2class-contourf.png\",fig)\n fig\n end\n\n plot_res_contourf()\n plot_res_contour()" + "text": "Code\n using GLMakie, Distributions\n using FileIO" }, { - "objectID": "statistics/9-2d-decision-boundary.html", - "href": "statistics/9-2d-decision-boundary.html", - "title": "9-2d-descision-boundary", + "objectID": "statistics/3-bio-norm-dist.html#load-package", + "href": "statistics/3-bio-norm-dist.html#load-package", + "title": "3-biovariate-normal-dist", "section": "", - "text": "简介\n\n\n\n\nref: probml page84 figure 2.13\ndataset:iris\nplots:使用 GLMakie:contourf 方法" - }, - { - "objectID": "statistics/9-2d-decision-boundary.html#load-package", - "href": "statistics/9-2d-decision-boundary.html#load-package", - "title": "9-2d-descision-boundary", - "section": "1. load package", - "text": "1. load package\n\n\nCode\n import MLJ:fit!,fitted_params\n using GLMakie,MLJ,CSV,DataFrames" + "text": "Code\n using GLMakie, Distributions\n using FileIO" }, { - "objectID": "statistics/9-2d-decision-boundary.html#process-data", - "href": "statistics/9-2d-decision-boundary.html#process-data", - "title": "9-2d-descision-boundary", - "section": "2 process data", - "text": "2 process data\n\n2.1 import iris datset\n\n\nCode\niris = load_iris();\n\n#selectrows(iris, 1:3) |> pretty\n\niris = DataFrames.DataFrame(iris);\nfirst(iris,5)|>display\ny, X = unpack(iris, ==(:target); rng=123);\n\nX=select!(X,3:4)\n\nbyCat = iris.target\ncateg = unique(byCat)\ncolors = [:orange,:lightgreen,:purple];\n\n\n5×5 DataFrame\n\n\n\nRow\nsepal_length\nsepal_width\npetal_length\npetal_width\ntarget\n\n\n\nFloat64\nFloat64\nFloat64\nFloat64\nCat…\n\n\n\n\n1\n5.1\n3.5\n1.4\n0.2\nsetosa\n\n\n2\n4.9\n3.0\n1.4\n0.2\nsetosa\n\n\n3\n4.7\n3.2\n1.3\n0.2\nsetosa\n\n\n4\n4.6\n3.1\n1.5\n0.2\nsetosa\n\n\n5\n5.0\n3.6\n1.4\n0.2\nsetosa\n\n\n\n\n\n\n\n\n2.2 make desc boundary data\n\n生成决策边界实际是利用训练模型对区间内的每个点都做出预测,利用两个属性的最大值和最小值 生成 grid 数据,这是 test数据\n\n\n\nCode\n# grid data\n n1 = n2 = 20\n tx = LinRange(0, 8, 40)\n ty = LinRange(-1, 4, 40)\n X_test = mapreduce(collect, hcat, Iterators.product(tx, ty))\n X_test = MLJ.table(X_test')\n\n\nTables.MatrixTable{LinearAlgebra.Adjoint{Float64, Matrix{Float64}}} with 1600 rows, 2 columns, and schema:\n :x1 Float64\n :x2 Float64" + "objectID": "statistics/3-bio-norm-dist.html#plot", + "href": "statistics/3-bio-norm-dist.html#plot", + "title": "3-biovariate-normal-dist", + "section": "2. plot", + "text": "2. plot\n\n\nCode\n μ = [0, 0]\n Σ₁ = [2 1.8; 1.8 2]\n Σ₂ = [1 0; 0 3]\n Σ₃ = [1 0; 0 1]\n simga_dict = Dict(\"full\" => Σ₁, \"diagonal\" => Σ₂, \"spherical\" => Σ₃)\n fig = Figure(resolution=(1200, 800))\n\n function load_img(src)\n load(\"./imgs/$(src).png\")\n end\n\n function binormal_plot(; μ=[0, 0], Σmatrix::Dict=simga_dict, up::Int=4, n::Int=100)\n xs = ys = range(-up, up, n)\n simga_keys = keys(Σmatrix)\n for (idx, val) in enumerate(simga_keys)\n local cormatrix = Σmatrix[val]\n local img=load_img(val)\n local d = MvNormal(μ, cormatrix)\n local zs = [pdf(d, [x, y]) for x in xs, y in ys]\n local ax1 = Axis3(fig[1, idx], title=\"$(val)\")\n local ax2= Axis(fig[2, idx],aspect = DataAspect(),height=60)\n local ax3 = Axis(fig[3, idx];)\n hidedecorations!(ax2)\n hidespines!(ax2)\n surface!(ax1, xs, ys, zs)\n image!(ax2, rotr90(img))\n contour!(ax3, xs, ys, zs, levels=10)\n end\n end\n\n binormal_plot()\n fig" }, { - "objectID": "statistics/9-2d-decision-boundary.html#logisitcs-model", - "href": "statistics/9-2d-decision-boundary.html#logisitcs-model", - "title": "9-2d-descision-boundary", - "section": "3. Logisitcs model", - "text": "3. Logisitcs model\n\n3.1 training model\n\n\nCode\n LogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels\n #X, y = make_blobs(centers = 2)\n mach = fit!(machine(LogisticClassifier(), X, y))\n predict(mach, X)\n fitted_params(mach)\n probs=predict(mach, X_test)|>Array #返回分类概率值\n probres=[broadcast(pdf, probs,cat) for cat in categ]\n\n\nimport MLJLinearModels ✔\n\n\n[ Info: For silent loading, specify `verbosity=0`. \n[ Info: Training machine(LogisticClassifier(lambda = 2.220446049250313e-16, …), …).\n┌ Info: Solver: MLJLinearModels.LBFGS{Optim.Options{Float64, Nothing}, NamedTuple{(), Tuple{}}}\n│ optim_options: Optim.Options{Float64, Nothing}\n└ lbfgs_options: NamedTuple{(), Tuple{}} NamedTuple()\n\n\n3-element Vector{Vector{Float64}}:\n [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 9.671834765600391e-132, 7.370874765820906e-135, 5.617320408186208e-138, 4.2809421636823794e-141, 3.2624925190464766e-144, 2.4863352574885955e-147, 1.8948282567827453e-150, 1.4440426374071983e-153, 1.1005003388488668e-156, 8.38687836794779e-160]\n [2.6707314921073497e-56, 1.0760270944086094e-53, 4.335270360659973e-51, 1.7466631832673993e-48, 7.037236485793986e-46, 2.835274587075858e-43, 1.142320852844153e-40, 4.6023652763324474e-38, 1.8542746623291115e-35, 7.470798854313656e-33 … 1.3627080472304016e-14, 4.184137656019423e-15, 1.2847218419308875e-15, 3.9446842977544977e-16, 1.2111986969540427e-16, 3.718934578207533e-17, 1.1418831965200673e-17, 3.5061042539858484e-18, 1.0765345419987361e-18, 3.305453963038536e-19]\n [1.6814079540479983e-80, 2.2062917118777524e-77, 2.8950280068446543e-74, 3.798766552624987e-71, 4.9846244275441556e-68, 6.540670593853871e-65, 8.582466430351137e-62, 1.1261648017760677e-58, 1.4777187549191659e-55, 1.9390170205959138e-52 … 0.9999999999999865, 0.9999999999999958, 0.9999999999999987, 0.9999999999999996, 0.9999999999999998, 1.0, 1.0, 1.0, 1.0, 1.0]\n\n\n\n\nCode\nfunction plot_res_contour()\n fig=Figure()\n ax = Axis(fig[1, 1], xlabel=\"petal-length\", ylabel=\"petal-width\", title=\"2d3class-contour\")\n \n \n for (idx,cat) in enumerate(categ)\n indc = findall(x -> x == cat, byCat)\n scatter!(ax,iris[:,3][indc],iris[:,4][indc];color=(colors[idx], 0.8), markersize=probres[idx].*10, strokewidth=1, strokecolor=:black, label=\"$cat\")\n \n end\n fig\nend\nplot_res_contour()" + "objectID": "statistics/4-pair-plots.html", + "href": "statistics/4-pair-plots.html", + "title": "4-pair-plots", + "section": "", + "text": "\"\"\"\n probml page 35 fig1.3\n findall(x -> x == c, byCat) 参考 beautiful makie 代码\n\"\"\"\n\n\nCode\nusing MLJ,DataFrames,GLMakie\n fontsize_theme = Theme(fontsize = 10)\n set_theme!(fontsize_theme)\n\n iris = load_iris(); \n iris = DataFrames.DataFrame(iris);\n #y, X = unpack(iris, ==(:target); rng=123);\n byCat = iris.target\n categ = unique(byCat)\n label=[\"sepal_length\",\"sepal_width\",\"petal_length\",\"petal_width\"]\n axs=[]\n colors1 = [:orange,:lightgreen,:purple]\n \n\n fig = Figure(resolution=(1400,1400))\n \n function plot_diag(i,j)\n \n ax = Axis(fig[i, i])\n push!(axs,ax)\n for (j, c) in enumerate(categ)\n indc = findall(x -> x == c, byCat)\n density!(ax, iris[:,i][indc]; color = (colors1[j], 0.5), label = \"$(c)\",\n strokewidth = 1.25, strokecolor = colors1[j])\n end\n end\n\n\n\"\"\"\n plot_cor(i,j)\n 生成非对角列的散点图\nTBW\n\"\"\"\nfunction plot_cor(i,j)\n ax = Axis(fig[i, j])\n #push!(axs,ax)\n for (k, c) in enumerate(categ)\n indc = findall(x -> x == c, byCat)\n #@show indc\n scatter!(ax, iris[:,i][indc],iris[:,j][indc];color=colors1[k])\n end\n end\n\n function plot_pair()\n [( i==j ? plot_diag(i,j) : plot_cor(i,j)) for i in 1:4,j in 1:4]\n end\n\n function add_xy_label()\n for i in 1:4\n axx=Axis(fig[4, i], xlabel =label[i],)\n axy=Axis(fig[i, 1], ylabel =label[i],)\n end\n end\n \n function add_legend()\n Legend(fig[2:3, 5], axs[1],\"Label\";width=100,height=200)\n end\n \n \n \n function main()\n plot_pair()\n add_xy_label()\n add_legend()\n fig\n #save(\"./imgs/iris-corner-plot2.png\",fig)\n end\n #save(\"./imgs/iris-corner-plot.png\",fig)\n\n\nmain (generic function with 1 method)\n\n\n\nPlotData\n\n\n\n\nCode\nwith_theme(fontsize_theme, fontsize = 25) do\n main()\nend\n\n\n\n\n\n\n\n\n\nCode\n first(iris,10)\n\n\n10×5 DataFrame\n\n\n\nRow\nsepal_length\nsepal_width\npetal_length\npetal_width\ntarget\n\n\n\nFloat64\nFloat64\nFloat64\nFloat64\nCat…\n\n\n\n\n1\n5.1\n3.5\n1.4\n0.2\nsetosa\n\n\n2\n4.9\n3.0\n1.4\n0.2\nsetosa\n\n\n3\n4.7\n3.2\n1.3\n0.2\nsetosa\n\n\n4\n4.6\n3.1\n1.5\n0.2\nsetosa\n\n\n5\n5.0\n3.6\n1.4\n0.2\nsetosa\n\n\n6\n5.4\n3.9\n1.7\n0.4\nsetosa\n\n\n7\n4.6\n3.4\n1.4\n0.3\nsetosa\n\n\n8\n5.0\n3.4\n1.5\n0.2\nsetosa\n\n\n9\n4.4\n2.9\n1.4\n0.2\nsetosa\n\n\n10\n4.9\n3.1\n1.5\n0.1\nsetosa" }, { "objectID": "statistics/8-dotplot.html", @@ -168,32 +154,74 @@ "text": "3. plot\n\n\nCode\n @rput df\n R\"\"\"\n library(ggplot2)\n ggplot(df, aes(x = x)) + \n geom_dotplot(dotsize = .75, stackratio = 1.2, \n fill = \"steelblue\") + \n scale_y_continuous(NULL, breaks = NULL) + \n labs(title = \"Stacked Dot Plot\", x = \"Data Values\") +\n theme_minimal()\n \"\"\"\n\n\n┌ Warning: RCall.jl: Bin width defaults to 1/30 of the range of the data. Pick better value with\n│ `binwidth`.\n└ @ RCall ~/.julia/packages/RCall/gOwEW/src/io.jl:172\n\n\n\n\n\nRObject{VecSxp}" }, { - "objectID": "statistics/4-pair-plots.html", - "href": "statistics/4-pair-plots.html", - "title": "4-pair-plots", + "objectID": "statistics/9-2d-decision-boundary.html", + "href": "statistics/9-2d-decision-boundary.html", + "title": "9-2d-descision-boundary", "section": "", - "text": "\"\"\"\n probml page 35 fig1.3\n findall(x -> x == c, byCat) 参考 beautiful makie 代码\n\"\"\"\n\n\nCode\nusing MLJ,DataFrames,GLMakie\n fontsize_theme = Theme(fontsize = 10)\n set_theme!(fontsize_theme)\n\n iris = load_iris(); \n iris = DataFrames.DataFrame(iris);\n #y, X = unpack(iris, ==(:target); rng=123);\n byCat = iris.target\n categ = unique(byCat)\n label=[\"sepal_length\",\"sepal_width\",\"petal_length\",\"petal_width\"]\n axs=[]\n colors1 = [:orange,:lightgreen,:purple]\n \n\n fig = Figure(resolution=(1400,1400))\n \n function plot_diag(i,j)\n \n ax = Axis(fig[i, i])\n push!(axs,ax)\n for (j, c) in enumerate(categ)\n indc = findall(x -> x == c, byCat)\n density!(ax, iris[:,i][indc]; color = (colors1[j], 0.5), label = \"$(c)\",\n strokewidth = 1.25, strokecolor = colors1[j])\n end\n end\n\n\n\"\"\"\n plot_cor(i,j)\n 生成非对角列的散点图\nTBW\n\"\"\"\nfunction plot_cor(i,j)\n ax = Axis(fig[i, j])\n #push!(axs,ax)\n for (k, c) in enumerate(categ)\n indc = findall(x -> x == c, byCat)\n #@show indc\n scatter!(ax, iris[:,i][indc],iris[:,j][indc];color=colors1[k])\n end\n end\n\n function plot_pair()\n [( i==j ? plot_diag(i,j) : plot_cor(i,j)) for i in 1:4,j in 1:4]\n end\n\n function add_xy_label()\n for i in 1:4\n axx=Axis(fig[4, i], xlabel =label[i],)\n axy=Axis(fig[i, 1], ylabel =label[i],)\n end\n end\n \n function add_legend()\n Legend(fig[2:3, 5], axs[1],\"Label\";width=100,height=200)\n end\n \n \n \n function main()\n plot_pair()\n add_xy_label()\n add_legend()\n fig\n #save(\"./imgs/iris-corner-plot2.png\",fig)\n end\n #save(\"./imgs/iris-corner-plot.png\",fig)\n\n\nmain (generic function with 1 method)\n\n\n\nPlotData\n\n\n\n\nCode\nwith_theme(fontsize_theme, fontsize = 25) do\n main()\nend\n\n\n\n\n\n\n\n\n\nCode\n first(iris,10)\n\n\n10×5 DataFrame\n\n\n\nRow\nsepal_length\nsepal_width\npetal_length\npetal_width\ntarget\n\n\n\nFloat64\nFloat64\nFloat64\nFloat64\nCat…\n\n\n\n\n1\n5.1\n3.5\n1.4\n0.2\nsetosa\n\n\n2\n4.9\n3.0\n1.4\n0.2\nsetosa\n\n\n3\n4.7\n3.2\n1.3\n0.2\nsetosa\n\n\n4\n4.6\n3.1\n1.5\n0.2\nsetosa\n\n\n5\n5.0\n3.6\n1.4\n0.2\nsetosa\n\n\n6\n5.4\n3.9\n1.7\n0.4\nsetosa\n\n\n7\n4.6\n3.4\n1.4\n0.3\nsetosa\n\n\n8\n5.0\n3.4\n1.5\n0.2\nsetosa\n\n\n9\n4.4\n2.9\n1.4\n0.2\nsetosa\n\n\n10\n4.9\n3.1\n1.5\n0.1\nsetosa" + "text": "简介\n\n\n\n\nref: probml page84 figure 2.13\ndataset:iris\nplots:使用 GLMakie:contourf 方法" }, { - "objectID": "statistics/3-bio-norm-dist.html", - "href": "statistics/3-bio-norm-dist.html", - "title": "3-biovariate-normal-dist", + "objectID": "statistics/9-2d-decision-boundary.html#load-package", + "href": "statistics/9-2d-decision-boundary.html#load-package", + "title": "9-2d-descision-boundary", + "section": "1. load package", + "text": "1. load package\n\n\nCode\n import MLJ:fit!,fitted_params\n using GLMakie,MLJ,CSV,DataFrames" + }, + { + "objectID": "statistics/9-2d-decision-boundary.html#process-data", + "href": "statistics/9-2d-decision-boundary.html#process-data", + "title": "9-2d-descision-boundary", + "section": "2 process data", + "text": "2 process data\n\n2.1 import iris datset\n\n\nCode\niris = load_iris();\n\n#selectrows(iris, 1:3) |> pretty\n\niris = DataFrames.DataFrame(iris);\nfirst(iris,5)|>display\ny, X = unpack(iris, ==(:target); rng=123);\n\nX=select!(X,3:4)\n\nbyCat = iris.target\ncateg = unique(byCat)\ncolors = [:orange,:lightgreen,:purple];\n\n\n5×5 DataFrame\n\n\n\nRow\nsepal_length\nsepal_width\npetal_length\npetal_width\ntarget\n\n\n\nFloat64\nFloat64\nFloat64\nFloat64\nCat…\n\n\n\n\n1\n5.1\n3.5\n1.4\n0.2\nsetosa\n\n\n2\n4.9\n3.0\n1.4\n0.2\nsetosa\n\n\n3\n4.7\n3.2\n1.3\n0.2\nsetosa\n\n\n4\n4.6\n3.1\n1.5\n0.2\nsetosa\n\n\n5\n5.0\n3.6\n1.4\n0.2\nsetosa\n\n\n\n\n\n\n\n\n2.2 make desc boundary data\n\n生成决策边界实际是利用训练模型对区间内的每个点都做出预测,利用两个属性的最大值和最小值 生成 grid 数据,这是 test数据\n\n\n\nCode\n# grid data\n n1 = n2 = 20\n tx = LinRange(0, 8, 40)\n ty = LinRange(-1, 4, 40)\n X_test = mapreduce(collect, hcat, Iterators.product(tx, ty))\n X_test = MLJ.table(X_test')\n\n\nTables.MatrixTable{LinearAlgebra.Adjoint{Float64, Matrix{Float64}}} with 1600 rows, 2 columns, and schema:\n :x1 Float64\n :x2 Float64" + }, + { + "objectID": "statistics/9-2d-decision-boundary.html#logisitcs-model", + "href": "statistics/9-2d-decision-boundary.html#logisitcs-model", + "title": "9-2d-descision-boundary", + "section": "3. Logisitcs model", + "text": "3. Logisitcs model\n\n3.1 training model\n\n\nCode\n LogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels\n #X, y = make_blobs(centers = 2)\n mach = fit!(machine(LogisticClassifier(), X, y))\n predict(mach, X)\n fitted_params(mach)\n probs=predict(mach, X_test)|>Array #返回分类概率值\n probres=[broadcast(pdf, probs,cat) for cat in categ]\n\n\nimport MLJLinearModels ✔\n\n\n[ Info: For silent loading, specify `verbosity=0`. \n[ Info: Training machine(LogisticClassifier(lambda = 2.220446049250313e-16, …), …).\n┌ Info: Solver: MLJLinearModels.LBFGS{Optim.Options{Float64, Nothing}, NamedTuple{(), Tuple{}}}\n│ optim_options: Optim.Options{Float64, Nothing}\n└ lbfgs_options: NamedTuple{(), Tuple{}} NamedTuple()\n\n\n3-element Vector{Vector{Float64}}:\n [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 9.671834765600391e-132, 7.370874765820906e-135, 5.617320408186208e-138, 4.2809421636823794e-141, 3.2624925190464766e-144, 2.4863352574885955e-147, 1.8948282567827453e-150, 1.4440426374071983e-153, 1.1005003388488668e-156, 8.38687836794779e-160]\n [2.6707314921073497e-56, 1.0760270944086094e-53, 4.335270360659973e-51, 1.7466631832673993e-48, 7.037236485793986e-46, 2.835274587075858e-43, 1.142320852844153e-40, 4.6023652763324474e-38, 1.8542746623291115e-35, 7.470798854313656e-33 … 1.3627080472304016e-14, 4.184137656019423e-15, 1.2847218419308875e-15, 3.9446842977544977e-16, 1.2111986969540427e-16, 3.718934578207533e-17, 1.1418831965200673e-17, 3.5061042539858484e-18, 1.0765345419987361e-18, 3.305453963038536e-19]\n [1.6814079540479983e-80, 2.2062917118777524e-77, 2.8950280068446543e-74, 3.798766552624987e-71, 4.9846244275441556e-68, 6.540670593853871e-65, 8.582466430351137e-62, 1.1261648017760677e-58, 1.4777187549191659e-55, 1.9390170205959138e-52 … 0.9999999999999865, 0.9999999999999958, 0.9999999999999987, 0.9999999999999996, 0.9999999999999998, 1.0, 1.0, 1.0, 1.0, 1.0]\n\n\n\n\nCode\nfunction plot_res_contour()\n fig=Figure()\n ax = Axis(fig[1, 1], xlabel=\"petal-length\", ylabel=\"petal-width\", title=\"2d3class-contour\")\n \n \n for (idx,cat) in enumerate(categ)\n indc = findall(x -> x == cat, byCat)\n scatter!(ax,iris[:,3][indc],iris[:,4][indc];color=(colors[idx], 0.8), markersize=probres[idx].*10, strokewidth=1, strokecolor=:black, label=\"$cat\")\n \n end\n fig\nend\nplot_res_contour()" + }, + { + "objectID": "statistics/5-2d-decision-boundary.html", + "href": "statistics/5-2d-decision-boundary.html", + "title": "5-2d-decision-boundary", "section": "", - "text": "Code\n using GLMakie, Distributions\n using FileIO" + "text": "Code\nimport MLJ:fit!,predict,predict_mode,predict_mean,machine\nusing MLJ,GLMakie,Random,DataFrames\n\niris = load_iris(); \niris =DataFrame(iris);\nnums=100\n\niris[!, :target] = [r.target == \"virginica\" ? 1.0 : 0.0 for r in eachrow(iris)]\niris=coerce(iris, :target=> Multiclass )\ngdf=groupby(iris, :target)\nX,y=iris[:,3:4],iris[:,:target]\n\ncats=levels(y)\n\nfunction boundary_data(df,;n=nums)\n n1=n2=n\n xlow,xhigh=extrema(df[:,1])\n ylow,yhigh=extrema(df[:,2])\n tx = LinRange(xlow,xhigh,n1)\n ty = LinRange(ylow,yhigh,n2)\n x_test = mapreduce(collect, hcat, Iterators.product(tx, ty));\n x_test=MLJ.table(x_test')\n return tx,ty,x_test\nend\ntx,ty,x_test=boundary_data(X)\n\n\nusing MLJ\nLogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels\n#X, y = make_blobs(centers = 2)\nmach = fit!(machine(LogisticClassifier(), X, y))\npredict(mach, X)\nfitted_params(mach)\nprobs=predict(mach, x_test)|>Array #返回分类概率值\nprobs_res=broadcast(pdf, probs, 1.0).|>(d->round(d,digits=2))|>d->reshape(d,nums,nums) #返回概率为1.0(\"virginica\")的概率值\n\n\n[ Info: For silent loading, specify `verbosity=0`. \n[ Info: Training machine(LogisticClassifier(lambda = 2.220446049250313e-16, …), …).\n┌ Info: Solver: MLJLinearModels.LBFGS{Optim.Options{Float64, Nothing}, NamedTuple{(), Tuple{}}}\n│ optim_options: Optim.Options{Float64, Nothing}\n└ lbfgs_options: NamedTuple{(), Tuple{}} NamedTuple()\n\n\nimport MLJLinearModels ✔\n\n\n100×100 Matrix{Float64}:\n 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n ⋮ ⋮ ⋱ ⋮ \n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 … 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.0 0.01 0.01 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.0 0.01 0.01 0.01 … 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.0 0.0 0.01 0.01 0.01 0.01 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.01 0.01 0.01 0.01 0.01 0.02 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.01 0.01 0.01 0.02 0.02 0.03 1.0 1.0 1.0 1.0 1.0 1.0 1.0\n 0.01 0.01 0.02 0.02 0.03 0.04 1.0 1.0 1.0 1.0 1.0 1.0 1.0" }, { - "objectID": "statistics/3-bio-norm-dist.html#load-package", - "href": "statistics/3-bio-norm-dist.html#load-package", - "title": "3-biovariate-normal-dist", + "objectID": "statistics/5-2d-decision-boundary.html#plot-results", + "href": "statistics/5-2d-decision-boundary.html#plot-results", + "title": "5-2d-decision-boundary", + "section": "plot results", + "text": "plot results\n\n\nCode\n fig=Figure(resolution=(1600,800))\n\n function plot_res_contour()\n\n ax = Axis(fig[1, 1], xlabel=\"petal-length\", ylabel=\"petal-width\", title=\"2d2class-contour\")\n contour!(ax, tx, ty, probs_res; labels=true)\n colors = [:red, :blue]\n for i in 1:2\n scatter!(ax, gdf[i][:, 3], gdf[i][:, 4], color=(colors[i], 0.8), marker=:circle, markersize=10, strokewidth=1, strokecolor=:black, label=gdf[i][1, 5] == 1 ? \"virginica\" : \"non-virginica\")\n end\n axislegend(ax, position=:lt)\n #save(\"./imgs/iris-logreg-2d-2class-contourf.png\",fig)\n fig\n end\n\n function plot_res_contourf()\n\n ax = Axis(fig[1, 2], xlabel=\"petal-length\", ylabel=\"petal-width\", title=\"2d2class-contourf\")\n contourf!(ax, tx, ty, probs_res; levels=6, colormap=(:heat, 0.5))\n #contourf!(ax,tx,ty,yhat;levels=length(cats),colormap=(:heat,0.5))\n colors = [:red, :blue]\n for i in 1:2\n scatter!(ax, gdf[i][:, 3], gdf[i][:, 4], color=(colors[i], 0.8), marker=:circle, markersize=10, strokewidth=1, strokecolor=:black, label=gdf[i][1, 5] == 1.0 ? \"virginica\" : \"non-virginica\")\n end\n axislegend(ax, position=:lt)\n #save(\"./imgs/iris-logreg-2d-2class-contourf.png\",fig)\n fig\n end\n\n plot_res_contourf()\n plot_res_contour()" + }, + { + "objectID": "statistics/13-understand-bayes.html", + "href": "statistics/13-understand-bayes.html", + "title": "13-understanding-bayes", "section": "", - "text": "Code\n using GLMakie, Distributions\n using FileIO" + "text": "简介\n\n\n\n下面我们考察一下一枚硬币真实参数为:p(head)=0.6的硬币在不同观察结果(obs)下的似然率估计\n估计两个点p(head)=[a=0.5,b=0.8]\n在这两个点的最大似然率的估计可以表示为赔率(odds), 注意这不是概率\n可以直观的在概率密度图上观察到两者的差异\n当 观测值逐渐增加以后, 两个点的条件概率会趋近于 0, 赔率会趋近于 1 \\[L(a,b)=\\frac{p(obs|a)}{p(obs|b)}\\]" }, { - "objectID": "statistics/3-bio-norm-dist.html#plot", - "href": "statistics/3-bio-norm-dist.html#plot", - "title": "3-biovariate-normal-dist", - "section": "2. plot", - "text": "2. plot\n\n\nCode\n μ = [0, 0]\n Σ₁ = [2 1.8; 1.8 2]\n Σ₂ = [1 0; 0 3]\n Σ₃ = [1 0; 0 1]\n simga_dict = Dict(\"full\" => Σ₁, \"diagonal\" => Σ₂, \"spherical\" => Σ₃)\n fig = Figure(resolution=(1200, 800))\n\n function load_img(src)\n load(\"./imgs/$(src).png\")\n end\n\n function binormal_plot(; μ=[0, 0], Σmatrix::Dict=simga_dict, up::Int=4, n::Int=100)\n xs = ys = range(-up, up, n)\n simga_keys = keys(Σmatrix)\n for (idx, val) in enumerate(simga_keys)\n local cormatrix = Σmatrix[val]\n local img=load_img(val)\n local d = MvNormal(μ, cormatrix)\n local zs = [pdf(d, [x, y]) for x in xs, y in ys]\n local ax1 = Axis3(fig[1, idx], title=\"$(val)\")\n local ax2= Axis(fig[2, idx],aspect = DataAspect(),height=60)\n local ax3 = Axis(fig[3, idx];)\n hidedecorations!(ax2)\n hidespines!(ax2)\n surface!(ax1, xs, ys, zs)\n image!(ax2, rotr90(img))\n contour!(ax3, xs, ys, zs, levels=10)\n end\n end\n\n binormal_plot()\n fig" + "objectID": "statistics/13-understand-bayes.html#load-package", + "href": "statistics/13-understand-bayes.html#load-package", + "title": "13-understanding-bayes", + "section": "1. load package", + "text": "1. load package\n\n using Distributions,GLMakie" + }, + { + "objectID": "statistics/13-understand-bayes.html#define-likelihood-function", + "href": "statistics/13-understand-bayes.html#define-likelihood-function", + "title": "13-understanding-bayes", + "section": "2. define likelihood function", + "text": "2. define likelihood function\n\nxc=[0.5,0.8]\nfunction plot_likelihood(;p=0.6,n=10,s=6)\n probrange=range(0.0,1.0,50)\n fig=Figure(resolution=(500,300))\n ax=Axis(fig[1,1])\n fun(p)=Binomial(n,p)|>d->pdf(d,s)\n data=Float64[fun(p) for p in probrange]\n \n\n pdf1=fun(xc[1]);pdf2=fun(xc[2])\n lines!(ax,probrange,data,label=\"n=$(n),success=$(s)\")\n scatter!(ax,xc,[pdf1,pdf2];color=[:red,:green],markersize=20)\n lines!(ax,[xc[1],xc[2],xc[2]],[pdf1,pdf1,pdf2];linestyle=:dot)\n axislegend(halign =:left, valign =:top)\n fig\nend\n\nplot_likelihood (generic function with 1 method)" + }, + { + "objectID": "statistics/13-understand-bayes.html#plot-results", + "href": "statistics/13-understand-bayes.html#plot-results", + "title": "13-understanding-bayes", + "section": "3. plot results", + "text": "3. plot results\n\n3.1 10 次硬币试验,6 次正面\n两个点的高度相除就是基于实验观察的赔率\n\nplot_likelihood()\n\n\n\n\n\n\n3.2 100次试验 60 次正面\n\nplot_likelihood(;p=0.6,n=100,s=60)\n\n\n\n\n\n\n3.3 500次试验,300 次正面\n\nplot_likelihood(;p=0.6,n=500,s=300)" }, { "objectID": "statistics/10-dirichlet-dist.html", @@ -256,7 +284,7 @@ "href": "statistics/12-likelihood-coin-toss.html#次试验不同正面观察结果下的概率密度图", "title": "12-likelihood of coin toss", "section": "4. 10次试验,不同正面观察结果下的概率密度图", - "text": "4. 10次试验,不同正面观察结果下的概率密度图\n\n\nCode\nprobrange=range(0.0,1.0,50)\nsuccess_arr=0:10 #不同正面次数\nfig=Figure(resolution=(600,600))\nfor i in success_arr\n params=CoinToss(success=i)\n local row= div(i,3)\n local col=rem(i,3)\n local ax=Axis(fig[row+1,col])\n \n local likelihood=input_data(params)\n max_likelihood=likelihood.(probrange)|>maximum\n x=undef\n for i in probrange\n if likelihood(i)==max_likelihood\n x=round(i,digits=2)\n \n end\n end\n lines!(ax,probrange,likelihood.(probrange),label=\"head=$(i)\")\n vlines!(ax,[x],label=\"p=$(x)\")\n ax.title=\"head=$(i),p=$(x)\"\n \nend\nfig" + "text": "4. 10次试验,不同正面观察结果下的概率密度图\n\n\nCode\nprobrange=range(0.0,1.0,50)\nsuccess_arr=0:10 #不同正面次数\nfig=Figure(resolution=(600,600))\nfor i in success_arr\n params=CoinToss(success=i)\n local row= div(i,3)\n local col=rem(i,3)\n local ax=Axis(fig[row+1,col])\n \n local likelihood=input_data(params)\n max_likelihood=likelihood.(probrange)|>maximum\n x=undef\n for i in probrange\n if likelihood(i)==max_likelihood\n x=round(i,digits=2)\n \n end\n end\n lines!(ax,probrange,likelihood.(probrange),label=\"head=$(i)\")\n vlines!(ax,[x],label=\"p=$(x)\")\n ax.title=\"head=$(i),p=$(x)\"\n \nend\nfig\n\n\n\n\n\n\n5. 最大似然率估计的正式方法\nDistributions.jl提供了fit_mle函数\n\n\nCode\n@doc(fit_mle)\n\n\nfit_mle(D, x)\nFit a distribution of type D to a given data set x.\n\nFor univariate distribution, x can be an array of arbitrary size.\nFor multivariate distribution, x should be a matrix, where each column is a sample.\n\nfit_mle(D, x, w)\nFit a distribution of type D to a weighted data set x, with weights given by w.\nHere, w should be an array with length n, where n is the number of samples contained in x.\nfit_mle(::Type{<:Beta}, x::AbstractArray{T})\nMaximum Likelihood Estimate of Beta Distribution via Newton’s Method\nfit_mle(::Type{<:Weibull}, x::AbstractArray{<:Real}; \nalpha0::Real = 1, maxiter::Int = 1000, tol::Real = 1e-16)\nCompute the maximum likelihood estimate of the Weibull distribution with Newton’s method.\n\n\n\n\nCode\nfor i in 1:10\n@info \"total=10,success=$i mle\"=> fit_mle(Binomial, 10,[i])\nend\n\n\n[ Info: \"total=10,success=1 mle\" => Binomial{Float64}(n=10, p=0.1)\n[ Info: \"total=10,success=2 mle\" => Binomial{Float64}(n=10, p=0.2)\n[ Info: \"total=10,success=3 mle\" => Binomial{Float64}(n=10, p=0.3)\n[ Info: \"total=10,success=4 mle\" => Binomial{Float64}(n=10, p=0.4)\n[ Info: \"total=10,success=5 mle\" => Binomial{Float64}(n=10, p=0.5)\n[ Info: \"total=10,success=6 mle\" => Binomial{Float64}(n=10, p=0.6)\n[ Info: \"total=10,success=7 mle\" => Binomial{Float64}(n=10, p=0.7)\n[ Info: \"total=10,success=8 mle\" => Binomial{Float64}(n=10, p=0.8)\n[ Info: \"total=10,success=9 mle\" => Binomial{Float64}(n=10, p=0.9)\n[ Info: \"total=10,success=10 mle\" => Binomial{Float64}(n=10, p=1.0)" }, { "objectID": "statistics/7-mvnormal-sampling.html", diff --git a/docs/statistics/1-normaldist.html b/docs/statistics/1-normaldist.html index 2037a58..7bc7bdf 100644 --- a/docs/statistics/1-normaldist.html +++ b/docs/statistics/1-normaldist.html @@ -207,6 +207,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/10-dirichlet-dist.html b/docs/statistics/10-dirichlet-dist.html index 5e006d7..7ee9b13 100644 --- a/docs/statistics/10-dirichlet-dist.html +++ b/docs/statistics/10-dirichlet-dist.html @@ -209,6 +209,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/12-likelihood-coin-toss.html b/docs/statistics/12-likelihood-coin-toss.html index b2eaf10..2c236fe 100644 --- a/docs/statistics/12-likelihood-coin-toss.html +++ b/docs/statistics/12-likelihood-coin-toss.html @@ -209,6 +209,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • @@ -302,7 +308,10 @@

    On this page

  • 1. load package
  • 2. 基本概率模型:二项式分布
  • 3. 硬币试验的不同观测值概率模型定义
    • -
  • 4. 10次试验,不同正面观察结果下的概率密度图
    • +
  • 4. 10次试验,不同正面观察结果下的概率密度图 +
    • @@ -365,7 +374,7 @@

    1. load package

    using Distributions,GLMakie,Random,Pipe
 Random.seed!(123)
    -
    +
    TaskLocalRNG()
    @@ -377,7 +386,7 @@

    2. 基 Code 
    @doc(Binomial)
    -
    +
    Binomial(n,p)

    A Binomial distribution characterizes the number of successes in a sequence of independent trials. It has two parameters: n, the number of trials, and p, the probability of success in an individual trial, with the distribution:

    \[ @@ -413,7 +422,7 @@

    return (p::Real)-> @pipe Binomial(total,p)|>pdf(_,success) end

    -
    +
    input_data (generic function with 1 method)
    @@ -448,12 +457,58 @@

    end fig

    -
    +

    +
    +

    5. 最大似然率估计的正式方法

    +

    Distributions.jl提供了fit_mle函数

    +
    +
    +Code +
    @doc(fit_mle)
    +
    +
    +
    fit_mle(D, x)
    +

    Fit a distribution of type D to a given data set x.

    +
      +
    • For univariate distribution, x can be an array of arbitrary size.
      • +
    • For multivariate distribution, x should be a matrix, where each column is a sample.
      • +
    +
    fit_mle(D, x, w)
    +

    Fit a distribution of type D to a weighted data set x, with weights given by w.

    +

    Here, w should be an array with length n, where n is the number of samples contained in x.

    +
    fit_mle(::Type{<:Beta}, x::AbstractArray{T})
    +

    Maximum Likelihood Estimate of Beta Distribution via Newton’s Method

    +
    fit_mle(::Type{<:Weibull}, x::AbstractArray{<:Real}; 
+alpha0::Real = 1, maxiter::Int = 1000, tol::Real = 1e-16)
    +

    Compute the maximum likelihood estimate of the Weibull distribution with Newton’s method.

    +
    +
    +
    +
    +Code +
    for i in 1:10
+@info "total=10,success=$i mle"=> fit_mle(Binomial, 10,[i])
+end
    +
    +
    +
    [ Info: "total=10,success=1 mle" => Binomial{Float64}(n=10, p=0.1)
+[ Info: "total=10,success=2 mle" => Binomial{Float64}(n=10, p=0.2)
+[ Info: "total=10,success=3 mle" => Binomial{Float64}(n=10, p=0.3)
+[ Info: "total=10,success=4 mle" => Binomial{Float64}(n=10, p=0.4)
+[ Info: "total=10,success=5 mle" => Binomial{Float64}(n=10, p=0.5)
+[ Info: "total=10,success=6 mle" => Binomial{Float64}(n=10, p=0.6)
+[ Info: "total=10,success=7 mle" => Binomial{Float64}(n=10, p=0.7)
+[ Info: "total=10,success=8 mle" => Binomial{Float64}(n=10, p=0.8)
+[ Info: "total=10,success=9 mle" => Binomial{Float64}(n=10, p=0.9)
+[ Info: "total=10,success=10 mle" => Binomial{Float64}(n=10, p=1.0)
    +
    +
    +
    diff --git a/docs/statistics/13-understand-bayes.html b/docs/statistics/13-understand-bayes.html new file mode 100644 index 0000000..fc7bebe --- /dev/null +++ b/docs/statistics/13-understand-bayes.html @@ -0,0 +1,665 @@ + + + + + + + + + + +GLMakie Gallery - 13-understanding-bayes + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    +
    + +
    + +
    + + +
    + + + +
    + +
    +
    +

    13-understanding-bayes

    +
    + + + +
    + +
    +
    Author
    +
    +

    math4mads

    +
    +
    + + + +
    + + +
    + +
    +
    +
    + +
    +
    +简介 +
    +
    +
    +

    下面我们考察一下一枚硬币真实参数为:p(head)=0.6的硬币在不同观察结果(obs)下的似然率估计

    +

    估计两个点p(head)=[a=0.5,b=0.8]

    +

    在这两个点的最大似然率的估计可以表示为赔率(odds), 注意这不是概率

    +

    可以直观的在概率密度图上观察到两者的差异

    +

    当 观测值逐渐增加以后, 两个点的条件概率会趋近于 0, 赔率会趋近于 1 \[L(a,b)=\frac{p(obs|a)}{p(obs|b)}\]

    +
    +
    +
    +

    1. load package

    +
    +
     using Distributions,GLMakie
    +
    +
    +
    +

    2. define likelihood function

    +
    +
    xc=[0.5,0.8]
+function plot_likelihood(;p=0.6,n=10,s=6)
+    probrange=range(0.0,1.0,50)
+    fig=Figure(resolution=(500,300))
+    ax=Axis(fig[1,1])
+    fun(p)=Binomial(n,p)|>d->pdf(d,s)
+    data=Float64[fun(p) for p in probrange]
+    
+
+    pdf1=fun(xc[1]);pdf2=fun(xc[2])
+    lines!(ax,probrange,data,label="n=$(n),success=$(s)")
+    scatter!(ax,xc,[pdf1,pdf2];color=[:red,:green],markersize=20)
+    lines!(ax,[xc[1],xc[2],xc[2]],[pdf1,pdf1,pdf2];linestyle=:dot)
+    axislegend(halign =:left, valign =:top)
+    fig
+end
    +
    +
    plot_likelihood (generic function with 1 method)
    +
    +
    +
    +
    +
    +

    3. plot results

    +
    +

    3.1 10 次硬币试验,6 次正面

    +

    两个点的高度相除就是基于实验观察的赔率

    +
    +
    plot_likelihood()
    +
    +

    +
    +
    +
    +
    +

    3.2 100次试验 60 次正面

    +
    +
    plot_likelihood(;p=0.6,n=100,s=60)
    +
    +

    +
    +
    +
    +
    +

    3.3 500次试验,300 次正面

    +
    +
    plot_likelihood(;p=0.6,n=500,s=300)
    +
    +

    +
    +
    + + +
    +
    + +
    + +
    + + + + \ No newline at end of file diff --git a/docs/statistics/13-understand-bayes_files/figure-html/cell-4-output-1.png b/docs/statistics/13-understand-bayes_files/figure-html/cell-4-output-1.png new file mode 100644 index 0000000000000000000000000000000000000000..3b8db3d9abf141bc6c148e4bcf5b8e865b6bbb1a GIT binary patch literal 18956 zcmch<1yojHn=bq!NT(nm4JzFr-6bd?2oh3)BB6A5DIJQWgaU#zh?KN+cc-Lum(;oO zoB8IOnSah&=d5*p*HU>m``!DAyYK6If*vW~$H64SL?9423i2{)2n5nT{AESI23LNJ z4X40=C@&?IBoTE9 zR!mNQevJH`P%iv_QgaInq}J!$@KdOqv7TO%)MxL-{(%8mRsXxB354)`4vqq^vw7;m z$*zz5uUFwUk92e}CDQu5;U|rr`g)|8k328K8G{Q^kZJw#a~dSjX&YBL@ukpw658}hrA3eAp{ekI)uNw2$ z^!2t!FY)rO`D^)pUgRXj-WV@&yF5RV%pFM-w8KV4*6~>t8t`Fxv!=4RUiIjaG4Ii? z)7NYsJ&uCBd>)OFc=OPu1fi(<=lNMY974_xwhTQlCtL(3+9g`6VgE#Q5SOk~=HXM!mPH4F01hs7g%9Tt^;M!iHu+#ivXRO`G&z1D&u!hrk%|549? za^IBZ(bLs14m<-lkJF0imSX+ibE4esLrf`9wF9$-hYpAHiiu@Q7 z5c3s+=Y9gX)Em{T!E+3L^|Y$%R? zCG)B_X8K7j+(7b2z$oORKj;aaXWCKeG{rXjysbB6YILpVJk5XxuX84fK6oxurBqcO zCGtsx4{un1e}Alj=ZPPl*!K3egx|Ms-@?fG>FMdsc@t=*-`a0YQx0$Te|h+k@AVJ+ z6*WV{**Gp^cMp$eL%Z$f`ZJ#0tVVsaz4Ds(s4#t&O^&b@1UhGuOi+nM#^h)TBy{*u z`*{yHr%ne{#67m_(M%s{Xx#HWJ5IYFarnFMwy3D6n3(6zTvLz89(L|XU)ueaF!GmH zR@C_Tg`hq-=5Kpq3O)pMsITG8M^iJ&6h4I*^=>F|>meZRLg~(SXFC=O5 z$dON2QW+Bz>6=M6C11bZXpd%zdjG!Ea>|@L-a2Hf_R{n6;`}4O)#gN5LA_xd_lx%X z)E>Onvl27m}?;dPSkd#f7 z+swnH_jY$*y?z~leQUR!r9l2m|5T+Dhe3_|ub*!$EiDnDa`VSK^J!^mGp=iF2qq>b zEE0ARht&^?^lMXFMdLy5-x0S&zBnbDg%Pc}zMBi;IyX?GNMbsV#ANsFtGRUZ@hfDj6p~b4?_K)u6r>H3JP*@+3YQK>Esc6 z*VbP8UdQS*;mI9bB^(2#5-BMj$(i!l3EZYQvP|)A?3_Dt{32JLtv6&FIUHOc99Zr} zaiOWK94__wc%=5?-ivmksJ8Eh^}c2U8I&T*SeliBM5caM`LpTChEN$We`A zxHVl}-c(_?On3V>lI{CwV&2F)_p>B+Dx+C%YXS zB*es?J53}o$DAj5KV0^fKxdd)SW?_J1(TOG^K|anF0{f@og5##|5?d|<#%8C5_=69 z`Dk-W*kM&&Q*(q$)b)S!tJ%%&WYHXT=7yg?!>PrSQ&P5OXSW9*t0>&Lp_S&Ri8XU^ z^d|WNcMlbXC%GYQKI2`_1~sp%Lw|r(c-)h?dvJt z@~K0R-mRIftt~tW=4j<)k&41X{5l_GG;6eY#ojr)nMl+;gHJpanmrH@vm#umf5eQO7Z#KgpD=lQT4 z%>vK8F8dltjL4_pz-%F0$%JIzBmnqIIPp`oD$UgsrBANdfd*rdIgDk-2`UCE+u zu!%M{R}zDaPfq?ZCBoLT&im4d?D>UMn^ zA7|r9^rjfrh258&i1UHhQxZCNQEi>LxVSDZF1}_y$x-bwt5Bzj}tKc;R1u1tj7#^ zco9j$PN~WNUl;XlZFe)2K3WSWDUX1_uCA_jc6KHtBuEF~3JMETi@NStEGK^ycDl0X z&t%)pLkGNsQZ&@6Ep_b+qcFP~_NhEER+=u457a^{`Fd3``u}ihVc~*%L*FC`5dWmk zY=agiR}5ZY(r4LaGD(K+UB}|k(2%JiCuyfiab6yqB|pePRaI4w*y*d`TOC~vpQTKE2#OC-Dir8i5VOg=5ewZyRp757e;o^s1bcq zGObVk&W)!}Cyr*k#5pj+guxLo&tZ%mZMMn`82nQL7yX1tgm8T2w>JT}+TnvuBA zOYVr-R&^mMtx1)#_jJ`MEiL^TUVl5M76C4vw1$R;n%a#`n%Iu1$w_Mx%r`%$rbKVw zM#jO!#AIe>K074^P>;3 zv1yV@N=jh86zLx`5-lz+nm&5`I4Ua2*unyo+Su6m@}$RW&XJDSb6IR{YY!+}s)s2Be~* z@tjw|-OJpE~Xy0_O&`25vb4Gc!W2pW|0njW3UqFL*o!XND#&-sDfy zYK`i-3#luj@%j_iU5<}us-+uL&Qg; zhU*_sbW3CPCvf6pf$Kj$cD_=Iv#giNE+$hy9JZ?Vy0JK_B9oHN`|eYRQ3!H_m#hQ( z`m!@_=T`uAIF+an4UOSYwgy;3%Yh6fPtS|p#g2xC1_z;}q$E8hr7m%h0&r(HZ{F;R z=V9RHzHvaWl1xEDas(z7tbR>usu38yzP>&J0)i+8mGF=dH?XQL5!42%s$5pnRlk4# z4h%#GQwd>I)%H)U#~F(69(;k2O2B1z0aA+fjSa!@kcfz*qa!B=2c?G(Co1e-zIv6L zmUg?vW@l#yY+ZhSzL}XB7>v-cu$w5)I^(#&ur^bE*Lkh3KIn~vtnao_0%D-lYWJToRf3{0s;^iZcbLfe!gpa?uUx=XS)GUHDX?XgM(x6 zBbTw9q-1@rb}=~B=m#%O4>y*3QaIkpg9YISU-P0P=Jj;dSXE79|kM2M${D_h0*RNkJ zn)!(^d+^>es{2Pr5JI_BI&L7`ciR}QHW-2<98|w2Ih34+byY5V&dHYxQ7Y3e%SoX2 zWZ0;APlOHYyute;S=-pu>seEjLF5G5=>>uCy?ghlsHj3nIH>q7DRO^~M~!`azb8gX z)IY4`XTsXq%*}XY7@uqa(PouOyU|)ppC>jSD08zmHS4!vBiyJ<|}VwM6jB|tprh@ z<>@tYXWlgD`U{pR_w^T{XlB<*<$3T435%aKz-D6gufJAmd+I@Vnlya@H}pw4k1^@5 zub%f07e7_oZI-F@`F;&M{g$p^Rm%UNhE!#}=hNYp}6- zs(w~v6=BlpLWka%ef6k%l7RfAR%UH}_Hhli`yoRUGtZuUb zRB80S8@q6Jy4TjtIy=)9qse=Oc+OSH^tv zwazue4Z@wLJD-bZb_d4lm=)LNss-<{izpc#dc>fmOnCqzbazU-#JUvrpi z4C+o2o`1*4SoeWxv=pa|`{U3pP4TXB^Rp`aXC1Q({xLl>l7e4~M`j>pM zfzrgW#6)5CX3D?yqzj7!^SG{KI}@HK1#-He=>vV17E#|W9e$73@pU_bAUHfLoLm{m zJX#-ntf6s}l2X*L4k3xZeLr8Xvh-z7a%dN=S>q!!w?Nk0*cC*aMTm@a(5GFL7McXbo3|3x8#?R}oh7oVpJg7o z_!{FYWdvg@Xa9O3gjgD<>w84^RivIPkp8hOD#Oq=HzmUE^*d8W0`=rsQkkeMI=?7$ z=XhN2ic)#r59rJ)*4U|S>hJ0m?o}zG)#*OHL-qhAnn3`y?|F9XY&w~K2C`v3+YOw% zzU_n*ex}?S&#l$M@)VESjcXo{Y|XeXY^Uo_)%`zw84kyQ4=)lb7DS;dO(X{U`}^uy z>LJ8zWDxpl8hAAAiCC6d&D4N6VEE21E=I@3@_SvJMTq?|?N3*rNxzkyEoH-Xffn#S zD>CVgx$s5KA1?JnpN3Srn(X#h#_Z!$2ET-rZ{8!9a2OA64>1tJR-C=*TKqZP)XD43 zmVn0KE5wi)3p0o33MZ$J=I^ez?^CDFZ`RDKs}>`vXXbsaI{wm7hob+PpZr0vC={GX1}J2^Ga^B6(adyYXYy9UC_{y#L>a97$hi7z9y! z&E>@@GAg##x934ppS6ljx*$p12dCNG+&o7tmu7Qub^_KKYE8Cwb}URxE|6n^4oqdL zP=nvc@1Jaa9@(L)AO19`a17^3;N^W+0^OSet+Y`#idn&NZM^kg8HSD)o;xPW;CLH9 z&}dZ5z>pKos$pO#p1)&`pm>~3XTq(-xvFp=Di}}6po%;jr*@Mo%FI_)C3~x_ZCxOU z4AJr-J1}06Fw|AeLq*x7MPU#jkfAQy`8>gJFlyD!mRndk6}67*3$qdZ07_c2`b<>R zn9YrQh;Swi0-kt!Qg4}Azj5p{dybJUwZ-D3|IM8T>3@LGeg^bP3vtwgP z(9JBr>l~lXH{(ny+hAYbe~+sm)GLV4VPoKYbBV$m`P=Bu57NPTl!~8T*U04Uxo+rT zuiE(P{t`#RfnNl%qaFxI=5MEXqu$U9l_BwSO{CF0R*^ssVwXz|krY*qWI1E_7B2Xc+c|Gh@j*933++Ja&G ze=-C_|8r7Sn?L(zNy;RUUApq|@jZL?3{o%pC-y3epTR7Me%=uf;y*&R78LRZwzG$`)D~+!}K+?SJ*MF zZEZy*B^p`krn9rN)^m**e)SMD_R?^2TL0N!F||I~UBnYVE7TL6oO>=hNypDWQeXc$ zHT7CdPHrym2r_NV{M;Oo#mFubQ@OC6jLbnGfsyR>+=m9#4&*Nvzxi8KS8IeWZ}$eu zNh{sz$!ZG5{C*8w6eNqb6%;W%q}VrtP+|-V{TYMH{rXHB^*+SJ`R7reNW10mqQ{IdS@h3|0KVQE16%i*Y!Yp<% z`t?bJp|^u7Ch5&3fi3ZsONGjv@(5x)0_wwC2)}yoCw0Au7@l~=P%bUs`WwvA?XUPg zD$?_El7(Tu@rAVi85F5NaofAPJT6X_GZo|K78d*`aIyoQLQr|1NY>5GO{@DA-^3s? z?a|9RZ{K?}1{V_@N%@I$vOE`fJ!zU87xNKATwH}}Jq|_MF3gt!Sa+q3+yx`)X=6Kb zva=yAc_#QqiZ+;zwcZc|Qj58(t^FFCl$1BJ_C<}!vdm#^81f>~lZ7aVIQMsUkZ9TJ z=gcb$P;1LL(&<;vihf{`AmpW_dNc7{5B+s<*A_G0MFcrsp3M%;$KRzCDRCB%DavOI zGUihxyw*QT%1S^&a&~OP&(Dt-DKJn`R`%SiSQZFxU<(coo`xE)YMLCIZh0gZ zs+5lJck<{N$_}#rX%<^`BF^mW??m2mX_9N4?9CoyeWB3l5OOx(c@&sh#QTi-F`8Lr z!SyRqGKJy}gy~|(v&irkC>FfIB1!040#!z4WMiWgv`tsed|K=20W}<`gats-=Cr9h zNyo=m`Z!aWCS5yNJVsmJBC&rgMB0Mr_g4ynl^;9y-5pcz5Vx7SQZRg-;F!Df;L6Zq zNa7Q->FobmE5M@?A@z=m8sMLCCnO^Bca4va59LCvtgLKqZo=9%ZX`R;Vc_?(>itly zSdwwn^|!dMBK-Xc;ka9_*t_bb=feR-w9e=P`m2*T>1Lc-tx%Y5q4w%*Z_m!lb8vFf zWRCut*=Gt#wW<{q<5a}Yc{gd^$xFx=f1lGbc~dX&rI?jVrU;?8Gw}!w9zu**| z8N*7SCX`f;I2rYF60uMW{rI50_q#x3M!keQi(}Czqm8!pcdCC5sdOS;ZJWhM>h~(b zGA}%84AnodB!n{#-@I$n`PNKIf{YxQHc!jl$iyM(W(j!)o`8AJw%i_qfu^b0S0&P+Kro z4?_fkZlq&y;`i&YsdVOYC1lzkUaM#*v_AidtjU+$y|HQ>LH~0{ZP}7W$q$xlMar`Y z(#zXvy2g|DXq>hO#ze-NU!0DuSsiTeeShrzofq5x`-2Wb9vr#7y`(CMO(VspWf~0M zN|UND+Kt6bF;k1ykL~N5m2Mpp%NrZhq#s;00NI(`o`0|itwdd5|8YF(uvY|q#bv7` zk=J*wktLz`X|W}kOC}&-0eW^D75)0o+!;6`E2A*a;z}(a^JT8r>XMSr`OCHND8+vH z$1v$jJ{ftab$>0r>PSY4-}@R?98ul%Cw+YspMe(x7sL<0zdcvzou1yygu!~kg#zeo zzP zc@&3J*+dzPeiTx4vwP6{(oAV3Bb^)UYN1+~wx($Q?m%AjXQ4af#5irLFYw=02EX6j z-Aia$NA;lg?b{s+G~1p(C_-`-StAUKq2acK&baA!uOyE#o7>$6xH{D)r` zUqMy|ISge7>&p7;fMz+Kd{EH{d&_DZPn*HMfFz0ejv@pVsh>!(aN z^hoDOq|L@>bg2D9(6^&=kBCe?VJDE-?0}84P=0PD?sopg65I;KmC7`Z4-Ll0 zGMG4S7rHTbUSIG8IoMh?F56MxOb<`S4Dv=}2zHQFpeTNpfPW?6fh7|zMYFxpsRi>t zU%$X!@Km^UY5%m01X3oiFV+%;gj8|=MI{0?a(IG3Hl_4l&C~7E;4QWBu|s>Pik9i+ z)$OGJst-m(hk7*?%l%yuXm3U*Lv;A==+^#`c`_c>j<|o93cuxHQl$?s4ysgc ztTs_Bb=P*OCTs(;8Fs_%eFJLc8PCl>Py?swB1NVhUOTdoy{vh(ZJC0o)u+rLYL)aRt-jNm*6QNNjo(EvJXXS z48;#iQ>rhdG^wrJ$RnEx$>-Ly6;5dd%d9@#rb6lA)r;}{JIX<)mT%7$oimbz<OsHduF=CvAk!nz8ZE) z|DQMDo{k+|<+y1=?=*5chW$?I*6%MpRD5y9d+8HA0ir z;juza{}^!W$LP014Bt+#3pyMP`S9MNcKZnB4Dpl>ugch3eKs^R6H@dd_K%Nu*@@NN zNy%BaDdZoMfw%i3bJ!ii51ygU13&CVOC}bdWqV=SUo-n8y=u68SIfl7K(CgDA(DQ0 zN|R(G2P79p)PM^2(#!T39+wT#l5*;Q=zz|8js1LG z^iPhr{#lFI3!l*))edm%Cmj#}JQY2&f55~__n)^wgEMF@*f<_3yuDmrv{jmyyWe48 zC2@v*%K1kgIV){&bH9I!l2jb#y2r?e^l}i~wfjprv@^w(by37cw#}+f4hswwe$5OgH6%D*8x{ z+lla*GyT@VW{tUEgXoR_vH@(mt6{-lTIcC>Cj=&ZMKDDIz=_R|q>OtlE&;Tn{Qke} zz;|nam1qvdC$G&S?f%}!vELsUJBv&_QvZszN$^kM*ZI~tw|by(C4AradVNKNU|vz5 zyDe0fbZfL%wmvmh&q*Yp!a&FVq$jQ|mP84;+R{Xf7NVa6_QOjsl-`{vWJQyEVg8UxBJy$yVvosbmskg{=AT~(4AwJR8v-^xb{KG@l6ZK zD+kUtoH0;g#?K<|mK|1SLlco8GYf{Oz)CzeQgyOHC7|}Fjqi_hI^|T9E7|(4cz4KQ zhO%dMajTcVmmuARbhLTJPt**`gRaLb+p+)Q7PwyvyPexBkG-|&Q}aPHqpDhe^#@J( z{M7xTa~jFD1bd};MUFt3N-$VZ!Rr3(gy3>FBzJFqaX_dldMw)?M+`TDp` z&8)oZ&bN0253&&bi6TomA+|%Xvcfe;LJRV@y=WyF6m8Se>~Cm1&X8Jm30Pu#S2ywiLU7t3JE_58s>` zC(9+qHK3VEacdE?pS8rMbYpR=@O=G%pZaKk`A#|cf7yu-GX*y{${%s&sCb9uVW(@) z34{g?M(D-vr*yo5kyw-bx1D(4ZWNRDK!~pyC;iqM;cwE}PveR2g?T5hJJtWEJ)rn& zUAx>q=aj`J;mGqWefl7FYqe_i-`&8)f4Bji@S82)7Oo5|ntU-otJLMtuljo71}(0TYVie z33M~sCwEM4Da|$FQqi#a$jr{s=^BhV9?Hw9a(JOlD)^h`5-XL^D0%xt#?U~(1K~-a z{79)KStKzXU!&Ot)k@5VD*X%Ga>D;m1K+Limw@iN9Bv~_U)=uw7#rR3&(XETL=u}p>n!uK;+x>zvXC-L z%cg(Er6#Ot<*o-o>-Epr^{O@`b?Cj&UR)Zsq;xM`R%g~0=Yb`K$mkZ7+jBluIX+jCiA-e4)V?4IA;v4pnd!6lX{A?Gkr&hh{JPc-H%H4s6ZEW82 z2P1!1#=@ix7r236mjwUJ49Yp6caQ2U&B@(osv}m zO#tlVI&DlVTMh*Lc8283RGvbSF;?>B9pXm@WAEROOiPDcQls_rl7h{jBe$xufeZ># zLLmVGcYuy~oleX+ZcT5jt%V4~$WjFDOvb*w0MrUiqtu%6(sLqBFG1EsKyt)zG%5`7 zx7tjMtqr&x3o?E==J-#V_VCW#yXX4@D$pglC1`sGp`y|Qm1RpPkOIW0c>_>Z2aqT( zFD=DHM_)%l@y36+5If@O#9HiC_!wjs--3&|9URo;jG9$wkRvk9fA&vWP?9GEfDiy6 zzAi5>LnFxJckgHD@sOA89UO#3M4)X$_Uwa?l_wr;S5df{;wODQ8M%nK`i%1DLFvc$ zcK6F17Ar1q3~rFv+_(!FTVZ1~izbu+C@Crb(d{c}2*4W|8Ukd?BwZoZh5BpJ$)BdY zyCPP6AD;=@Ax9E#%6*{~<0>}zByHuDH1?|bwh`Y)q_{4VNV(MayV4E1#8l<3LYchn z8aUp9j!1@lANw%l$m?uV8GlM)$1t+67$5XX&__Ylg za}OOT(i5^mPx;FuMnjCZ_g^*!5_tUCuV;h4J$>{8nA$VJ@cCD(KHr~Lytqf+txt-- zF7VAyd|psQNGQ*~zmFqIs z-)1|IP~*PLmWduD*c7ivJy~Hc>?MNx`-vO zCx|8V^z@*kXKMPD2AyWXM~@|2Cpk{=j_!XnBKL6t(w!# zr&iU`@7{fWlH;^7e$V}QThb2|>8cKRIQ9}6>!?uWit zXTOgk_%c1rXIZ!(60pmbU@)EVKi+*joT7cHx)qTef7Yn zB*SMZ*&46u33q>sDo3zI^7%pIWfs?PRpr)H_h!=h`nd!ZHZpBj-TT`sinAGxhu4bq zAF8nCW#>pC69_yBxqPRcP;uBfenjm>y-?4Izx$O^fnypVz<^`;%w9_KWYuJ->$diN zYdlATG0~h;Rv$tMDoEQSJ}}PZG>ojjx)UYr;;U8Q+wS{wRccZSF+G@C=O4XTuHoq7 z_~o-2+hacw1gaEO7Q_b^&P#s)?|N?kAab}7g_X%q0mTdx+h0UTuK0LyYh%1@x$^g@uJCOTi37>lZLMdDLf=ej|A(&@7;zI z>TS(GCkKQANPnjgI7rrIaN)@|s&Jqvj2HRD(7ZCMHvL>pPcy7?Z_qg$hZ<|QjuhQR z##1D!?6kj&BjAl-b~;U-1V@U#m4l1LrcXsVmS9xR*1egTY3cqFymzKyB=eRp!Mry zr87sXzj1{U1K-`{3wL%ewficO&27rjHN(@3~9c%kr4)A^M|JrXN- zLDo*`90Y9{>5Y7Y!prs_>FfWz#OX?NnH?H!ncYY$bc#plUt|)iqX=gG zl}Wq~evDfrSj?r0rkU@#@ofVD;mszTZ~~bNU} z)4Pvhe|%xT8KUO_!bgJ#kx~Xj@bTHKM%1kW@y!-HngzTa3Vyz3kHOU_@#L zsf!O`GI4|Tc5_FWkT8fdx8Gdw1a?8=la$WskB-JCt@716xKbDGk=C(IwH7f9Cb$@O z@3^QIYT$hz*cUY2zisV80Y#mmtbuKLPiH$PJAsxUQ*@KP>Z`n>-nzJqj1Z$(p~2z*d02(A(|M7(2Cebz3y(v zH7&iISoxrazl)&%ex&65<52ba=pUMuR1GrkWLR})OctCyM*~b-b)~IeXi{Vijkw8cQMOlgU2z)Hf7~J_?tG7_(*C^M z^lD?sq}^0DW;tUW6a;&mc6XtcIZ2KxMb-u!2?HwSl-_A$(>N-B-W#}agMD^Thxn16 z=;V6p_xC1geeEt`VfNoc-^$!h>?qL9WFk(+@z5Q0Jm@Tp(C}V9z11pe2>$I>*W78X zzEG5AjpM7QP^`i4|B4cR`t*L!=NZb&O|9XH>T}6WBxG8=b(r9BW2cu0fq=l^b$?zT z24Bhm+265u((xjC(FB+!ugDX4OS3~iD{A$~h&Ojl@Wk=^c~QdM@+x}(;9)>jbJ!Hg z8Cm;Nq7@mU#*4uutUa3*L~etp}n` z)$UWHzoEp1Zf3;ccyO_!IK7Nk;J6)1?2NmxNs3-!@Gv61Kp^H`nO83a$293HC&G*d zzgI;qUamr=2rZz(dy*1lcoFItSHq+F>@{y=fWKRKCDG?Ald@<9hq%^{+DKdVIPr}t z68jTsFKz3`o*V+XLW<5iGN4J4$>V@)W|=QpKlTF;l>PtSNSmgeSj}<(In*DUz6;z3 z_`B9&p|(@j%*%lB5WTg!o*mYL!}# zTL>j-y!?IeALfe2pKpeTq^0T_Y;9smQN64235O50q6va@{r8{-ynu@hpc=?}rjJ*2 z`3+7%0&fB(LWGVYK$CiJ%tq`IUhKg`ckSw49+IXPnA#&Q_Q}4C*8$&AL%OTacw#PO7i=B1N005b6 zBWJZpv{G}$H>q6Oxqu0WoS_=0`uc_XhQy z*ln8hu|20p-7a*$;}0^g3{x?~#HS%>O%jiESQwn|tOKcK`2g>a8;6(AvV=Sq0m({2 zr0N;u;UMD;KMzc&z*t#-Dco&hdl{Zi!a+vy=PI=7PqnTSuD%xjBBpaXnrtz9XK3ux zP3TSN+ggqdb}?yUvcxvaUgayaY>`kisP)?3)@WHiiy;$d4&I^AVL_Ghmevf(c&(7d z8&^Y4d5^R4mj8Dx(K82VH}y=~6or(Xf9TdX%_Ei1IojGg;jk37U?F@GVRjl9A&PGH zOMi1IDpFO%>gWq(&ZJH_01xjE)*uZ{p)UV3cER$wyY%54-Iaw*Uo3uyX?Kc-bGq5M zpwz~a5+2-})Axrna+q1L?A360h%g7SqUqL-XRX9v=2r-M^`8j>f*ubq*WHW-4V5w+ zN43w<&uH=@w23q+dy@B_W@|V0{~~}?9p4$(Un^E033sMZR`HLR3)~-_mUc$*!BmM< zo+awB4i5|uEE4F*9vaI1N(+g==vC_VeFxp+y^m_im$WaMelv0iU)RAw?N2yMa5aX} zK!mb&h~ZoukdhbbgP+4EFxbS|DdJ$)4r{t|`?lUa+wZv8=?<+dIbkZDT7{Rv1PfjAr+bIhyuKr-9 zwYk@v@h{KZKbQZo>?V|i-ddqPcH@t@P9r0ewz-*tb7y*vI?0ov!a}aTO>;qHS`_C4 zE`#U~T{BUKgTeDwN4qKgV-K61wQi#tvE|2(Y*szmhA^f+=MC_C?RJd;Mt*9 zG(p<~+~>uP#j3?b57mz+#}S^A5N8((VW+>z;fW{YdC12Aft#3kTDj1;5>JlRY@mE~ ziZX^&D>a5gk)AN$pCbQIfJ5R0i}U4J)fI1ukoMFEb4JD}5WUmWIqx_*>uPfW0P=TC zg-gRsD6xX`{-B2BA1Qv<_>3ju1j;dPmm`Br;^`DWboTXU5aLNUhJEX81MsM)ODr^B zym%3SOZgU;k~ykv`#w;;fEoq47a2Z2Ct$GwEgp&XMCY;7q3g*n3)z5t;=3&A~8t9L8px|UH*K25g*qSkPp5eT9d~mn>7$^QR!UQ5ktkgDt|<4>Ia2} zs^<&=fOfLT_&%^EW$5ZR*!6z!LRvOHKF;@gd>u$h0G|L*5H1D=Qju<6pZkf)3nP=A zdKTHB_$E!v9HJXFVZ|E&n1kT($Dh#-Ai}kJs49#JS0ls6+FF)_Me3f;88&kKISP`` zTT2#Jw+VWUxU}Z}9IQdYN)C9|xu%d8++)CrpjS>5cm+^;>9@0vDISN~0B;8<;rZEF zzJ9fk^=uvBe;Zp`{LXF`GsX;)7C+49n92R2^~jU9cCJ6HeeYzkM6NicJ0mJ~Dm3|m zviduCv#dHl;Gh0DU3FC4biPyV46MUijs-6#R-R@(9;hX%e|-q73O0PdyU-vnE-q$J zNoI{|Lx>y-tb{mFkAPl^0A?Lf2TlOlbCZ~O)m|vl)r6o# z?y*Z=QTmUoN=I!KUqSAnyQhn}WYf^kY#@`7rD^9FZj!yJ^|N(OugNL%hzrTn>6U;H z4>PD&xbCx)x3scZU3bX!k5P8ZXO1fa9-7dDgfJW1`jP9zgg{n0H0ol!>u&&bCOkYG z2Jw|4wnJK4x+h8aC3ifEq=X-^3;-=4bN_ytq@rRcFtxrF7dr@1>aY@c2VUNw4MxZW zdNWG-n{tPk5Xcpqa(}#HlipbCOa(G3Zyd8l*kVcZu8uKGBDY!Vi7{(lRtHN?7e_(0 zs`gVY4b3OUx%oUukv!9~Y)2(!w+5R98?|q)r|%Y(*_1{zLG(k57xGc=qAWZafC+bA z#9E3*9YJ>z$gBlCK{{_5m2MvmjX~N0C0Zoh2QQz-hwgAUKfJC_w;j4n+Osq1_3NE+ zp$yP1zEP2lPD}tf3&7cLq$*># zpmuzDKBeoN8#(WCMnF(TXoIENczHuM&>KB46Pj`e@r*Bo zUNx5|$br5HATSBCX$rjXk2Y^YlTfiKZ;(mgwWG@tG;vPV5gY2dnVFd_x7wK<;fVQT zJG9-}!37ncYs?}3-^8)+Bq6}`gvahP{dwZb<&C`Q%1V`%cbct)41MV=Xxv=sdbprBZz zAgvzdcP;9fHY+aECHq}gsjh@IfbpJ2KCvy8#!j;5dbKk8d}`w@WKvS&xH?2LqDxPZ zLrC*IetA(5m%@8|^zcibvWc*dpIqu%T3cI7EJkJ5o?BiwQdaUt{2fg0m+u){;mv$XOL7W{ z^dYDUyuVRo8Q!j-#A5|34+e$);wYbDZWZ^=&PJ$fTt9KbiuK=*E8LIc?EL9(Q^vw+ zc-%%HFl?{R1;CeTy<}06d?+t3|L~zOuyTQnGTS#g8W$526BmaL^nWjL7)O$DS$g^cFe!?)SWSvu06P8t zj_Bd2F#;70?csw5xlb&lrTvsM-GIBt4<;7)0btTHGBSXkW&ZbOI6|OxCgL`}FU{b{ zh!s?xcUq`5HSXQL`vM?QWmAsiB5wl%x=eUP0F4v`(EYM$*PM)uW>elL8X8cIy#aiD zON?;D-8*+EM4SsVGh0^db@HAlyq|r>9k28B>94!!KbN{brnoG{uZ`-tee53r~nC@ag+z>VN4dGS#4i;D-0iW(amy(0~OZ*HPSzJGsB%J-R} z;it?@24I8vJXBT&mcid|09`NQGHbsuv`WfzQV7~EZ2xKi zq+j`60XSE|p8^-wn>OHI$0yYke3J!EHIS5#B)WOiCvsLnMa7`V=(^?2@JB_!vBbY2 zA}cFPbH}+~)H=r+ri>v-0c$0Niw0ATY>-5!y|ZAgcn`>Q8yg$SnS0^FoAL?@^>>go_@O;^)Z_+DBHkVAI; zD$*!GRnEW2S+EN3NRAH0PF{Rq^1xO;~Q2rCriKtgD;MbY5~s_u<^im#uuS@+yeAU78VwuUmxmK zIx+=+gBo*ob~e20rHRSU834%kBns}0y-ZnLTA~v7bO%lB`gt$S2l6p|BwR* z7ne%dkyTLe1gJ1jRtLrqp!43;=I6634LGikB2mNn6LUQ##SRV*Ng^(`z$gFmvB2H5Rd z1?hNr2*m&{`R?}pWZ)Iah7cn;FTA^dby5c4>FnVUi1&AV@C(2*6EY}+{`AE2ph?Ke zHUTXMsFm~m{c6XjKmbomO>OV!0JN(_>jFe{e3p{}e0*SI>7Ae-1N4Hca03WUz!n4I zC3LPpF*|yC?nhcHQs(950kajLf!l!*5h@QKTD*S!rWV8)FoM3m*PN$-Vtzk@s@nCy z94L_JJLXnaPL)#-U-<>FcI++ph&!&wO<8IZT^C~qfMi_UF)$IG937d1-vft;gal2> z@{z7?m1!TH(f~9kI59LdypDeb2-+?A2Lv?S2E?o(=l7D5#pUI{%LL7guU{Df#ERFf zrK>v$TL+{FSpt1HDFBh7l7d?R7A&CQSuu3b-b)0IfPTHuhP8RC37jG?ugmT2!c@lX zzOk_tYk{k7OYmH8Z)BfoZb8A-OC)QNV-18WDMR$%50IT3NtS%mp))?Qu! z8Wea7evx%6_Cm+ttGtynrMF=l6%-U^JdToJOk`o(zc)5Yv{>OGg*tf%2yrnBKmZ=c zEW+C2&)ue_tzWS}JUVJ=Y^+Pm$$=#x50g6%2@Ld(fYVtN-os~fjKow%5DU$GU7r!a zK5uVZ8yV601{kjbB-#Z`Ku?Nz%EymbpMWWumS()z9=#o#kYHY-^HfJCI3(oPm9qgd zSb4bsyzlmHKQ2cg^8xn_`H{u}D1p50hBwjBF*9gZ(}4JI876dcX~QT+V;d}U?j+UhDu8+=8gRSC|dp&}=r zD17GEWaP5ncOQ>Z$bPoY2Po?}*x2jS)xyAA0K)6lp3BI|U6EA5PChL(L|UFSDTX0J zMuybOy8(>Da_C^q?!b8h`sLQo0O1X+POzVNqWitYJz`i=XF;juz_bN=Z8n?`VaU0E zaF7Rc2T~(~SLqa5xJa+2ehiuN4i62%^~24|t?}@{2tA(IhSn7n4UR+(`Xh8PD m`|mFGfAF;bx0jb|k~e2hMyTWJDA51{y2d=U;PjJ73~bog4rMWPxC+h31?*>xk^?-(suuC zJ5!0n(&XJqZ;1?lQK{I6zTHilT3C$1vGYHqa2c6act~V$l|%`GNo3-PFyQU%rvLwM zPy8%4?{M81O-N3T`1v{SmlG|n2wGmAGcAsciho`nZ7ZYG+0oVw0m~k{iiwc)sEp*{ zG8-Lpb8{6HB2u!5=5ybAm0S$Q)(2Wz>D{@kS6ZDI=G{GNjdb67q$va*ZcP>Cv6D9T zD9(`@bMdId=w8~7RC+BFMn3Y{DI# z!?ht|*bEiYkbA*RGg$cK%Lk{zlkT^~#s3cI7}uUGCbve>xa~Ak%sFz`Z2!2sRkIi4 zcRn}%ynCYBZR^ixjT_$EjRW~Q`FWmGK1ZMPPDh@lFDxzj9Zk6Pe^k~lwd`Fhswa}c zWn{sVy;%pB>r56dGWcHn;RB+f)T&>*K=0t!ogxG8o#tgG(pNR-yB$l-Xp~~^+b@O; zkA_WS2LAk^BnUCAb~W@l7-(*8Mqm?A^7x(~!wk><`s%)w#$}v0Gg16mXjtu@#?8hi zog$w_C1g7+d3mB*KgHd_Cn_K$T#M9;9Yqs z6UW4OmEcv6bs24I^HW}zM_O80y^ge@D#&Ql$c*^M94$_?rt4!4^{M>Nx(m6rg zqU5jPVv=(iAiWOPRV7ZCxVeexKDK7W|2t)?jz~uV(Ns8(Q7XC5k+ewOUQ06nrWq1g zozGQ*+lJ4lS7&rxYZr@*>1!4{x4PMiyCrTd)kqNOl1WwZYiXH&dpBWEA$Ju|olH(k ziz$%t>d455Z+-pe&+YaVw>BFBFocqoi-xkUC3k-8Uje@mLv(O(ke4T#+Le0Ot-QQk zC0R&0s}~W=v;wy5z1z-Ge{s52QzJ&wSF(}!LfAbp%ef%_Or(f9>0q@hiK0%-ePngx z*@AwNzpN9I3MI>gD9L^DN!R zwUy3|!G3SC)3v(RZEMoLLSIWOvEwBXS7O8*5)Ol6OvU$K37P#wT|DTDy*Bpr@xnfq zi)+=Z+CFN6b6dtyoL~Li=crJooBks_v*O!z7JE4X@MePjRoqFM4IGaSynhL{p!oILWK2^{ue<(;A$=cCOD z;c>?xjt<{?zu)h8!PDIEfAi+eQkvIh*WE3#qw(j0)XJ=6)I25)&%4E+2_;J_c5BZO ztrt4fUqspl%sC>n-V%24aGtrz>BbLP^#r2U;zMR;=Fz0r7I-Ur2Z!#Cj*jkbg%sBC z%Ok%_94Qh^MC3QgzbmRG7M-Gw)BY3oH;|4DrUXvlMn)^0<|b>sl9Q6acc1RB+)^UN zUD2k5L zp)qF1cb*v%)w({ax~dRwlfbXs*QF_ydGQwmhbTqJPETGQ zThw->c%f5xO3Y=srzbeUJ1@_8O~!37zVfw+z%NUPky~&d-380Y})C3yec^PGK)913%wO1)|anSBs}HaZ!5F~me7~21ZUMzRj^6D@DZl!~}S?SMRT38SV!(UiffWvnY%P2cKI?7eOUL1PE z*Bw@ZItqJYJV1LYE#z>l>(>(U95+>dc}W0{1o6G6q;t}3-*{s&H^wR)dNbwVGTljn zHWC-drqOfE%1Hu#7v8BVp*I#>tLsl!?*ujSn6;6(sH>?#j4Tgt0tTfBEy( zHDTLfM~JO{e-DPHz$k)a?7lm$&+@m42+o1HAR- zMD>w_-`UvJTk@asJO{LN%%44zJl}2#^IAyIne^W4ynaXfUFV`AQ?c{HuS8?tkp<_4 z`T5`AcwwCpu>QXi?il-?m?w55#K-fyZH_C&Gc#^N+%W#`Y;I_nI$@tK@_Z3&v3*a> zzd*m@;-uRTHc7a95BydOELwbODxt~RV1eJoUOM>H&b@Og_@rC>$n4p(IfdHYUv~g} ziCv!Wp*hkDzgYbPPP4SE>}Aso1aQ<^(eat#BC>kj4+eCm9GeK?%?W=rjGj&(tquDF z1^I=?+Zx!zyALppE`*Tg@NmpQSXdYZuOh#}G_qusl-gGZ^2f%;qU3Lh z2Fqjs#2XqMG=e>dWBPkglxvypz#oqN6f!SZ#JEPX{*@AIi;Mv8$3) zUFmiwN)47_Oj}CiHtup7e@S6<49M1QtemhRn(prE>Z)UKozEFS{F~8$#0ztBUd|yS zngap2HEkVLl$8OB4$HLUd*9P&UjZ*qyd|=F!-xh329kvxP$&5fn+9TjX-V>QIV&L{ zK{kY-k@rnPLc!y2PvJ0COt{kcx5u-@&7abc50$;0uGtr#=Z(;MyS z`3l=JJNp#TYq{SVoQJ5)s}$-_oYlGs{|x5Lxygt zLu&`_Vnm@BbiXb5yZ1tg=wyT`!dvQ))6^;UE|jImY|TWqx}x;GxlGdos`$S8}Z~*cTR60^o|CbYyKR`_s3LEST4&-2S z;tYVm`G$iZ{s(Y1s(z>O-Ahd~_O%Ct5aj_`xvcb=cf9?(s#6cvi1MOa-tN(e1~E2qjDW$a1Kb--aogBJ{GwDW%b*;VSAxj@lH> z;A-bzZv||Kz(hxL+4ii>WW67kJ+!8!MTi}=JiT{Ps#gGO#zGW22M0^9%4G#y=a->1 z93onz&*?tIDAHDKW8Ism(^Q!Ye#3?-{(nPGFii0=R{!Q zZ@CPso;HWW`OePHCZG|r9sHtujrk#^p!FYEEeM<6Ai4oCdeRiq2@d^7FrJl_6?g4E zX*5IlO@ta*?!)(dFIF`XLdkX)=O{Y&7xl@c|ey%YNrhe_tQuR7AoayYsKEv223>f$%>1^$802$*G}oK&Q&X zLsu;=G)P7jnP!j9dNXQjX#I_fV(#69{D|&XM+aJ4@fX%#XIu5@uU`kKFe#z)fa`zH z|7;WdA9yVXfD14iA^UNpWzUE2;9ycjo_D8-Lx}T(Lzb(YT<$d2^5@T=qPo*73Ow&0 z-ja{}I%50&{rhE;GQHbM&eapl@$EB?ZxR!;WbN$i1gJZI|Gp<&;Zc95pXoqVaCV6* z<%#Fs+9^iCg#?i}t_U`^^R9rfFt=e9S?HX{1-IV`esqCOk&%dq2r42&a6|*=teh-F zNFxLhlbDhR7S{nn^KwrHV%+{{Dc#R^@%=Mv0cu2}35maUz7Hzg@8?R2gI|I4kq*$T zd7`1Vx_U7rW*yx1wEn0~x;$T!3|pJ3_hTX*sqxs$Ni@15r8VHOw+QLx$=;H}&DW4` z8+)$gAT+g7g&or1y0foX%ycih{VtyOD+@Dj-ncw&4Tl&Di8!R&cB7>$WkW_ZRKkGr zA^U*&*Z(~nOBS;8+$bB8Z&6M`I<>_z0?x1~@fayK&v=9zS`EnTx*RI$27g6OfL*4? z>EG9<4j#(D0E|48_{PPS-{siWQAaVsG++WnW#x_jT;i${6Q;vvw4rGNjc}u&^k-;k^4bujZL}4d8-@V(MYn|B_jsl;2#1(jU>QRERGH0? zO;Ig2!<{>8C>EjIR@cbL$k+D*pSz}Kco=Cj$bopRIRBf!TT|-=WKF;yfYUI7^9R-d zi*P&_lb@dt1MG4#(AFMAd6W9fvwAEntgX5;d+?nMxQP4r?=!2URGKsfKIYJV^yn*? z4FExxge%RRB{`Qc~^|R?cRG(vxC7fei-A36)!b zlSxiazWelN*zE^z#b944jNG39^K&+JS;%_d$?01r>w~u(9{artD5^z~Y4@o9(opRU z5H6MG1;I4yHRFB%2 zJKE!}V|Hw8Y|CBj`rcFYy9Rs=gp&G2WT}R^Bf4X*I}HLP&IANZPPJY~@+_IGd@YE_ zv)dX3*g@L)ax4{-p7}#-KJ%{?08sYz_osWVeg@x8;S@p68Q%O9!v4|OU$Ak&1+;*J zW;w_oXlPvXKR-Wb0Y(ZyUb5(myu3V>+`XC58)cAfYazP;hIua?+h&iCH; z@0zjX@pKS^`Qf{-jx#?X1_IrrPS#vgQ}Zm%<0AkGumXGvE&~YN03GGjes^?u!O1&4 z+`wlL7pH6Gu^+F5h<3QsdZT_dFP}+CbnQzCIk%B54<#by#h`At#MydQj6{{^fpx8s zB&5aFn^pJ|{hfQ_gn4;K_fd>sbN)tqA{ zE0m>mx#&{wnWY4x!94&N8Z6FYJMJ-Al0h9my3 z-aMv4Wfl=!RiAQl1Wrl?4^=+ zvbK`P1HLMeO1Jg$#Qww!^j87O=q56{sh4urQJ@9sLK-Cxn{u3Q{5PC81h|Zn;$lb4 zmd6$r%%a<3F)Hg$$Ut!-1)kO!+!iElBV{P2l|a&;Kbo{r+WE@L%5v(s_nsQl+tS57 zDP9p!{5eg~ziKydrhgvPCr9EeK#=VPCvkK9`^m`(!-Te6$x*dz+UKDD*B{WVLL3Fk ztLuCq>kg(R#8@fc-P@C{DzU?v2vBq}i^jrZFI`E&-sqw}+K(eZT`D-%xb4}Qelh2L z<6Jbn`@J@vIm-{$Mm4Vk?Hv763isspaZ_&rj@|>lhuU)Y0^Or&-r0MvtyNhVw-xsL z$KC=;*cZ=SdD15;EkDG3wvs?eSYmJFV%H+wFc*zc!Gt0BoV1+YB{XZbjTS%p!KDC3 z>k~J_u{K+NToH5(#`(mrB`v8|Z%>i?NG=a;!z#5=-7##L=#T)C-3P1=hZgV6gMK<1 zaD+cz?@1gimov?4c+N2I2*cLsq6pXie%@5_8@+X=(y`%#%iM#pJYtVPD|x(N+Gs{L zOMxU5q6g}RAv+#rv($oDa?_a|KGlgowj7$lBYRosg+r;0LBH zf6O>)XlPhlv!kf@{YSm)b4af4ylKq&#|s=KH6`Wn*8u?}&c70b&p3+TJ&AKfC(s*g zwK}dUd`b7>?CEDLS$7ayCE^x?sB{R#j1=a+?X zq1@|xXe<_ws=hT-Tx`a&wB8JgRY5T;Z|{;6x}GU4zH|CjXg$4HxQM=?K(^0%tOWX&;nf>*-L# z(e^0e?%<#)BVds&)Y=Kn8rYNqGYfVBvq+Jtb1cnk78dB>*rrtp;V8*Am(hwK{%;Ll zCe>t|9-6>v$N#kI(Ifj-+>9e5`ULeY!GZv7zsjDE(@eY%o+e4_kIdpQrt&atLcZV) znL7Ach;fCrhoh+Dq%0@b{D8Z*6hOxM`Z`3hr++@FEkDxHsQ@qu4(Nw{g|X|1IdCdm zxGPlbdfv35V)=pb_nMA5N(>fu<({(7iI9-C#(sc!h05L`_wV}s+u*zN0*W&MU_WLn zL|d4dIr{^jdmgz7p?Mqv4`ddQ`Q3f;0}}y>RsBh)gBX3c?z2=aBclu8Sr#F3R`2&c zkX9sqG&wKW@iauYr9J+cmg`_?8P-j0IZQ?;f>)Myc80}f?QkibLIbw&S>OshCOwvF z;d^f5PB=tIe>7A*#L|c*jz4_;+B3HU=yw5C$`DyfX`=K3R5Us^hDzRnz?5R(0U!h^ zf*}OcOxD0yGo8VCW&Zm@azEj=8cE(6@~S;~VGK;sj8`E!gCdi0x$S#!m-3O6PAyc= zY=^(jz*LuzM0DrAmpC`Pd93=ffBcX^c-CG3P5~m-$iu_K+IkgI-JdIBG{OOKs7i#e_!J)8@|xycG#aXjzt|GO!$E09pIv9rsBN(W)r zLae8cX*J>b&;(hTHPol>-~X2~2S`Ew!%lVu)-;Ne_bH@WUPqe)YX?ARxQM1QWaBvE z{1_;>xTjv%$Fwp=wUu9bIQ1&!hlzgACub*;2{l0d9jgAJSACX1L z`=yh`MT@khtlnjobnjy80rsHAY3J|2fWGB6mhe4yci9Dmxo3QNN+Bu&)SrdL$BqS= zz#lVT$2JC0%zr@a9C-W*`wG_ZK!$0|b(?)bPnne~OU~aNciVbnm7aQ8>ltuiw5{CN z(yAEa#0_K(SA4`k(qk#0h}rZxGmqW>F$Xy=rUalUK-BxdB7OY02y8|TaDzg2qh_e_ zS_)X0nSHocG1lpn!how7o0lg;i*huTj3&4Esl{v?QHk)3%y*U z=RI^}m}v7ZIeB>ufVZH4!jQHk5-od^lJp7R`={R?^HAQyfg+Dm910&5b3|k=@3{9g zF72T~1Q+YZ4|F`+s34{nqm+8`WCkGlN0rpmjfyFv8xpj1^3DIFWC1Y$^#?~JEde&L zHvfinFH7y;Dq~7F;C~B~+yn66REPnc0N8n7g0LVZ-7_x%f5Uepv%k87m4m@Ya@|Uduv`j z5l!S;ElwVRg^orNxBe^1+qG6`;6ro!TN^FcA>mS-TNb*Qj0oAA*hI8o0TO_vz^KqI zZEeS)+Iq5--c5-yeRKx-m#t8;LL}u*ufh$^$6ux5m=KVHvO=Yioqg6r0>~cN{2{0Y zMpZ)A^S*O&eccvyzvNn7i*{Dt6k8ANoOT-P{ovolxH5>4duHb5ohc%FuUOI)A3W%2 zZ`UoePUv0&Y_)~}&Y+}AG-BxlavA(U9)MgAxW3N=P)y}B{{f3ec_SJEpe)}*^zIq6juQ4XlhPJC*89ByB>s%j zM^k?Ftl{7gZUC8fw!Z>Ye-aQwH&u#;p3{Z`xNi{$5>Qoa=O<)eDRFVTJ+j0TK=_3v z^FLD&pf&_N3FVXi{(e-xm+VkS70gr(v9plVapi-6O8`d+$CYv4MB=aiZdS~J0;vF> ziC@ZK2L&v7{Ct&p3@eESP|u0N*n0=um^gh5M~7YkIh}i2I+R|bgeR8|Jg>|3T9?7e z{p5tB_jd*YlOe;vC-8{o5D0$q(o%>3$ABaXHI01VX+aysOccCNnZPL#Et{q#P4V?t zX8(r*+ns<|R0LyF8=Mz>@H9l}vlh;+|5XbhLlJwAv4ctH`?;FdbtF-Ct0_bfQa~%3 zS!Pfw5p*utm8`ef$IN1W6|nCPk*<&;>htx^I)ey~2=prXG-asfB-KdA|CI6+!4zaBi(pm zxm;yI#bv%?RmI(?w5R}* zzUK28ORWgMM%`Gvh2WAlVmr3Bf}VCpX4EF3Be#T{R95*l)JvsYl1Y7!7>As_myF$)fVgJHNQFw0wkO z(1uEpxJ3$g<^I7&qrAhXp28-!*ESlX=&wi1rK07qR~k`MSn==JB9z;;6HOr%v7l-| zOZ3x>#>&R|chhi-qMKJl$fQFT+?9ZkyL*BLxU+%HFlkCGF zP6Kne1N+1Guvm%I>He|+M-$VA+X)I4_nD%;^^k!y-}&TBx=S4YvIOzL!HK{uhoZK&c{Ml_0`PvhUW3+_(uzQnf`O_A}?Zl(+S zQ-)RYn@>>{DqkynXNDJ_@)mWkr^t+ui+s$(F}S!oViB7b4l85+E4fCVvqRtk@+IHH zIOz^Zff>XpVXPf*guPO+?ZQ=pGys-oFAM(LBGh6EK4oW5@);r5&3%G==|aue&$Q*R z^BIBDCyN?Lh%lMR@Tz2PPnQ!i!9BX;TiX~>Ira}VC?$I`p)4Lbk0<0VgXK!6@?9){ zcXKj?+!)!!zaUhg{Z1nbP34xam-1z+o%iT_&|d5?JE307URtGOwLMP+Mc z*E@~=n>p4GwX(j2?r-dREpu=BZC*Z%Fa6)XBxS>3Fh`V<;Lyx1JzZ&#UtC)?KoQNJ ziRjxy9VaOt{SPBFkoRokvwdaT+O)Df3?VBa_#Yo8blix%w`Y6Bmk~=S6~OZgoMPCC zZqV_ppkMElNDVPhHb0jE#FF>4zC)6)6&=1uOHn^{<6p4r-b0ovGWAVPa&;6lDIuxT zyD12;Imz*HmAJsmkkz(W^~_SlJ&Ai5_jqz0Yvna1>HOZwTMpe~`w7vQR*lDxQB2Iv zTIQH!|xo*^W+icK$BNY`BpRFE;CU_O*(gQ0v8XYz*V^JbwuNhct^H&oI)FhLOzFd89Bg0L|N zLUFerbU|iTF-9m8_n<^rbBLUylC~&_&J8u3&5f zo3TBog$}?Wn6aI=@W=IP*mw7P0f8^@!iGfOGYAzIy*Z# zIXQt{S&mIiq=;w+*bxdBH))T%o$<)BF%Sz%r2Ns$q^|($xR6FPBcnC6uH&iS^^cv) zaScq*C>{Zq`K>K4Z#EII{b*2cP{9RyGWzZVjl0r;np$#!l%&0bgUq;r0K9N&jkfkq zWHLxZ@vwtt93SWQXMi98rDgF4%4!WtEgmpEJ0Y6r8Q3{**yi3>{VVIBS(25-ASELT z)z{OjM3OigWcFGD&zF}6AoSU=)_}>ME2FWw8H~Zk(GrUk9_x8uwz6XU&-*mB>}+i# z<+CzbZ%RX=`_`N(aOMv?H+Q1r6L+6=(IQ|wB+|D4{;;#LeFsTOFQXt)Y7KBod3k_9nZo9fYGi`cg#Rl5MZJiE zp|SC56<8C_KfYd~(!%}g$2B$(-CxzxN#{$-7X%6VUDWcN4LGkx#Z?|!*M2n4(SARXlV{JAk`zyOEsh#BkY zRRPHj9D*ZlZraATdVwt!fjH+LIud$cjK@yS;zp=y7*Z=#ZS_Pvoo`>m&8Q%mDoj@Azxr>3V z10lQ#W!=K6Rdy7m3`L&-O92o~T>nQa-2hq|v~V!hGmuEwPSyErjFv&|p0In5fk~Iv zlF(fiMHS=*peGlOP*__MkG`OeO$%m}?$1?Keel2wxOF)67F{VGKoRP5kb_!FIaTyJ zqSXlji-r<7BpP{If$#;=M!>0;P%Ci^=I<8cr;ZAYWl|qmAK;T{y*Vck`%L3}W*IH9 z-Z&2pwQl1gPDol)TbTSZlcb~{FgwR;C{@otZ!3{}^CJ~g3O0bRn^&%*5DIWe&gu!m zUi=T^bS;pCNJ71875&M1td%E#6KZ+~?@$&=NS4@}vYT**>@X z(qnrH6xeuF#LP^P!W6nWTBrk0?6N{D;$x#gTav-sDCucruIS^pI}C0u?ivm%d| z-HM_$#I(5lkPyFLAnhPBav4_+(kTkxXAs8*_GX^V;gGsi=B@4zt6Sad3R-2q7H+&X z$6n$6Hz1Gr?u;SPeZ3SRfE^IrlK_aII3Ff8ElmNh5{_b{JJY4FGl!4|qq`Q|imE8G zNU3I~EkOHn+29VyWbBp&IXZ6Y2%P=HDY6j!RN|ObKV^tNHW`|6`QZczE?rx|OZIB* z+kkxn)J(ER!7ttfY?}KXJ$d)Ea8%wwvB;c390PUR&$c^ww23xaj!A4^R{)byyD@a< z`TvhE*0kHaZFb2D?td7qMUE^iZWdMA$KJ9%@X>a<+yXJ#0KI3VeHGvxH{*$*(uz<%`rWvkr~YBA}En}PcX z)D)i>XX&sM+eSq8tJf*-yu68{|F@Y)T#>i|R*&lb+9+SA`di%=E4J1zmhXYG!&pHn zlvp@uqkJ@`jMjT-N(pjvG_wW^`vGD}dt!rm;Mr^k!{*elCIt2;D3SI@smqlRXFo%2 z{?}q=hZ1f9hED(f1XjoI+TVdk3vUrXCyrZ3U3ESA)b-VlPYQNQ%N5%y)jN6ybzl2G zU@5u6K(A@*jlF#mSKP^`Vv6OO(@c?PgAXPfYt{+_<$e9`RN2SNA#@UQ&6eFS_gV8G z^$X&E$j}}QtMin|BXFT}^4rT1WUYCUd9C&*%bkFK)4!%MId8Kcg8h-o$@BSg8BuL{ zO*?)c25OsrZ=MQ8H5@rDJ9+Y30-)0DUDrS(VrmCmsPfRtP(MWJS9u$~>PVZO@=# z{mN$VSJ^_>Ce*>81UZT`I)K;TEIC`<# zOP8I@RKlSRM>38MwixrxV0jKrnTR&lM1Wy$HdRh-sMvpVuK_upn~zr!8C3Xt=~$_s zNF@QmoGeM3Acw}}0h5W=iRacFoKMw`Yjj%X3Nd>(InwX=VPNX69C2I?%%(NJc*L#k z8tpradq=#;iDf+bqdyf#D!oBW+_2d@ryuK5mntPxFsfpnhE5jtuYa;ihMLp=xHPO& zlG?AbU&|8eQR-^p_6kw@SjgCHl!T3oS8Kv>XtJZQZ)}_GoSmmi5giiBiR)nsP0BZfv(jFlp-im zsEcst0SHU4%JI??O#yU7DGH-W zh(i62s;3J6Me+J6Nx)K2OUus5L1Vph2U7knARbEf>{!Wub$9!%e`^KRHik@9!Y^`T zSuaPEq^lF14FXwtnWg1+mLKt$lz5qLl1Fr=OV;P)m_yYD1g21l?}KV6())NDBug`n zpn;k~RbfG%W9-<71-v4^#c%iHnNX;?u_0hg#z4(hy;uWY=CPU@FUZ`WtnCXZ6{wtG zmm6NxqiFlef*-2czTTl?wo-^Xqc%fL@E(AkBcOWE_tGyxyY~nIwQvZdIw-!?s*ec; z#DG8X+5EoL^^P}a4#jhTT3`)`o_dIRPrnTVMiqmbeKT!ErxbNr2AQWL3N;UZn@f4# zg~m9cU}dMDr(jw^(Ez;)EZc67)tU4uRw55;s5%U!cCR83-H5!w3Lkm zS1~S_2=)ji5d%?JyA$nf7u49GCh85~8|unG!>$Hl6UyMEs^Bo70^H`UkF#s1hc~Oc z^1_zwktT{&D=)c;!N2i8Bvs^ls(4RXIf_ZiwVIC3+1WWMIl0!Th)HP~=twboRgn6N zW`F$XOgq%m!(;LL9aP`4g1VRVzde|REANSgk)$Pf6HdyV0h>%UT7~Nq2d7DZgAaaJLT5}47tI1c{7ekvMZnS!ZGEt z)Y6t-gWwu7NE+&yAA%Tg8%=nqqcTC&OdchVB`t;1@-8nA3MQ$xP-sVKT@N?L(tS@? zfO>b?1xBmD;5%t5$T?7iW=;+kUiiJ6P*Wq5i=LBeb9yD81uw#64XkP_>NsPcO#_JS zAvCwl3U6~@bD*?D6A^?u_jd!aZvK|o4^j*a9ejY#a9b!5T3Wr*T3RiXAifKpL#6>3 zMzu%mh{aR_><9&oBA^ypFgHj~!~x3EZbI&n#%qc#>H+pA68&5>71jB$=rOtCg6gh_ zejq@-fb5!a#0~VX3|Y0lKK{vwN+p+SB6gCrvuyLSj&spPR9Xd3u*<$80-z#WcpxSm z08rag5|o#g4%-U-1YJL@3Frz{ILISZa__otGlMe7 z0P%?peNBuUfruZUCjsT97cvxSA%>RvV;I4kcU(}=8%(< z3y_)y3DNJaE+e2Jn`RJ<6R~+u-3;|#8-808ST0BT#IZ3)$mZ|J;l#N>Vp8vS3DP;q zz2E#$@o+`bqU1DiK0(%o%Ru@{Y`;equ?e#1mSKRPPAH8eRQaJHgOi)vxXgOMV<|Q9 z^=p)%(srzT1F9YHT^7;y&d!a2d=pH#8dj4|qzC$_}$M*^isz+5uX zh1}z5fdttZ&U8AIr#)-y3JdK)?gA+u_J_Px5HuI|S_0yg_$M8Js)z_Af$;C`+qX~? z$DwyeNP(^vNt-5ULoPal+)Tebe+GoAl%m_Ch24QNEbzj;sn_rB3#C?OwfLxAYC_LW z><18GpwufMSQi zZw!CKQ244>z>u9h+&}#)&w5G|qpCzRNYiX+6jJ=H@soqwM6XS`ExJSz%$G2@(x8y&>&5D|96;z!Ud0zn_=c1%Nu) z#qXt%VU^OT+FYHdJGr7Wl+l=;jx`P{atv}0_-}USS%%Y{rk2_DzGL5}RV4lQaZ^R; zsVHGkf7J2?l1P>>lG|Q{7_ROlZwLvD6n**YW~1b-%Kq|= zw|?#}JkmN(#K3V}D`4eror?)t>iwHTRrb<#B05Rv-{Oe${Yr4(RrhVF*t-?OXHF`E1zE;SH3Fr;DDX-kq$MX-%7Z|Tn_?7gV zS1t-W2xkaC1ie=qone-puwwIdG4u0>YFUlbjR%1lQl$xqu*&;Q1zy=-+G#cZSn{|# z)|#-A*xg3l{GWE|D|Ye4oqhH7MYZE~ay(>V<)nzKI#!i0V%`Y8Hw~vI#3nAJAi&Bj zz?Ax_MQuNE%UpT!U8a0cW6&CwY=Gq*G!mzT5d9-X-FT6YpL$KUr#%}VF^t|v+%w3T z85CC^G>c}l7N8WXH2QbxSS~>SC)`{V(Iifpf1NGsF(?A_7PTVv)4A0KiJiMBdFRJ)Z{Zs*p4#gAHdsRQh($pyWo$1tRy{pB0>-%UbL6ZkYa zJyvb7M%w}V;XYbE!AbEOwZ!NHG~Ql`HmRmCv8$?M<<$#;RH$Y6;2Gqw<$Q)^%DtX` z>UZ~00KcA}yp!zz@f!utIlvBPt;Sp{40vQ=YyZDZmvm3q|C)O{5p?c>(LJzck7l6v zKj7WoIQH9-c$|*d!jjKHe-vK$x|nPF_FaULX>qh)=#bvkTr$I3=7tRqToE-bW^|WC zMP)Zt*sA1zRrx2bBiy6fMK9=#Y~)hZU83H+K;s@k+oq0xaelXEFdUh)OwG22Zd8mh zEy%of?B-r@-RXX>p!y(IjTugZ&BMW+t!=$Kvu`qOvcKTN!2E~aeVLn@DbxPRCO4$0 zOR*jIw@UlF z3-)C_YWM>BDM>)}#ox{%hh@8l|eAz3p&eH+KO zf!AU)g6RLn@0T&`6b9|QSo7IWL8+6--_;Oy_4xexdF9nxjxAJ*4IrTA^% z?_Bo1aA$ZYgFkDFj?3s~S~_}YR7g*jN6cUOAdYaNC9JOQ>g1;Wu0!K zm7;)PoF}T#`C_VmTm)vpxNDar`O?AjWe3T--0Ul9>Y z+lRz14`@xVk)ks@h*n?Te_{QQI5spEA)`oh@Nux>$cSkBk%c5(7iCq!jF0!}v2}S$ zwFvV^oTk4_9h`|0W`WhBVs!r|H+np6x$`;yBTx1(Sr%-7dRBOrflDB_)=`-+t?4m_9$6+@b!_WKk&j5%HW@ z`Y@tJD&;4#Pm#WgnE zKXI}Th%ucg2V0lwPmjAModQhXz3aJl6l9x!_mnSJz+K`54`o@|WBj^yupw%A;|^)B zNRZgGA9U_W&{N#?v1!b|%f(|NBcqC>oBLDC!rw^v^D`z%-^M~(&9~~Wvq68}S^9A* zdP%3PjypN^SCIc6-4SjX#@|76!^6WtW?x)f92^`BqXj<20!Vep@6jDkb{BwE zNa|cfHAX`j56Jq7ayxxZO)}i?muA!X5eX4T#^Q#0M^B%23#VO)+3oi~7pScncW!M} z9!t86*jS>Hg>&K3?!_{dT;Z0cT#gp&KV&Z*x8r9oH2Hfx0gpkPi*BACcR3Q;dq9-k zO7CS04L5oPuF%>Dzd*x-1_(??^=my}Kv53b2Cv_^VXUPk2t>iu<;4;7mECT2LiHGa z(Ov_^GOFbrN`*jVaOs!ZPE@^MnD>nfh#TYVsQX@i;bLWM-1$)n&(ze*F*3kK{#~Qp z#__S0U;AZqIle-K198N@649&ZW6nqY;Laj$l$!=d-B!|u1#}CFTUS4?%Ca@wKaia0|b(@a}JV)`~w+; z4)k%6W26_KzkFNoKJhaWbKzCg!CQ_5ncGX5tJkk%-Ho0jB_%Dl8yo!dC%e~j8v49n zz?wsoJ@RKuix4!wQcwsB3PN)j3p4W%Mm9Dp=qG^Irl3X>1S4HMT3hUnTQP?$wi%U> z65p8F^ZbHlYnyMwNu1@(&uy|s4h`i@eRco3B}5I`M9D{VD%_Gku8i^i>ELY;)A$5C z*6Z5lMjdDECGo+>xn^*fg%SoXYu&ki+XG=iG6zjJebB`;MazE3vRP^cR<~8NxhFKR zxkvK~yrbDIydgtFyPbwsMrtwip`cdI>$$o!Sq*yj156irIk}gt;lL{r!McO=8Y*ni zTZv+`mhcJXx@ASqdv!L zc5Mxd%5$n-;cM|;R5o?bU=+?>PTtfW2V$G@CVwv#V z;O9nraU@Zjt7X>r@864yN$r5+(R~YLSW%_8Y>S8u?UL*BlJpo2+inB|2-D-cA%4@8 zb0XW3uh3hcPgBi_wD2x|oG!g@ls{i~U;552_3fhT9&WFRw;y3!%H=f7q?@%{$POtI zF}Te!7QcBE^X3t8mf-TLooUwhvPN4Qp0dIBdTYuBFPQafKbJL*vnyFzFCeI3Drw~7 zN^uK%Sz0FKm9ZB1@l5R&B^mK%cTuAUNntSu)ZCEo9x~^se} z`D0ynKW7%%g?eqq=R}*U{<)`x=Wgli`1o5L#UQ3<(JdGhi3z=^sJ~%KFfU0Xt3Zj- z(6GG)1_0bvZY~QH0qOcN(CUsURv*5=yK$26CtV4^mhKh`JSjoUPD;+2?x-ZazWBu-^CM;cRRDJXBrgW#JdXE}kwjoU2O~9HT zK8#fVjADfi3~2qKSra%TB6=~<86k=}Fiaw$+IX9ZX?|e=dUEjt5hma@8yg#6zka>B zxv8e620Dp!(HHFN@DDqD=(ogE;y&;^=wxc1n8Q|ZAi;wR$ z^O6FZAFtBW)02{do}&u5M`&R|2QCQMqXi`DK?`Et&&;H> z23jgKTi(nc5PiP*Frv8+t_f5)4;vevt9;h$C68fv0E_cVU(VLZXU)9G$Vi~dUoad!luA54KAtAJ zpF_G~wc#i%dKdIM`HhT0IH3fE!UI!66WG%QhlA{83@SJMA>H(GVNX@b!Eh1+U?*b` z?H_vYl_81}2lW6*cNZ61KnrDfQ%U{7gNWp08V~|xKQPm^x3W^pAAkpj#6$5}Mn*VCjOibHQKd<^}U1n6{&MPFOO(sXucUvgA;Okc&QBl1D?W}nC z-t635FY2pBwTIMI^ScWjum&I^jgZemV>E;Ya^9yu!$2l}dUh5R%VlM4P0eT42CZl} zZ{2db1n5jM>kQx=*GUBMd~hZ0e=|VnfczC63?Qe5$!G{|Dn=d)*Wqa>GFA7%vlu{fS*fu34uFPm9w01#1(Zs~W+}xbFn3$t52&*ksKWloryZe*iL^Kye zcM*`vaBMO#pko*M$e=X^VsUh2WUvOFBVHU6>VZ1o8mXgNECu8h6qXhjo}HY0{^@B>;5FB)-fd;u{&Ftv0TwSo$GV#04qSsD;bL7#0A^z2^35UvDnS5N)52OI*AfFQuf$5=T8V*s8xE;Usg zTE79-L634h^oT>{OpTSyvdN+(;msSUk3ofnmo(<^GG_Ldh-fkv?0~qT3 zWX~fF8pS`<*ZbMqZ$iTzKeT7Qc|(SLg-ArT*3Hk)!%VIO_5hy_ND+j^(;%vynVGRP zg-0vE5CQ0b3qmY*CRf1R;5j6})5K}A$Ijr1Ibd-{J_qnj8&e)iz?()OaDs(F!@&L! z1&>^UZf$thkpt+T7_w`D`vy=tI;sba-qG+VflWZlG{W|g(b3pE@c1NX@v;Sw3Qcp4 zq-ezT@NH-qgDz$x63k(7=o$o3JQ%u|7_A+&LuNt|5?tcfH873GkB4FSQdQ#l11Qlq zGcz;jY=&>Mu-r-b1vvfRSdx;G;DWKSv4AEkN=t)7Lg0xq>2pXD^vdtw=_I@<8=w@y z>vd@BGCltN`TZO|l(Pai#LK$_uYqk3WN4tm6g!8<6j|HY$jZo|AV$_BR%nXtHPY36 z9}^REKPNDi0oT{p7y2}dIrv+os=mRD;l>|6Fd#uh93`)ztsQa*k38Za?|=fq`JpxVgEJm{?eJqAqQg0`M9cj|)ON z(&9AO3T26xjQn1oZlDB73d2MiM+ba3heU##Xs{H>tj)~4nW_TMLrO?WV(-|6$L7GI z(UFs*9%m2~o8Dqwh7Q7KBHC6^Ei1)MGojP_qNMj`x*;f4mwTFmS0-!HL*`PJ@SBPi zt(4a=(d-oCFr_jceW(SaiTYnEIoFUT!zhkV9n%b!!+o%HE7Np|%xjW1bV$MkE1APo zrp-&@MCKH75;RxZymc9hGoi%0hNTFS_6tg>3BcZfv>%(Yfa>A@bIZo&oIv2 zsMSulc6DJ@8kO2IP7nluh$pXeJB7#ibsc}U$>1UQj-uh?{(gFLGOcoSbQC*C5PTi9 zG^UKayiZ__G04TKIkV;OLZ{1Qvi{y)bN5KokibPl9X4bdAGFRw^fEbA&x23RgTL32g!cYElQz z0l;bVjYsW~`CzIaVu;}2#7-exK7J7qmb~z0m?JnG7u9S|0EZ(;NDy(D_!3-F$%1Hz zQ|@&kUnG@o0|r;qS*I8b2MJoy@D3cjH7X;ZZvZKdJQj0gdkhYB^P&vW60eybK#ZNe z3DAWY-XdK=9o^l_OH1(oas2>cgy+tCq-seQ(&W*% z0VGqNv0U?r6s@epD)0RgRfGBLO}^= z*a2?(3PHWpb~ZZ@F(o>hU7;g!CRU7cZFSX?&R{V7uQiJT8U{UqiK8`Vd>Uy+fEzA+ z4?Vp0$y?!`U1^`z>lbM=MNWHsfir@WDWiq8Vzt}ZAz9SV zC;nzzh}*dM$E(7v(FZ>#s^%E5^}YP?&m{U^Z(sQ7=u^9*95ma!;3%ZT$6dYLA1k@{ E3zv;>+5i9m literal 0 HcmV?d00001 diff --git a/docs/statistics/13-understand-bayes_files/figure-html/cell-6-output-1.png b/docs/statistics/13-understand-bayes_files/figure-html/cell-6-output-1.png new file mode 100644 index 0000000000000000000000000000000000000000..31049b24fcf3b68e1fc21042c5d905142f670019 GIT binary patch literal 18476 zcmb`v2RxVk`#1cxN3s+864}`+U!}4~;UbY;GBQeLR@s}7y&@wsl*q`IO?I-!m#mPG z?B_V^`u+dU>wfP0_1yRCd9K%Wxj4V4^ZbnCa~$vEeY_7L8tO`?Nf}8|DAZ{cW%*ku z6wWsM$486{SE@(9q{IL4AIPZ5piu83PaT>Nz|T02x0Ga2#h=m3@B@pf@+~zK%AFI1 z^7coe_TiHE5(?!kghH*{L!mAwqEKh;QmZwk;8v81yv%J^obQc9u5mj@*AI7}O;3BI zJ`+(_r)MR`p(f(>_AS>9!#iUhic1i1ajEy(7kvVSmRGl4b6>nae4YCi&&5ho?#hC> zf>3(-C|1$%DCb|5i;hN%dy7>~g7^;J?9$gt&c5e6wxntMT6 zS(#~DgmSvM9=gpZisqgkU7J12Jicyxi#>Xv!e*qH`KtTjUayy8mT)TPjhi=J)+dtE z(<8-_au_1pn8=BUC5_|FTkZ>gocwH$h9CMZhb-Q2N=1sTjg@zm@@4zYiuJ!UZ96;#RnlVmG8tUp1EW|AHmAdiHKV}-ge}AY+`ywYgn!2kr zikwiEnDqhrf9GqUsH~gYUQ@s+ZEc0CIaaOrg>wf7HJnPT(*|@ir{Fq5hrciMi1;qn4sMWklDRJ@WEnMv^-WXM<#ZX3b63Vvouw|~0Om1JyylG+TUjA{nw&(78 zL0zRI9kB`x@?%G$hZnqOC&#;wp0-6ho1b<0jLT|uA9=Ja$>dwU$R^bMgRg-v z6As4qk6iA^f80m2urrbQM`y?dmnUs!y>|I{C$B$)2}P3=4j6uTIQny| z_d3?DIfz23^IK2Kv@7=c+IU5}<7Y{U2|GDY&ui921{Lc)64MhEkMH)TU2O|z!0sN{ zD<14Ercb-C7Qv_cN5?k>Yqv|g1h5!#9&K&ylGkw=$f=wyndw&*@x|~lXCrOQgh4OOBVr^OBM)}TxEmNQHpnz|bqA5#{`N6L> zxV3EbVV?NC2JhMCAiUqd#vhk_T(5rf=1o&mlR_9Bzh3G5TE(;gH#h37MNjf3w(Tq} z?fl~6^`q^r{hwWeN|DTKy9cLEow|ARX1d3E6@#dS>dVV#q}9}V z@-(wPJR1L%GJW`MwWQwo?Y%}Y)F5&u(t8_IwQ@n^rvqe<_IC{4TlBznu*aP+fokI1 z6lvKQhK(qK>>p@};>siEu#wt|s;7qb9+U%?Wb^p+fhKx%rex3G!jye|eJQT}_Qqm{ zO!U5(Mqd54dvNEC?u(e1_s|1ca;T906`#l$Sr`XV~_v)@tt45m- z9?`e%r1wwO#>_iLvY(SL89A7b`R9W9aDjH%SwYWlNe>P^>#dSY zYJPUkwS*NG6pU3ozV7yO-ln9M&ij9HooA!&U_DH*PFF`KJ}F7^$@i}*w$&<>_^8*f zUpp>-=I7@p&@eQVdNfw%>x=t8y(q@V9Qk;na&x7y*0X2#$W5byiiX*9_6IiGKyMs_g+7FI_0)lWL&!+%3F*%`2HTp%*-tO{Pp7!ucOt^ zFRz>;CAHHh`Q#kL&YkJ4FOW3ZS2+1>~55?)?O3;?xn%Q8D>sJeZ2jwsJPmdG*WDf}qzZSPNUE zlmiUe9@}dE;K4s>QhM6MZd>2t;^K<*%L@t$tS2h%VT3VU|LABtPt~}Yo0}Kv6j9^j zoPIm{@Vw`5*4{1UkTxuABbcF{o?hcUyYUZ?VmQ@gi<+8b(-tcpjl;?8>g+sT^g1Rg z|9o|KLYMCYtPSjQst&Gw509wcjAmQ#`lD~2b*3E@=6F{g_(ZjuSQ7jqoU5qc(UO&T zQ72WhHafe=6RyORoHIN$bhwp%LW^$#2M%#c)y}J^!^5rKk1iXESBJzb#@qZ zh?}?1O+$VwK1n) z?B~`eKe_;VaNV3Uu64(WVqbpgGL`@4%|w~iAea}Opz%cE9YItyd+pZeE3h-s?9YyN zM)KNb?a0uWw~xj@pg!16x`Sg*NJzMXj(Cv%>=G}p&B^gmh6eYQE9c6km?=!Z&ouhI zP-bV2Ja~F-sI$F&22)&GN}Zc>Pu8(Zl@>GfBHiSir}6sg>h8(ToaLQ#Jby7va>NW?HUQdej7Dm8UFE+zSkbbs{C$z4J`Kq-)e5U zCn-LzTp*cD=%j7>RgVwmdgB$3r@>0K{<&KIX#7(4;jfDNmy(X(;`t56%B-60C%^}D zswLmn)y04Zg2RIXuYhW^2e+6aW_8QNWOlwICRa5Hcc#Gy4;;ATb7-0l{sUOA30Q_%^93sYp3AU9j!g~I(OB*tjFTYM$J}n{qa6{ z$&f`iu~(UyzQW(O_B2+TP6^__b=|z*n%KK2S9D6ShOz-xjWMjH_4&D=%Vv7yLCxTO zJm*5VzV)WpKO8NEzz!4qCs}N{ps;vl^!*~uKT^V8$A{)!@o<)R)+apxEh3xTyXXOf z)u1pXD+h=6bi~L7hyJqDlJ*@V@u$Ck+!Uyik6+N`3)Js08z|*V=z7f^${qSnm#myR z1RwQrGGZ{)LN~)dZLvHft&XY2uX_LZXsgk~^9~CsfMmKr>@7V#hOm~W*V20|`xDl_WMp;IwQn0E76dWx`$QGZj??NFep*tiH2!uX64a$ra{viE*~c`6AKA< zk3UWy*}~~v9VxjBhy&n_ZSps;g79XA3shaDC=8msS?K;ddh4~X@0v^bf|^a}75V7o ziXdp}#=}G>Z4uZ5(Ly7LT?KXgyX|@eKpn8L)h|UJ03+HWm_EAiteg_4oQd>0#Db@} z^mwWUkOty}zJGmJc&B82tbCQ?`k}Xbs`?`Mh2!JhdITE!5uTR4vAaHLBrpOe-B938 z%B4q7)+a*ri71$FYl~mMew~f2d39&~>pQ%1a66a?oAD3)>bwOeiqX|dCV4@Nm~^=@ z9jakfEn0j^2R<^Zg8cld-6;Hqo{JY(z^TIEN8ad(I?OAh0!;UlcehtZdS~rWVErgFz2ooidvDwls9wqheEeX0sUN1lzu9*F1$9`9 zh7Gtx<>fC$DJ;Y=xV4EY(OWrcDRg9ni`_|V#7$%6Hn8{Mu@E%DjRSJ6nq2R_cJLh> zw^AoytMsQoJNC8~OTeWrbjFp`?OUhE7M9did+e?oe0aDAXB69rTDpmU{g%CdK z;7mZ;1_!YPYmc^t01OUo3f!sA?b`skSh62ZR*O@Iq&k1mE9FBR+1;m`3f!~~A16`Q z_U36K+Z@-c7A-6+eyvX)#QbwD>C&TQ?oiaV-!p#iKVL3{Ibe!DiL(_j8m9!aZGQcf zSHR^wC_u_GuS-9$-GE2&1^^zc_RA|S2wX0zm`YAhpN1p>y$)gi<7)9o+uEHknpuj^ zPL6gV#O%1kWsj1-P!x_hfQ7wZYvW$~%ejDkez;2QPB>Did%;vg;ZL4N@9hMm=sgxJ_{wB*lJ$Ir7ETR;wjIclqr{y{QQZ#`UqIK{X-#E_zvg71gE@C zjvl1BeVcIXh2+QPL3bi+=zVbKXJO5c4>qoEb&CStVgT?10qNl-1BiBzsPE1AQP7<| zdzOw4wrm0(0da*{+&b1=?{VE8r?+G8RSCxD$fN(C*SWYJY|gdhDN*`EpaO{3^YY5b zvJ*TE;=@4)S32ZZZLehfP=}k5UQfXroOhk^A)(=XIR0S`$uxvaWFfi2QB4wA0L${) z?pK0jXTO=^I-q1;imTgE<7=%}965~kJwL|e(IHAE_C3Kiih*kQtZhwE5#VTQV3A-~ z;X6KW$ZZ7ce&6rURDfKjW!BhrHUKubH&s$Ae&O1){f7X9!K>Bnd@pRnp zj~mw_u&9ZrTRj;XOK{d92z%C_gm9%*eoMBA5^-J`G85SYgsX7j1z1o$qz-t5WYVE+ z*bw;(sId=^;2cf>B!1eRC>(DNIBIz?7XS`<18jEfAs=idP&6yB!RgnYxnt13_MSqLH{(x2gOY+{6#fo@ zrm?Yrkh{P%!Gsnt>lO7!Yu!nrXCb#S>!6mQ{JFKY*qeUW(2(pNq%=T(7?zraunenuyKJ2A#JE?@7OYou2IS%g;D?gvbjR7#L9d+sj*EfX5Jv4eFj&^;jG< zpFb8dZ6(V_{9dorR`*1e<6YNT0toaY^~XDqK({S(W-$syachF&L?rqL>9yLLlK3&;R67Xl8e+sH+-Pd>1#VCW*-wjjxj zkBybj(p1o_F&12@gkuj;&`?VY<}qO_vGN<)2kW6cX_gr~jG@3;b-CH>VuqYMRQBx5 z=d@0XW-_v5Q_&B+^Vqdj&;)eor;oT&n7kMx_{uzOBO@j&Fojq4@g(+ zb6C`hY0)9<2t}uq$$i7tR)Ddmly6Y1pV)#w0^4WiD?I#`FM%&HC&b+_`qeod?&9Pe z5(maozOFX`VtVlbV#3*MXmUq1mkPD@=)({47cvSiu4rg$4~oJ2^+Pjt;{jG`u7dA* z3U0`+aVC~WtGMbLBYGWNxf`%Yno1v^tN|kd38_6h(3cAi8RD~;^>_nY$<;USr6cZ$ zG4cxvDzVWAYHDkn>&3%-MLaqS2_gnfGk!YY?a*(&C9&yk>~#1@hzzT+v~)mBt(AU& z_POHPlPd#_c1@TFIW{|V(B)Qe=?}lY1x{<~+3w`IOAmFc9A5Kv<%xrd14o1FjJTs{ z_PwXq&>@qyN~?qUT3kTMR*XNpcAl|OhzxZt{U1tvF$PB^S?b?XULSfF<*%;MP-83X zLz@@3d{P7Xx1 zh|J6jYJ*4|JUo1giZ@3RpUG-xM+`a-=CH7^0INe@4uJ8y9pFU>kCkx5CUy_sj>#H1F@MR)BZIBO>QBs*e2=4MJ$R5UL7D#oo`K zfm(3nqdRE3er!@Pj=WmgQawo}yzK?|g~w+D?iu1i;uVzx_{|+?DBe58Z;Fcrcz76u z@B25;q8hdi0jQAx8HfT<;08W6LGW%JP#VbnzkkdceSD(O2{?i9G~+e|&Oq4=pj612 zr7CYG2(V-ipK5ik=%(+TRc7~bY3EDqlR=iY#D$cci?}IUB}OBtckv^{T1|@Htvfm1 zl()LxzPTsc<)AFpHmfY}eH{s{9zcLX>5>7kC zA^u-?7mso@3XobRYEiTMP_+K{Q;H?oWvB?mx>-p|+lgoxA)Q1xeiq6RLy;JK^3(4_ ze2r**Ntpc6r-0n2rBq!K2oO$FYNpkEds9}Z*eA-=gvGy1tho$`{f%f0xe_aBsuuNs zej)p{RUm`XWP0z zB5Rw|{FG?TE;h*F7nKE8Q&wgtJ4n`*6%0}Ed5bt(F$2qd-k+Bw_TjkF2AXbi-85@Sz+S@FjVUFnW>TYd zA48YWjj_wf{_o2@0wk9-8{So#wc-oJ_DQ#ELJ3D-UnoP8BS(n2Z7!BTWyTL#j5Bzw z*{=!as&86QW9|8#sA$MvqaXSA_E-XZnq&hL+_qWWhpKPGWQE`)WI+|Bj2l_*+FGSq z3_7OUHRzKs$?3mOQuwpn%qhW&lmt^eryBn=&P?iR1{K7r74JzcM(c5B7@&w`89zxA zTWQ)42U`W+q3S0_6TvA7&>k$qPlMG)mM$+CSx6gGOLWDBV||7})tfVtlmjS2S>>=$ z_8nxoFEF-cz0E5qtj-r=_80Sp1s&I!AQ|O=5b`9_R%gwb2~s(%E3I&Tismrt?@s4Iu>yaJM^r3A$s`#Pa`2Ay zuv7{JbIEG7=xFa&WFz!uVZ|9+k%xl^Pr(1k{mR}@1S<(EF7wxvE>e!XjAjvKevXYM z(40ZcQY;)9r7uYjhp6UEKZZP{!Mp!4#X0L0@izUNG@ZrM&`ZLZ2c4xL^F*xOlGz3!|?~wEP z=Q*i~)~#ZQOR~i6t)iv7u@p}Lenbh`2}w=uNSnYB&dh#H1ONNddrfcoTA_BTj2J6d z>X3Hx&x;oZW;H+c1YJlr5itfo18&1gG9%kIF}|%eioD{vf;L-z$|X(Ylc+i@+}6>l+xwORG^v)c@f1XkYQQS`zzeuT+Kbd;Q&K!2mhP^NpY<=h z(OCHz2OQGnno#;X&nqq5AF<2_(+`w6EOe^Ga4sw^9vmG*LMm{#?&L?J2w}Byo#A5M<5_8qjFC-|@s9LJd z7>dYUmLd&XRa70sVem(j$X(9tUxNFH$mkJbE}?4&@{;F}Q{RPH%3m?z+l*jxNRbT7 zt+j!TAZNNP=Css@#CXV$#K}y*ywL+Z0Gx#0lP6Ea7;7QZXLg%^0adC}NC;W7<48V< zCvB?hNnf70q&eJUIRN%|Pq{7d^tnrq#>~jfqSjNTAFMU}aj+HF~_uTdQGa+DUPW9B^NCj*BL(zY|Fohe# z9_Zvso0Fu2$@hyQ6v{RQY9Leoi>1?lv;f?6XU?z&Ebc+h3DH)HIt24leF!M0zO<_z zwsm_JLo4u!7~^D>;{(80u-TiC=M$B~{P7!XEBvCqa7FboMc<=w%Wmcu?J8|j0n7Jh zWdh6}%ie+Ae5IqVt}a9$-f|!It-NX0PKtT=;3yB0m6oli3ZN{I@q-GNh7Hu(fLOBw z0syITU#`0$!vyr<<^izZNcEi`O+LQ%BmyhlZZAv*xS2Qf0x_<7L8wM|Yo!?<<1IHa z#>eb#5m$6+$OCAj2K~+TX6lp3wV$Q7do~oYWz} zCGQ`3W`7?m=arBcnYDWr0iZgv?Y(Im`4E)-Ktus#55Zfrj_51vm8MWD0t0#ku$9=V zh-{FI%elVGM{E?_RKPu1#A7Qf0oQ5%hxUUE0EdN^^x0OA1h9iIfLG<~!fmiOggg{V zTu7A!h&|y})YHSwHYTttp!}zkgQ$G!xoAlKniPbHQfER*pgv`6m($#}|F2$E2s!LDt7 zqJf&M5u6{SR0f2uN9|PP^&ps~_WUp%_zrhw24q(k8hI zF{c>Bn?SLvhd-K=bJUVV6b8#CCMFisz3M6}>3QDiLZ*;y=sc_q&-<;#6Pk(OX0N@j zyQuF|wI2;CY@sI1KrZ3FZAm%M`RO)W4g-WhG2t}F4n8D8;L#7nCME3w9|45cr~B_9 z?L?AHOQu0MOE|A?&MpLV=B%46E99Y!(8c~pFo1EF;7v+w4cJ|ml zvcaF+WXoIJ_vMGRo2dyXZ-GD~I{*@bHhbFu)@RJ;!U_P%GDV|gR}icfbx#(FDWciE zS-}vzZ3tO>|9uqY2yN>X^T2I4+W5YnL~HLTwm@X@3Fe zRNvQ1uhD;96KXli=~@=i4FuTB#FEqy*KhWn=1uJz_n<|e{m9;~itMBV^4j<)vF^78 z9*-_Wgc;F3yl@^u@81x|nvv9>JzE}v5g^4fMN@{a{^xqZ@?To-Y{;g12myjn|2#jK z^SQsR`x&aXfX!4blMDpu_(No;Zc;H4RgHW;TloXM7_Y~jpAz5V2z~+$lidf;jBgNO6h<>=0W_)S!KpZqlhU|%4Vk3!NXZWTef(j8+>JdmC!m|euZjYNv9YcTq#?s#b zy+8G-Jh))7+AI|XrTB1)R`GwIqpeV$ZSf6bo{CCoA6HSKHX3eXiOVJGK zsNTk+w~#3%#y2>!WP=7v`TZ4Ib&F921zNlN6QtutOoM6&vb;qEvvXsJ*QQe2Po>66 zrqV}VP`nM0Mg%LPjHrXY#SL{AoTqSjqqLvWsHe!l0EvYI#?}YH{Z0GRW5U|HgnGTOPaU3 z3P5_>K8uW}P%71lx4TFvRVbIeJ)DBm>EGvEq~H==i#DRQYeJX?Rs`z*xxc+kgG%!~ z6E%%%s*nQh_he}d+o>_Jc`i&_6djKW>3^gmSVC0in{ms@67E(w- zkP|@O28O2CSjnuu;+YU0BQt|ZXab~L+5kB60nzJwB>JW&cFcGqoqptDMxs-tc!0F zvX&!r&6A2Q3zN+;`lK>2WqNa^)kkA)!tg9{ev6$U6tSVM2rMMlf2JYK%7k z`31!)v)|I!--ox;(b4%v=``b#D@5y&t5Cb9vKyx)K1u*l1=x3>6BH$}P|X88I0n^9 z6fG_78z4#khk^re744;lhK4o3q(FTQ3Jieb1_8ChzqpB)oi7qwsogKoWEwmaVBn(YwO^C^@SLc* zz_v4-h?NI8_F9jMOm3ikWQsu1Q=<2tlkqv=ouPu_Vsb)1Aj$LM5)z2+0X#%OSv%Ym z#3~{3xXzZzoW#CX$;vT$KR+I(?gcUJw>IIAEj4PD@M( zx&F8}vdbZ$CQWZy+JX#Z(qjs z?G|KCE_UHcYU{~X3};ML%s|!4Gcp}?Nm@RF#XE;d3jJ=c=D0SxpLAr zrAJ3o(?=gvY)~>VH_xz_;>nZKr0puDH&u>7CM48~fpoDg0E>~Kq&Gz-c)P{fXYtVX zoGBUCKL=!2uvK_e&w(g~CH)_W()YdB3>!y>tTcR!9@kaD!DIyz_WZRPNOCa~K=dPvU(}kR zEYiiw6W*r6VP|+Hvb5+zCSwl!4;K7;%6Cg>(dv}`z`bUYr}s*jchZVI9Lm!K7uTm& zc%k5%8*Liau5`gIPS8ash{io$4wn@L&Qu*-krdY}Eua zQnZE3+c)H6tX6d}_wI7z$Rc#pU|60zOYy-=k>v1cjw|dsbjse82-&plYtNrLbR@u_ z#yKcpG1gM+gZZ;+Y9$wQIK=p1b9GE*qK=vt=He1Gg%^_f#&&)@L=vv`~q z`WRXW_tD%xEN@f^j0%u?;E_b>M8cFl-JX0CT2MRY(f6qoeLv7>0^-)+-=WjqNd}11 zX?b5gX#oRK0R4AG-`$g?;BZ2=I|je5$YAuWDt!*e5419GD}aIf!@}al3;+glMy3iM z3Q#C%Hx~&-=^z}~r;Awy7Ri$8Bx)soBRXfaZ03M#WGO;ewh*S!yy%7pSJ=gBe;FHK ztk9wCTcQvD{VXmh=80ze+Yn3lHK|Y&{R2Oed7m;1FbPf_W|cp2T&&hG;bDS*={?-P zs|UwXq**|c9Y;-YoYn2o!^bS-1JwvRiza|4!!oe{uyhxINzKCsDz75rzL9LAqd44z z=4#7eXb*pk)%it0lViS#CK}L3e({#7xb=Pg2rWfNs}1dIN97lD7}b+}%5-29IkY#H z8vWB&qx%hi)_f?))2A6p6~VJALd+|lv$*B0V?0>#+et{4UVJaBt5U68|cW_@g<_HT9C&rOJl5P0#rW+ ze~Vs=stK)>pd|YH0cZ43!tvV;-^1`jvMb5JVn0*2LM9!GjJ-R^fcHxr!QHyLTy8r{ z0TJwYW!$qye=8rbLoUP)h4^c#ZSTIAb{ig*v{Xd*O=UeTFXSv{fZBpp9P(aMBfS~F zqBAZh#g6etvRj>7V;8d)fCGubkJ1T&xroIVPd_j{biO*4Il2H7_1B@&9(Skx`yqZO z?qI-Cp4yqZHceM!=WE2S_ti6ObzXjsmdHK@u!*oFW6ZPWNe%@1Dd`FwJ z6d$gMZ==CAuc$C~*olxeuzASK;7h8&TAx;7fBRH2Uf_)9@#tBypUf0MvPn#MzhKMY z#za?0F!P;iZ97xZmTz)9^26%`xk#_l+6Py5KY23QUF6{xe74Yb=i2B_M`OE88g{&D z%xG+RfaEHdQ`ZWiLlgZ#h3C;epj#kiLDb{{u?Xt^O1FUj`eN`nEuei+pdB2$arwHEPbLqO(tF+nYrp#9^TTIj z*Q?&Oc}_2(n{nDysy)O+c@X4YTx!WfP`L#$J}A`q?$#Xwymh?%5;M{0M+99Hb)Z85 zCD=2djzJcF2)q!~-ze?@q5x?yh|xgh^$R2{HTN2Q7rGM8o;l+arPbExM=wEs_gTO` z->E{Kj2yzpbs-amLz^8RI#1T~-8{`)k~T5sm1e&x3PS_y4kydBQ>N+frmjWx4Wu1! z40;{Slmwc;N_&>^tl~`#LBnrmtc*Gep^Ghn#sVY^!rt!=RoC{HvOv`ZeHUh+zrujv z=_=BNglM}S^rnG|i1E_HuUV1T_BXsvQWp2>PmVT{b811W0Lj4=WCe&o^TFrTfKz9E zVaUxC*uoM!S{Ss_`IZ{DJae3qV&2tM^ZIY->*~Ub->R9J3|+A{tLKmHR<%aby6i1w zrq@o2=igH|i}ZPCr`2r49DIK4l?iJ0U!gP$`*oktE1d7C>4Ij4KHXNT{0pNq{CE2g+Abu@&G8#`Ewuz;e8$q5MJT5u;t8neTw6 zz}y^`w=teyx0<%C%#6ggo@ugLLEF?t?IG6Z$3gw=iQ|+ANp*(YwEnB(@#_?y_qIoj z*Q6wV^~iO(*-wV~$fK0X1_y--$APylV0W?^NeE2LU5CJ>Rb%-{a8kdT+(C_Fs$ z#%+1+&gYNcbMZ2a#5XG~e!u()N7=q0`H4!xz)=Q~&f|_J$_$tKO=*lo9B=rl(XzLS zrpXPS;*d*U-=}^4V>W`~rhoty0TnbsWn*EHqS@V~n1?z$I!Flp-syJ3*mMBLw_%}K z;RB57`yXDD zzrW8(|(W?!Qwy*U!R)kYxMEI`0iJxNoO;}1k zg4efOYx?~z$4+983PU`z&dm*~r+9int2j#sQ^LPWqTRc{Ws}{cD15z3T<7zqc?Y%q zY=i#bj>=UF$6JA~gN6q@@TaA}%(|oY>s-Q*3bUEtZ=C-+`^|iM{rA1^XJ0V$HYI$O zKW6MrsCV&v&+Qh(^5Am7SKVLGTld>$iY)79c=K7ic=HS#9nzIx#Vp9~^W?O(yGjuTUv)6NKvDZP=oDf!#D`2b z-#$bIokW`Zam)o4kB&@#HIjNe{Iw)(4V}uKm~5MGq7`{A5T{{iE{9Hp?4>B#QuXUV zwX3#}uxW**=(h_weHMDP8qauXr_zKTrOIz>&Dz{2+3a0xVEuIn8)$+ODzD5BAd_1p z!>j&(f7o`@!`pg+Uh23 z$aFK1fB18L&DH3k`*8)NEUYd6`xlew7w+NAptvZKJCufOZiSIcb3WqyQ56p zL?Ve=Zab7V4z+j7HtSxcR=1>{ipg~2ilcGu5_yh3*BP{>hcAbn$jO^!~|Y$5H6>>ydO5e>dsc$I6zi zQTpGjGmm9&b0&I|>t}z*H$^q`>(gY4G~zq;T`=dO$D_CHPo|6`;e6#s$8)8TzQHc+ zi+`S~5DrtB@+e6wD$t;eQ{LQ{Wa*Vcr@9cK&_t#<>qQz>X(IkxPke>FXX%PaG?VV( zdb1bNQJH(0eAg6m3ijW=8z~|a%_8F)8{PSpR0u^K&7u;MsEru>4xgd7s0W;2-vje92;R}8la<}}%+X^wqH<|;^ktY^Q%SIkkQ`Vfc&<5mr ze}vW(het0XBO_Hgj)A8Y6LxZ%Spr8|Wb+xeH z5}av}I$fATbT)`;_74ttX+uGelsVWEcGg@(Pg)vvEEFzd#w6YTlKv{=)Rg7JPbX{{ zVfd~d*b*-VBEnX?b2NPcKEe~N_)aGG$SQ8{hxn?*FaBCpBOAaMrQOs zd{te{KQQ<)Eg)uynrYZ@6Pwrkr`Z483smcm+}Wk)d3oFcvhTc3j+d4!QP2l){poMJ zTeqSO1&CNws8OwnB_$>OmT&U%+Pk{2sd!&O&1Yw4=jvLSkTCLNwi(LJ2`^sYcibr*B7?&}Z7+W*I|LlxS8;>yo0lS^3zq*sL*-#uhd-gqt?t|#K& z5gK}C$*29(f@y*hMKEF2DReEU$}WDKtOIoO;ll^$OI`xjGAu02FJ{>Ki~^N9w0A+n z=kN1euWwn2B+>Efy?+0}`&x$Lr?J=eOt;=smrAdbD+yYaj{#$0ROQfOT@L2q;p*C; z!Wjc4Jz{U}P&<(%RnZ5ZYutAV^oxG&>?m-_P&zm|F7~91^!3ePv043=((PBdMW|1X zC2rjKQQ|c_#FJl@x1 znf+$^=0xE@;nN=i3%?Yi4yBH2&P=f!`TJ4+JixShwY0L@+}@_e$AO>NRiM?8nVH#__+LLo#OhC3 z(foI9C}oN5ub`o)e}X*dTz+x!HLv3u#jKF9u&p0Iq$#XzY(RJ65!5UUjYh9tu}}{5 z0Fe%IckIH#85$bw+S>DJAjHm8-r}a6nwnzC)-o^{heiNl>mlyDckgb$diQP<+Vs@a z)DR(nt1D*zOppQe@&PBtFCd_ke_N+e{KADyzn*yL;L_s1&;<4w?4lWsGuJ@y zsd;!D0MZ{A9Asl->sKip_NB~)UWVfz%@p7hE~T8l`gVM+s;UZvtH62_5<>q1(f|Vt zUq?sB5*OMwBRiX2SXh{i&Cikf{P}cf&4UI}ctT)6Kt>9ra9@Up0oD)@5#i~UnyFY> zSpg`1^5jxG9;ibe5&F$SOY^P#+m~M=4I!X^5K8GGAt71%2Hj@}TWe-E58`;4;KKO$ zcqrdPNuj8y==4?S5e8nf1R9PmT)2QTGdBm(z$0jv_&cWE`g0Ptsn7-qeH7nPo2{S? z2!0Fn@!&9*a!N}b`d-P*&d%DA!kbQ?zUnYf+>Y04TYy?giWVjQV@NJt22z#NzW zGM1dWx?>m)xLRdMPn=Xu?n2+uz*eSF~r=+BekB>up&f?I}5U6%s zJv{srCwn5wJAbvb*Hn|leeAD9GqU+?DGyb|d2+X)wrB1*~4%}q~FZ@Dj% z8yFY}HUo#h63oj^98_AMZ5~^DQe_C^mxfRbdR93(KmDS>n3ntf%TBPBV9vYO%k`P4 zsO;z4se+nWiE+K*4+C^8`~)vqc|%$GD-_JYg9V*FN7w+;i7PGw{QRfqeYj#e!WFYX zIY0UkcCrDeaB#d%PGZ$B*q28iq&s8KSDt*o7v2ao#!VR+8Mw!PB8iGkQMc#|bg`?K zjNn1*Geqvl=g;Fh%#cGhTeX-bnxqOg02>OW5wK@8bS6cy%AubAo(+O_V0Tcl#6(Af z(kgG*dfM|~4dmIxzMyU+l%Y&>og<;dg}xF|L@N!tFfquc_xC|DJraqQ6uESXt^s-j zLDG93RAJEKs+uUo3aK`D7ZBUTu`ff}vT9 ze7F9ZfPi}8uqTWI$O64uckz9SAOlxu!+{5%<R9} z&dx^Bi&<*O$TWf-cXf3!n1Yvujt~o7U3y|1m`)=6!3F3tt(bBz1c#c{58o>YB;Z-v z*zof4O;kJcdPDmORI#AB#SWC)(5dGoG>&pviR*Av<2ZD1vT>rXX@U*y$utE)*$`dK~UHMENSN9;4p_`CPXIn z!Q6lTQC0?z6wUseUibn#JHEeDXK!yh)HtD4`}Rfc`}Z?TN(58G6|oZ2Pa(DtNax?S z9V_cF6L~4+7A*{NHi$NOvWCykNvy-c`L(qrmkZ?;^tGqD85#{dtgIQ(7QWP<1xo-$ zov`q5`T9A#v*GV_xhjF9p*JMB;uAGZy?}oQaa@X&{`g#!19-K@6nA9Qms)#!AA@-N nx2fI#Ef)SSLgfF)%O^xP 12-likelihood of coin toss + +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/3-bio-norm-dist.html b/docs/statistics/3-bio-norm-dist.html index b977ace..6a7b0f8 100644 --- a/docs/statistics/3-bio-norm-dist.html +++ b/docs/statistics/3-bio-norm-dist.html @@ -207,6 +207,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/4-pair-plots.html b/docs/statistics/4-pair-plots.html index aacb61d..b8865e4 100644 --- a/docs/statistics/4-pair-plots.html +++ b/docs/statistics/4-pair-plots.html @@ -210,6 +210,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/5-2d-decision-boundary.html b/docs/statistics/5-2d-decision-boundary.html index f7b5fe8..160e4f3 100644 --- a/docs/statistics/5-2d-decision-boundary.html +++ b/docs/statistics/5-2d-decision-boundary.html @@ -207,6 +207,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/6-mutli-variable-distributions.html b/docs/statistics/6-mutli-variable-distributions.html index c39dbbd..d1cd8b4 100644 --- a/docs/statistics/6-mutli-variable-distributions.html +++ b/docs/statistics/6-mutli-variable-distributions.html @@ -207,6 +207,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/7-mvnormal-sampling.html b/docs/statistics/7-mvnormal-sampling.html index 0e1ae59..58eaaed 100644 --- a/docs/statistics/7-mvnormal-sampling.html +++ b/docs/statistics/7-mvnormal-sampling.html @@ -207,6 +207,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/8-dotplot.html b/docs/statistics/8-dotplot.html index 7945cde..82764c6 100644 --- a/docs/statistics/8-dotplot.html +++ b/docs/statistics/8-dotplot.html @@ -207,6 +207,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/docs/statistics/9-2d-decision-boundary.html b/docs/statistics/9-2d-decision-boundary.html index af9d1d8..9cbf30b 100644 --- a/docs/statistics/9-2d-decision-boundary.html +++ b/docs/statistics/9-2d-decision-boundary.html @@ -210,6 +210,12 @@ 12-likelihood of coin toss
    +
    • +
  • +
    + + 13-understanding-bayes +
  • diff --git a/statistics/12-likelihood-coin-toss.qmd b/statistics/12-likelihood-coin-toss.qmd index b910969..969bf9b 100644 --- a/statistics/12-likelihood-coin-toss.qmd +++ b/statistics/12-likelihood-coin-toss.qmd @@ -80,3 +80,15 @@ for i in success_arr end fig ``` + +### 5. 最大似然率估计的正式方法 +`Distributions.jl`提供了`fit_mle`函数 +```{julia} +@doc(fit_mle) +``` + +```{julia} +for i in 1:10 +@info "total=10,success=$i mle"=> fit_mle(Binomial, 10,[i]) +end +``` \ No newline at end of file diff --git a/statistics/13-understand-bayes.qmd b/statistics/13-understand-bayes.qmd new file mode 100644 index 0000000..501fe89 --- /dev/null +++ b/statistics/13-understand-bayes.qmd @@ -0,0 +1,64 @@ +--- +title: "13-understanding-bayes" +author: "math4mads" +jupyter: julia-1.9 +--- + +:::{.callout-note title="简介"} + 下面我们考察一下一枚硬币真实参数为:`p(head)=0.6`的硬币在不同观察结果(`obs`)下的似然率估计 + + 估计两个点`p(head)=[a=0.5,b=0.8]` + + 在这两个点的最大似然率的估计可以表示为赔率(odds), 注意这不是概率 + + 可以直观的在概率密度图上观察到两者的差异 + + 当 观测值逐渐增加以后, 两个点的条件概率会趋近于 0, 赔率会趋近于 1 + $$L(a,b)=\frac{p(obs|a)}{p(obs|b)}$$ +::: + +## 1. load package +```{julia} + using Distributions,GLMakie +``` + +## 2. define likelihood function +```{julia} +xc=[0.5,0.8] +function plot_likelihood(;p=0.6,n=10,s=6) + probrange=range(0.0,1.0,50) + fig=Figure(resolution=(500,300)) + ax=Axis(fig[1,1]) + fun(p)=Binomial(n,p)|>d->pdf(d,s) + data=Float64[fun(p) for p in probrange] + + + pdf1=fun(xc[1]);pdf2=fun(xc[2]) + lines!(ax,probrange,data,label="n=$(n),success=$(s)") + scatter!(ax,xc,[pdf1,pdf2];color=[:red,:green],markersize=20) + lines!(ax,[xc[1],xc[2],xc[2]],[pdf1,pdf1,pdf2];linestyle=:dot) + axislegend(halign =:left, valign =:top) + fig +end +``` + +--- +## 3. plot results + +### 3.1 10 次硬币试验,6 次正面 +两个点的高度相除就是基于实验观察的赔率 +```{julia} +plot_likelihood() +``` + +### 3.2 100次试验 60 次正面 + +```{julia} +plot_likelihood(;p=0.6,n=100,s=60) +``` + +### 3.3 500次试验,300 次正面 + +```{julia} +plot_likelihood(;p=0.6,n=500,s=300) +``` \ No newline at end of file