Merge pull request #6 from mcreel/devel

Devel

Project.toml

@@ -1,33 +1,29 @@
 name = "SimulatedNeuralMoments"
 uuid = "43b7c8c0-5622-45b8-90a3-338fc50d2232"
 authors = ["Michael Creel <[email protected]> and contributors"]
-version = "1.0.7"
+version = "2.0.0"
 
 [deps]
-AdvancedMH = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170"
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
-PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
-AdvancedMH = "0.6, 0.7"
 BSON = "0.3"
+Distributions = "0.25"
 Flux = "0.13"
 MCMCChains = "6.0"
-PrettyTables = "2.0"
-StatsPlots = "0.15"
 StatsBase = "0.33"
-Turing = "0.22, 0.23, 0.24"
+StatsPlots = "0.15"
 julia = "1.5 - 1.9"
 
 [extras]

README.md

@@ -12,7 +12,6 @@ With a trained net, we can feed in statistics computed from real data, and obtai
 
 ABC/MSM draws from the posterior are done using Turing.jl, concretely MCMC using Metropolis-Hastings sampling. The likelihood is an estimate of the asymptotic Gaussian distribution of the transformed statistics, which is justified due to the central limit theorem applying to the statistics. The proposal density is also the same Gaussian density, estimated at the transformed statistics corresponding to the real data. Because the proposal and the likelihood are identical, at the true parameter values, random walk MH sampling is effective, because MCMC concentrates around the true parameter values, as the sample size becomes large.
 
-Please see the docs, in the link at the top, and the READMEs for the two examples to see how this all works. To run the two examples, using small training and testing samples, do ``using Pkg; Pkg.test("SimulatedNeuralMoments")``
-* [Stochastic volatility model](https://github.com/mcreel/SimulatedNeuralMoments.jl/blob/main/examples/SV/README.md)
-* [Mixture of normals model](https://github.com/mcreel/SimulatedNeuralMoments.jl/blob/main/examples/MN/README.md)
+Please see the docs, in the link at the top, and the README for the two examples to see how this all works.
+* [example models](https://github.com/mcreel/SimulatedNeuralMoments.jl/blob/main/examples/README.md)
 

docs/API.md

@@ -8,37 +8,40 @@ three functions.
 # the type that holds the model specifics
 struct SNMmodel
     modelname::String # name of model
+    samplesize::Int64 # number of observations the model generates 
     lb # vector of lower bounds. All parameters must have a finite lower bound.
     ub # vector of upper bounds. All parameters must have a finite upper bound.
     insupport::Function # function that checks if a draw is valid (in bounds)
+    prior:: Function # the prior, used to do MCMC
     priordraw::Function # function that returns a draw from the prior
     auxstat::Function # function that returns an array of draws of statistic, or the statistic corresponding to the real data
 end
 ```
 One uses this to define the statistical model, the parameters of which we wish to
 estimate. This is done as
 ```
-model = SNMmodel(modelname, lb, ub, insupport, priordraw, auxstat)
+model = SNMmodel(modelname, samplesize, lb, ub, insupport, prior, priordraw, auxstat)
 ```
 
 ### The elements of a SNMmodel type must follow:
-* modelname is a string variable, which describes the model. 
+* modelname is a string variable, which describes the model.  
+* sample size is and integer, the number of data points the model generates
 * lb and ub are vectors of Float64. These are the lower and upper bounds of the support of the parameters (the parameter space). At present, only bounded parameter spaces are supported. 
 * insupport, prior and priordraw and auxstat are all strings, names of functions.
   * insupport(theta) returns a boolean, indicating if the theta is in the support.
+  * prior(theta) computes the prior density at theta
   * priordraw() returns a vector of Float64, a draw from the prior.
-  * auxstat is called as
+* auxstat is called as
     * ```auxstat(theta, reps)```, where reps is Int64. It returns reps draws of the statistics drawn at theta, in a reps dimensional vector of vectors.
     * or as ```auxstat(data)``` where data is an array. This returns the statistics computed at the real data.
 
-The files MNexample.jl and MNlib.jl in the examples/MN directory show how to use this type.
+The file example.jl the examples directory show how to use this type.
 
 ## The function MakeNeuralMoments
-```MakeNeuralMoments(model::SNMmodel, transf; TrainTestSize=1, Epochs=1000)```
+```MakeNeuralMoments(model::SNMmodel; TrainTestSize=1, Epochs=1000)```
  The function generates simulated data from the model and trains a neural net to recognize transformed parameters, given the statistic. It returns the trained net, and an array of information for transforming a raw statistic vector from auxstat to use it as an input to the net.
 ### This function has two required arguments:
 * model, which is of type SNMmodel, as discussed above. 
-* tranfs: This function transforms draws from the prior to lie in ℛⁿ. It is expected that this tranformation will be computed using Bijectors.jl, as in ```transf = bijector(@Prior)```, for example. See the MN or SV examples.
 ### optional arguments
 * TrainTestSize is an optional argument. If omitted, then 20000 samples per parameter
   will be used to train and test the net. The default gives fairly quick training for
@@ -53,24 +56,53 @@ The files MNexample.jl and MNlib.jl in the examples/MN directory show how to use
 * nninfo: the infomation needed to transform a raw vector of statistics, from auxstat, to be used as an input to the neural net.
 
 ## The function NeuralMoments
-```NeuralMoments(z, nnmodel, nninfo```
-This function takes a raw statistic, the trained neural net, and the array of transforming information, and returns a predicted transformed parameter vector.
+```NeuralMoments(z, model, nnmodel, nninfo```
+This function takes a raw statistic, the SNM model structure, the trained neural net, and the array of transforming information, and returns a predicted transformed parameter vector.
 ### The inputs are
 * z: a vector of raw statistics from the model, the output of auxstat(y), for example.
+* model: the SNMmodel structure that describes the statistical model being estimated
 * nnmodel: the trained neural net, a Flux.jl chain. It is the first of the outputs of MakeNeuralMoments.
 * nninfo: a tuple of vectors, used to transform raw statistics, z, to then pass as inputs to nnmodel. nninfo is the second output of MakeNeuralMoments.
 ### Output
-a vector which is the predicted value, from the neural net,  of the transformed parameters that generated the inputted z.
+a vector which is the predicted value, from the neural net,  of the transformed parameters that generated the argument z. 
 
 ## The function mΣ
 ```mΣ(θ, reps, model::SNMmodel, nnmodel, nninfo)```
 This function computes the mean and covariance matrix of ```reps``` evaluations of  ```NeuralMoments```, where the inputs, ```z```, to ```NeuralMoments``` are drawn at the trial parameter value ```θ```.
 ### The inputs are
-* ```θ```: a trial parameter value (untransformed)
+* ```θ```: a trial parameter value
 * reps: Int64. The number of simulation draws to use to compute the mean vector and the covariance matrix
 * model: an SNMmodel structure that describes the statistical model being estimated
 * nnmodel: the trained neural net, a Flux.jl chain. It is the first of the outputs of MakeNeuralMoments.
 * nninfo: a tuple of vectors, used to transform raw statistics, z, to then pass as inputs to nnmodel. nninfo is the second output of MakeNeuralMoments.
 ### Outputs
 * mbar: the mean vector of the reps draws from the NeuralMoments function
-* Σ: the covariance matrix of the reps draws from the NeuralMoments function 
+* Σ: the covariance matrix of the reps draws from the NeuralMoments function
+
+## the function snmobj
+```snmobj(θ, m, reps, model::SNMmodel, nnmodel, nninfo)```
+This function computes the MSM-CUE log likelihood at the parameter value θ.
+### The inputs are
+* ```θ```: a trial parameter value
+* m: the output of NeuralMoments, when the input z is computed using the real (not simulated) data
+* reps: an integer number of replications of the simulated statistics. 
+* model: an SNMmodel structure that describes the statistical model being estimated
+* nnmodel: the trained neural net, a Flux.jl chain. It is the first of the outputs of MakeNeuralMoments.
+* nninfo: a tuple of vectors, used to transform raw statistics, z, to then pass as inputs to nnmodel. nninfo is the second output of MakeNeuralMoments.
+### Outputs
+the log likelihood, a Float64
+
+## the function mcmc
+```mcmc(θ, length, lnL, model::SNMmodel, nnmodel, nninfo, proposal, burnin=100, verbosity=10)```
+### the inputs are
+* θ: the parameter vector used to initialize the chain
+* length: the desired length of chain, not counting burnin
+* lnL: the MSM-CUE log likelihood, as a closure (see the example)
+* model: the SNMmodel type that describes the statistical model being estimated
+* nnmodel: the trained neural net, a Flux.jl chain. It is the first of the outputs of MakeNeuralMoments.
+* nninfo: a tuple of vectors, used to transform raw statistics, z, to then pass as inputs to nnmodel. nninfo is the second output of MakeNeuralMoments.
+* proposal: a function that returns a proposed parameter vector
+* burnin: the number of burnin draws to be discarded before saving draws to the chain
+* verbosity: an integer: report intermediate results every so often
+### Outputs
+* chain: a reps X k+2 array. The first k columns hold the k parameters for each draw, the penultimate column holds the log-likelihood value, and the last column holds an indicator as to whether the value is a new accepted value (1) or a copy of the previous value (0). This is simple to facilitate computing the acceptance rate

examples/.gitignore

@@ -0,0 +1 @@
+Manifest.toml