{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "b58ad65c-ebf7-43fe-bc22-8f2eb0db99fb", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to D:\\OneDrive\\public\\0043\\pvalue_and_posterior.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"D:\\\\OneDrive\\\\public\\\\0043\\\\pvalue_and_posterior.gif\")" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Distributions\n", "using StatsPlots\n", "default(fmt=:png)\n", "\n", "safediv(x, y) = x == 0 ? zero(x*y) : x/y\n", "\n", "function pvalue(k, n, p)\n", " z = safediv(k - n*p, āˆš(n*p*(1 - p)))\n", " 2ccdf(Normal(0, 1), abs(z))\n", "end\n", "\n", "posterior(k, n; a=1, b=1) = Beta(a+k, b+n-k)\n", "\n", "N = 200\n", "pā‚€ = 0.3\n", "sample = rand(Bernoulli(pā‚€), N)\n", "anim = @animate for n in [1:N; fill(N, 20)]\n", " k = sum(sample[1:n])\n", " P = plot(p -> pvalue(k, n, p), 0, 1; label=\"\")\n", " plot!(xtick=0:0.1:1, ytick=0:0.1:1)\n", " plot!(xguide=\"p\", yguide=\"P-value\")\n", " title!(\"P-value function (n=$n, k=$k)\")\n", " Q = plot(p -> pdf(posterior(k, n), p), 0, 1; label=\"\")\n", " plot!(xtick=0:0.1:1, ytick=0:20)\n", " plot!(xguide=\"p\", yguide=\"probability density\")\n", " title!(\"posterior pdf (n=$n, k=$k)\")\n", " plot(P, Q; size=(800, 300))\n", " plot!(titlefontsize=12)\n", " plot!(leftmargin=4Plots.mm, bottommargin=4Plots.mm)\n", "end\n", "gif(anim, \"pvalue_and_posterior.gif\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "b43cbd5d-8303-4315-ae86-828927aa1a60", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to D:\\OneDrive\\public\\0043\\pvalue_vs_bayesian.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"D:\\\\OneDrive\\\\public\\\\0043\\\\pvalue_vs_bayesian.gif\")" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Distributions\n", "using StatsPlots\n", "default(fmt=:png)\n", "\n", "safediv(x, y) = x == 0 ? zero(x*y) : x/y\n", "\n", "function pvalue(k, n, p)\n", " z = safediv(k - n*p, āˆš(n*p*(1 - p)))\n", " 2ccdf(Normal(0, 1), abs(z))\n", "end\n", "\n", "posterior(k, n; a=1, b=1) = Beta(a+k, b+n-k)\n", "\n", "function pvalue_bayes(k, n, p; a=1/3, b=1/3)\n", " post = posterior(k, n; a, b)\n", " min(1, 2cdf(post, p), 2ccdf(post, p))\n", "end\n", "\n", "N = 200\n", "pā‚€ = 0.3\n", "sample = rand(Bernoulli(pā‚€), N)\n", "anim = @animate for n in [1:N; fill(N, 20)]\n", " k = sum(sample[1:n])\n", " plot(p -> pvalue(k, n, p), 0, 1; label=\"ordibary P-value\")\n", " plot!(p -> pvalue_bayes(k, n, p), 0, 1; label=\"Baysian P-value\", ls=:dash)\n", " plot!(xtick=0:0.1:1, ytick=0:0.1:1)\n", " plot!(xguide=\"p\", yguide=\"P-value\")\n", " title!(\"n=$n, k=$k\")\n", "endc\n", "gif(anim, \"pvalue_vs_bayesian.gif\")" ] }, { "cell_type": "code", "execution_count": null, "id": "c54c3899-c62a-4fd7-8923-daf9945233eb", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.9.1", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.1" } }, "nbformat": 4, "nbformat_minor": 5 }