{
"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
}