Pymc3 model where the results of a switch are directly observed












2















I have just started learning pymc3 so I might be thinking about this completely the wrong way.



Assume that we observe a vector of 10 booleans.



The process of interest generates (observed) booleans with a Bernoulli distribution with a parameter theta1. So I define a Beta prior over theta1 and define a variable with length 10 that is a sample from Bernoulli(theta1).



However, this true sample is disturbed by sometimes switching the true data to 0, with a probability theta2. So I define a switch to 0 with a probability Bernoulli(theta2).



The switched values are the observed ones. I am not sure how to tell the model that I observed the switched variables, i.e. I am not sure how to fit the model to the observed data.



This is what I have for now, and I am kind of stuck:



# observed data (already switched)
observed_data = np.random.binomial(1, 0.5, size=10)

with pm.Model() as skeptic_model:
# uniform probability of the bernoulli parameter
true_model_prior = pm.Beta("true_model_prior", 1, 1)
true_data = pm.Bernoulli("true_data", p=true_model_prior, shape=data.shape)
disturbed_data = pm.math.switch(pm.Bernoulli("disturbed", 0.1), true_data, 0)









share|improve this question























  • Before using pymc3, you have to create that second array, the one with the switched values. How do you do it?

    – David
    Nov 13 '18 at 6:41
















2















I have just started learning pymc3 so I might be thinking about this completely the wrong way.



Assume that we observe a vector of 10 booleans.



The process of interest generates (observed) booleans with a Bernoulli distribution with a parameter theta1. So I define a Beta prior over theta1 and define a variable with length 10 that is a sample from Bernoulli(theta1).



However, this true sample is disturbed by sometimes switching the true data to 0, with a probability theta2. So I define a switch to 0 with a probability Bernoulli(theta2).



The switched values are the observed ones. I am not sure how to tell the model that I observed the switched variables, i.e. I am not sure how to fit the model to the observed data.



This is what I have for now, and I am kind of stuck:



# observed data (already switched)
observed_data = np.random.binomial(1, 0.5, size=10)

with pm.Model() as skeptic_model:
# uniform probability of the bernoulli parameter
true_model_prior = pm.Beta("true_model_prior", 1, 1)
true_data = pm.Bernoulli("true_data", p=true_model_prior, shape=data.shape)
disturbed_data = pm.math.switch(pm.Bernoulli("disturbed", 0.1), true_data, 0)









share|improve this question























  • Before using pymc3, you have to create that second array, the one with the switched values. How do you do it?

    – David
    Nov 13 '18 at 6:41














2












2








2


1






I have just started learning pymc3 so I might be thinking about this completely the wrong way.



Assume that we observe a vector of 10 booleans.



The process of interest generates (observed) booleans with a Bernoulli distribution with a parameter theta1. So I define a Beta prior over theta1 and define a variable with length 10 that is a sample from Bernoulli(theta1).



However, this true sample is disturbed by sometimes switching the true data to 0, with a probability theta2. So I define a switch to 0 with a probability Bernoulli(theta2).



The switched values are the observed ones. I am not sure how to tell the model that I observed the switched variables, i.e. I am not sure how to fit the model to the observed data.



This is what I have for now, and I am kind of stuck:



# observed data (already switched)
observed_data = np.random.binomial(1, 0.5, size=10)

with pm.Model() as skeptic_model:
# uniform probability of the bernoulli parameter
true_model_prior = pm.Beta("true_model_prior", 1, 1)
true_data = pm.Bernoulli("true_data", p=true_model_prior, shape=data.shape)
disturbed_data = pm.math.switch(pm.Bernoulli("disturbed", 0.1), true_data, 0)









share|improve this question














I have just started learning pymc3 so I might be thinking about this completely the wrong way.



Assume that we observe a vector of 10 booleans.



The process of interest generates (observed) booleans with a Bernoulli distribution with a parameter theta1. So I define a Beta prior over theta1 and define a variable with length 10 that is a sample from Bernoulli(theta1).



However, this true sample is disturbed by sometimes switching the true data to 0, with a probability theta2. So I define a switch to 0 with a probability Bernoulli(theta2).



The switched values are the observed ones. I am not sure how to tell the model that I observed the switched variables, i.e. I am not sure how to fit the model to the observed data.



This is what I have for now, and I am kind of stuck:



# observed data (already switched)
observed_data = np.random.binomial(1, 0.5, size=10)

with pm.Model() as skeptic_model:
# uniform probability of the bernoulli parameter
true_model_prior = pm.Beta("true_model_prior", 1, 1)
true_data = pm.Bernoulli("true_data", p=true_model_prior, shape=data.shape)
disturbed_data = pm.math.switch(pm.Bernoulli("disturbed", 0.1), true_data, 0)






python-3.x pymc3






share|improve this question













share|improve this question











share|improve this question




share|improve this question










asked Nov 10 '18 at 17:07









whatamesswhatamess

639




639













  • Before using pymc3, you have to create that second array, the one with the switched values. How do you do it?

    – David
    Nov 13 '18 at 6:41



















  • Before using pymc3, you have to create that second array, the one with the switched values. How do you do it?

    – David
    Nov 13 '18 at 6:41

















Before using pymc3, you have to create that second array, the one with the switched values. How do you do it?

– David
Nov 13 '18 at 6:41





Before using pymc3, you have to create that second array, the one with the switched values. How do you do it?

– David
Nov 13 '18 at 6:41












1 Answer
1






active

oldest

votes


















1














Your model can be reframed as a product of Bernoulli random variables, and therefore as a single Bernoulli random variable with a multiplicative p. Namely, the following model is equivalent to yours:



# observed data (already considered zero-inflated)
Y = np.random.binomial(1, 0.5, size=10)

with pm.Model() as zero_inflated_beta_bernoulli:
# true_model_prior
p = pm.Beta('p', alpha=1, beta=1)

# dropout rate
d = 0.1

# disturbed_data;
y = pm.Bernoulli('y', p = (1-d)*p, observed=Y)


You could let the dropout rate also be a random variable,



# dropout rate
d = pm.Beta('d', mu=0.1, sd=0.02)


However, it should be noted that this model really can't distinguish between dropouts and original outcomes, so the posteriors are sensitive to the priors.






share|improve this answer

























    Your Answer






    StackExchange.ifUsing("editor", function () {
    StackExchange.using("externalEditor", function () {
    StackExchange.using("snippets", function () {
    StackExchange.snippets.init();
    });
    });
    }, "code-snippets");

    StackExchange.ready(function() {
    var channelOptions = {
    tags: "".split(" "),
    id: "1"
    };
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function() {
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled) {
    StackExchange.using("snippets", function() {
    createEditor();
    });
    }
    else {
    createEditor();
    }
    });

    function createEditor() {
    StackExchange.prepareEditor({
    heartbeatType: 'answer',
    autoActivateHeartbeat: false,
    convertImagesToLinks: true,
    noModals: true,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: 10,
    bindNavPrevention: true,
    postfix: "",
    imageUploader: {
    brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
    contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
    allowUrls: true
    },
    onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    });


    }
    });














    draft saved

    draft discarded


















    StackExchange.ready(
    function () {
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53241342%2fpymc3-model-where-the-results-of-a-switch-are-directly-observed%23new-answer', 'question_page');
    }
    );

    Post as a guest















    Required, but never shown

























    1 Answer
    1






    active

    oldest

    votes








    1 Answer
    1






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes









    1














    Your model can be reframed as a product of Bernoulli random variables, and therefore as a single Bernoulli random variable with a multiplicative p. Namely, the following model is equivalent to yours:



    # observed data (already considered zero-inflated)
    Y = np.random.binomial(1, 0.5, size=10)

    with pm.Model() as zero_inflated_beta_bernoulli:
    # true_model_prior
    p = pm.Beta('p', alpha=1, beta=1)

    # dropout rate
    d = 0.1

    # disturbed_data;
    y = pm.Bernoulli('y', p = (1-d)*p, observed=Y)


    You could let the dropout rate also be a random variable,



    # dropout rate
    d = pm.Beta('d', mu=0.1, sd=0.02)


    However, it should be noted that this model really can't distinguish between dropouts and original outcomes, so the posteriors are sensitive to the priors.






    share|improve this answer






























      1














      Your model can be reframed as a product of Bernoulli random variables, and therefore as a single Bernoulli random variable with a multiplicative p. Namely, the following model is equivalent to yours:



      # observed data (already considered zero-inflated)
      Y = np.random.binomial(1, 0.5, size=10)

      with pm.Model() as zero_inflated_beta_bernoulli:
      # true_model_prior
      p = pm.Beta('p', alpha=1, beta=1)

      # dropout rate
      d = 0.1

      # disturbed_data;
      y = pm.Bernoulli('y', p = (1-d)*p, observed=Y)


      You could let the dropout rate also be a random variable,



      # dropout rate
      d = pm.Beta('d', mu=0.1, sd=0.02)


      However, it should be noted that this model really can't distinguish between dropouts and original outcomes, so the posteriors are sensitive to the priors.






      share|improve this answer




























        1












        1








        1







        Your model can be reframed as a product of Bernoulli random variables, and therefore as a single Bernoulli random variable with a multiplicative p. Namely, the following model is equivalent to yours:



        # observed data (already considered zero-inflated)
        Y = np.random.binomial(1, 0.5, size=10)

        with pm.Model() as zero_inflated_beta_bernoulli:
        # true_model_prior
        p = pm.Beta('p', alpha=1, beta=1)

        # dropout rate
        d = 0.1

        # disturbed_data;
        y = pm.Bernoulli('y', p = (1-d)*p, observed=Y)


        You could let the dropout rate also be a random variable,



        # dropout rate
        d = pm.Beta('d', mu=0.1, sd=0.02)


        However, it should be noted that this model really can't distinguish between dropouts and original outcomes, so the posteriors are sensitive to the priors.






        share|improve this answer















        Your model can be reframed as a product of Bernoulli random variables, and therefore as a single Bernoulli random variable with a multiplicative p. Namely, the following model is equivalent to yours:



        # observed data (already considered zero-inflated)
        Y = np.random.binomial(1, 0.5, size=10)

        with pm.Model() as zero_inflated_beta_bernoulli:
        # true_model_prior
        p = pm.Beta('p', alpha=1, beta=1)

        # dropout rate
        d = 0.1

        # disturbed_data;
        y = pm.Bernoulli('y', p = (1-d)*p, observed=Y)


        You could let the dropout rate also be a random variable,



        # dropout rate
        d = pm.Beta('d', mu=0.1, sd=0.02)


        However, it should be noted that this model really can't distinguish between dropouts and original outcomes, so the posteriors are sensitive to the priors.







        share|improve this answer














        share|improve this answer



        share|improve this answer








        edited Nov 14 '18 at 22:59

























        answered Nov 13 '18 at 18:31









        mervmerv

        25.1k673109




        25.1k673109






























            draft saved

            draft discarded




















































            Thanks for contributing an answer to Stack Overflow!


            • Please be sure to answer the question. Provide details and share your research!

            But avoid



            • Asking for help, clarification, or responding to other answers.

            • Making statements based on opinion; back them up with references or personal experience.


            To learn more, see our tips on writing great answers.




            draft saved


            draft discarded














            StackExchange.ready(
            function () {
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53241342%2fpymc3-model-where-the-results-of-a-switch-are-directly-observed%23new-answer', 'question_page');
            }
            );

            Post as a guest















            Required, but never shown





















































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown

































            Required, but never shown














            Required, but never shown












            Required, but never shown







            Required, but never shown







            Popular posts from this blog

            Full-time equivalent

            さくらももこ

            13 indicted, 8 arrested in Calif. drug cartel investigation