Created
December 12, 2019 09:27
-
-
Save aakhmetz/a2a754a2d6f31209010cc7bd8eeeb64a to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "┌ Warning: Environment variable CMDSTAN_HOME not set. Use set_cmdstan_home!.\n", | |
| "└ @ CmdStan /home/aakhmetz/.julia/packages/CmdStan/Ohkw8/src/CmdStan.jl:27\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "using CmdStan" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "\"/tmp/jl_bjO0OI\"" | |
| ] | |
| }, | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "set_cmdstan_home!(homedir() * \"/cmdstan-git/\")\n", | |
| "prjdir = dirname(@__DIR__)\n", | |
| "tmpdir = mktempdir()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "(1000, 1000, 4)" | |
| ] | |
| }, | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "nwarmup, nsamples, nchains = 1000, 1000, 4" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "\"\\ndata {\\n int<lower = 1> J; // number of cases\\n real<lower = 0> onset_day_reported[J]; //day of illness onset\\n}\\n\\nparameters {\\n real<lower = 0> incubation[J];\\n\\n real mu;\\n real<lower = 0> sigma;\\n}\\n\\nmodel {\\n mu ~ normal(0, 5);\\n sigma ~ cauchy(0, 5);\\n\\n incubation ~ lognormal(mu,sigma);\\n\\n for (j in 1:J) {\\n target += normal_lpdf(incubation[j] | onset_day_reported[j]+0.5, 0.5);\\n }\\n\\n\\n}\\n\\ngenerated quantities {\\n real incubation_mean = exp(mu+sigma^2/2.0);\\n real incubation_sd = sqrt((exp(sigma^2)-1.0)*exp(2.0*mu+sigma^2));\\n}\\n\"" | |
| ] | |
| }, | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "stan_code = \"\n", | |
| "data {\n", | |
| " int<lower = 1> J; // number of cases\n", | |
| " real<lower = 0> onset_day_reported[J]; //day of illness onset\n", | |
| "}\n", | |
| "\n", | |
| "parameters {\n", | |
| " real<lower = 0> incubation[J];\n", | |
| "\n", | |
| " real mu;\n", | |
| " real<lower = 0> sigma;\n", | |
| "}\n", | |
| "\n", | |
| "model {\n", | |
| " mu ~ normal(0, 5);\n", | |
| " sigma ~ cauchy(0, 5);\n", | |
| "\n", | |
| " incubation ~ lognormal(mu,sigma);\n", | |
| "\n", | |
| " for (j in 1:J) {\n", | |
| " target += normal_lpdf(incubation[j] | onset_day_reported[j]+0.5, 0.5);\n", | |
| " }\n", | |
| "\n", | |
| "\n", | |
| "}\n", | |
| "\n", | |
| "generated quantities {\n", | |
| " real incubation_mean = exp(mu+sigma^2/2.0);\n", | |
| " real incubation_sd = sqrt((exp(sigma^2)-1.0)*exp(2.0*mu+sigma^2));\n", | |
| "}\n", | |
| "\"" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "df_inc = DataFrame(Dict(\"n\" => [1, 5, 7, 9, 17, 18, 16, 6, 8, 10, 7, 1, 2, 4, 4, 1],\n", | |
| " \"day\" => [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 25]))\n", | |
| "days_inc = vcat(fill.(df_inc[!,:day],df_inc[!,:n])...);" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Dict{String,Any} with 2 entries:\n", | |
| " \"J\" => 116\n", | |
| " \"onset_day_reported\" => [6, 7, 7, 7, 7, 7, 8, 8, 8, 8 … 18, 19, 19, 19, 19,…" | |
| ] | |
| }, | |
| "execution_count": 9, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "days_inc = vcat(fill.(df_inc[!,:day],df_inc[!,:n])...)\n", | |
| "stan_data = Dict(\"J\" => length(days_inc), \"onset_day_reported\" => days_inc)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Dict{String,Float64} with 116 entries:\n", | |
| " \"incubation[111]\" => 19.5\n", | |
| " \"incubation[86]\" => 14.5\n", | |
| " \"incubation[87]\" => 14.5\n", | |
| " \"incubation[75]\" => 13.5\n", | |
| " \"incubation[101]\" => 16.5\n", | |
| " \"incubation[57]\" => 11.5\n", | |
| " \"incubation[97]\" => 15.5\n", | |
| " \"incubation[54]\" => 11.5\n", | |
| " \"incubation[77]\" => 13.5\n", | |
| " \"incubation[52]\" => 11.5\n", | |
| " \"incubation[49]\" => 11.5\n", | |
| " \"incubation[12]\" => 8.5\n", | |
| " \"incubation[7]\" => 8.5\n", | |
| " \"incubation[90]\" => 15.5\n", | |
| " \"incubation[35]\" => 10.5\n", | |
| " \"incubation[79]\" => 13.5\n", | |
| " \"incubation[13]\" => 8.5\n", | |
| " \"incubation[74]\" => 13.5\n", | |
| " \"incubation[16]\" => 9.5\n", | |
| " \"incubation[76]\" => 13.5\n", | |
| " \"incubation[27]\" => 10.5\n", | |
| " \"incubation[33]\" => 10.5\n", | |
| " \"incubation[78]\" => 13.5\n", | |
| " \"incubation[10]\" => 8.5\n", | |
| " \"incubation[103]\" => 16.5\n", | |
| " ⋮ => ⋮" | |
| ] | |
| }, | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "# stan_init = Dict(\"mu\" => 2.5, \"sigma\" => 0.2, \"incubation\" => days_inc .+ 0.5)\n", | |
| "stan_init = Dict(\"incubation[$i]\"=>Float64(v+.5) for (i,v) in enumerate(days_inc))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "\n", | |
| "File /tmp/jl_bjO0OI/noname.stan will be updated.\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "ename": "MethodError", | |
| "evalue": "MethodError: Cannot `convert` an object of type Dict{String,Float64} to an object of type Array{Dict{String,Any},1}\nClosest candidates are:\n convert(::Type{Array{T,N}}, !Matched::FillArrays.Zeros{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:320\n convert(::Type{Array{T,N}}, !Matched::FillArrays.Ones{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:320\n convert(::Type{Array{T,N}}, !Matched::FillArrays.AbstractFill{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:312\n ...", | |
| "output_type": "error", | |
| "traceback": [ | |
| "MethodError: Cannot `convert` an object of type Dict{String,Float64} to an object of type Array{Dict{String,Any},1}\nClosest candidates are:\n convert(::Type{Array{T,N}}, !Matched::FillArrays.Zeros{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:320\n convert(::Type{Array{T,N}}, !Matched::FillArrays.Ones{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:320\n convert(::Type{Array{T,N}}, !Matched::FillArrays.AbstractFill{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:312\n ...", | |
| "", | |
| "Stacktrace:", | |
| " [1] Stanmodel(::String, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::String, ::String, ::Array{String,1}, ::Array{Dict{String,Any},1}, ::String, ::Array{Cmd,1}, ::Sample, ::CmdStan.Random, ::Dict{String,Float64}, ::String, ::CmdStan.Output, ::Bool, ::String, ::String, ::Symbol) at /home/aakhmetz/.julia/packages/CmdStan/Ohkw8/src/main/stanmodel.jl:33", | |
| " [2] #Stanmodel#3(::String, ::Int64, ::Int64, ::Int64, ::Int64, ::String, ::Array{String,1}, ::Array{Dict{String,Any},1}, ::CmdStan.Random, ::Dict{String,Float64}, ::CmdStan.Output, ::Bool, ::String, ::String, ::Symbol, ::Type{Stanmodel}, ::Sample) at /home/aakhmetz/.julia/packages/CmdStan/Ohkw8/src/main/stanmodel.jl:184", | |
| " [3] (::Core.var\"#kw#Type\")(::NamedTuple{(:model, :init, :nchains, :num_warmup, :num_samples, :tmpdir),Tuple{String,Dict{String,Float64},Int64,Int64,Int64,String}}, ::Type{Stanmodel}, ::Sample) at ./none:0 (repeats 2 times)", | |
| " [4] top-level scope at In[10]:1" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "stan_model = Stanmodel(\n", | |
| " model = stan_code, \n", | |
| " init = stan_init,\n", | |
| " nchains = nchains,\n", | |
| " num_warmup = nwarmup,\n", | |
| " num_samples = nsamples,\n", | |
| " tmpdir=tmpdir);\n", | |
| "\n", | |
| "_, stan_chns, _ = stan(stan_model, stan_data, summary = false);" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Julia 1.3.0", | |
| "language": "julia", | |
| "name": "julia-1.3" | |
| }, | |
| "language_info": { | |
| "file_extension": ".jl", | |
| "mimetype": "application/julia", | |
| "name": "julia", | |
| "version": "1.3.0" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment