Make MCS @defmcs accept constrained distributions?

Using the MCS functionality:
I am trying to assign a distribution to a parameter alpha = Parameter(index=[1:5, 1:5]). I want to enforce that sum(alpha[:,:]) = 1 for all runs. I can see two ways to do this:

  1. Assign e.g. a normal distribution to all elements of alpha, and divide by the sum of the values taken: alpha[:,:] = Normal(1,0) / sum(alpha[:,:]). However this is not supported by @defmcs.
  2. Assign a Dirichlet distribution to alpha, which gives an array of values of sum=1: alpha[:,:] = Dirichlet(25,1). Yet this gives me a 25-element Array{Float64,1}, while I would like a 5*5 Array{Float64,2}. A way to get to the latter shape would be to use reshape(Dirichlet(25,1), 5,5). However this is not supported by @defmcs either.

Do you have suggestions as to how to solve this?

In FUND we have a similar situation, and we just manually solve it in the code, i.e. our parameters won’t add up to 1, but then we have a variable where things add up to 1, see here.

@plevin and @lrennels might have an idea that is less manual :slight_smile:

I’ve prototyped an approach using a pseudo distribution that does reshaping:

struct ReshapedDistribution
    dims::Vector{Int}
    dist::Distribution
end

function Base.rand(rd::ReshapedDistribution, draws::Int=1)
    values = rand(rd.dist, draws)
    dims = (draws == 1 ? rd.dims : [rd.dims..., draws])
    return reshape(values, dims...)
end

For example, rd = ReshapedDistribution([5, 5], Dirichlet(25,1)) behaves the way I think you want – each call to rand(rd) produces a 5x5 matrix that sums to 1.

The next step is to modify @defmcs to assign these draws to the indicated variable.

Before I continue with this, does it seem workable to you?

This looks great and would work for what I’m trying to do, I think. Thanks!

I’m wondering how https://juliastats.github.io/Distributions.jl/stable/matrix.html fits into the story here? Maybe not?

I don’t think it does? The only 2 distributions offered do not fit sum()=1.

Helene, I’ve completed the changes to support your request. Although I’ve created and tested a simple example, I haven’t tested it on a real model. Would you be willing to check out the mcs-reshape-dist branch of Mimi and try it with your model? This branch contains minor changes on top of the current master, which may be ahead of the version you’re using.

There is an example of usage in test/mcs/test_reshaping.jl and in the copy of the 2-region model at test/mcs/test-model. I added a 2-dimensional variable tester that is unused in the model, but given values in the MCS. You’d want to modify your model in a similar way (but actually use the data!)

@hbenveniste I think @plevin can make these changes for you and test your model, he’ll just need access to https://github.com/PrincetonUniversity/mig-fund-quintiles (assuming that’s where you’re doing your work)

Great! Just invited @plevin to contribute. Thanks a lot !

Helene, how do I run your MCS case? I don’t see any “main” code. Do you just run the top-level functions interactively?

Not exactly sure I’m understanding correctly, but yes you just include the functions mentioned at the beginning of the file, then call scc_mcs_nicefund for the number of iterations you want.

I’ve made the changes in your repo in a new branch called mcs-reshape-test.
I added a test directory to your repo, with runtests.jl and test_mcs.jl. You can run the simple test by including “test/runtests.jl”.

I used the Mimi branch mcs-reshape-dist, which hasn’t been merged into Mimi’s master yet. I think we’re almost ready to do that–I just ran a test of your MCS with N = 100, which completed, but saving the data failed, owning to the change in alpha to storing a Matrix of values for each trial. We’ll need to discuss how to handle that.

Once we merge these changes, you’ll have to change one line of code (see MCsim_nicefund.jl line 53) owing to the change in API from set_models!() to set_model!().

So it’s almost ready, and it might already be usable except for saving the trial data using the built-in method.

Thanks a lot! So this is probably a stupid question, but how do I use the Mimi branch mcs-reshape-dist? I don’t have a github repo for Mimi, only for Fund. In my understanding, as long as this branch hasn’t merged into the Mimi master branch, I can’t access it right?

There are no stupid questions! In Pkg mode, entered with ], you may change the branch in your current environment by typing add Mimi#branchname. When you are done you can just type free Mimi to bring it back to the most recently tagged version. Also make sure to ]up once in a while, we just released a new version of Mimi.

Thanks, but it doesn’t seem to be working. Whether I try to run the test files that Richard added, or directly my MCsim_nicefund.jl file, I keep having as error message
ERROR: UndefVarError: ReshapedDistribution not defined. Not sure what I’m doing wrong?

Hi, Helene. The changes I made to your code rely on pending changes to Mimi that exist only in the branch “mcs-reshape-dist”, so you would need to switch to that branch before testing my changes.

I hadn’t planned on your merging my changes into your master branch before testing was completed. I’m working on it now, though and hope to have it running for you by tomorrow. Then we can merge “mcs-reshape-dist” into Mimi’s master and you can merge the rest of my changes into your master branch and all should be consistent.

Sorry for the confusion!

Helene, was the function scc_mcs_nicefund working before I started making changes?

In particular, the computed variable scc is a Float64, which the CSV writer doesn’t like. (I believe it expects an iterable of iterables.)

Hi Richard. The function scc_mcs_nicefund is working apart from the point you raised on CSV.write.

I am confused about this branches thing. As Lisa suggested, I did switch to the mcs-reshape-dist branch, but I still get the error message on ReshapedDistribution not being defined.

Thanks!

The line that fails is:
CSV.write(joinpath(output_dir, "SCC.csv"), scc)

It fails because scc is a scalar float, where CSV.write expects to iterate over rows of values. Were you expecting this to be other than a Float64?

I’m going to check in what I’ve done, on the master (not mcs-rehape-test) branch of your repo and speak with Lisa about possibly following up with you. (Since you already attempted to merge mcs-reshape-test into master, I fixed it there. There were some conflicts to resolve…)

I’m happy to help out from here. We’re going to merge this into master tonight so you can just do add Mimi#master and just use the master branch of Mimi, and the master branch of your repository after you pull Rich’s changes.

Note that if you hav to read in data that was saved as a matrix, i.e. a reshaped parameter alpha, you’ll need to deserialize the csv file with deserialize("results/mcs/alpha.dat"), but we can approach that when you hit it.