simple inference problem given prior distributions and sum of variables observed

[Nikolay Evdokimov]

Hello. I am new to Stan and I try to solve simple inference problem.

I would like to solve a simple one before applying Stan to making inference in more complex models.

About the problem. The task is to estimate posterior destibution for variables given some observed values.

Given:

a ~ normal(mu=1, std =1);

b ~ normal(mu=2, std = 1.5);

In our single observation we get information that a+b = sumValue.  Let’s fix it at 4.

What is pdf for a variable now? I mean pdf(a | a+b==sumValue).

I have solved this problem with Matlab and Infer.Net library, but I can’t reproduce it with Stan.

I have read the manual, and especially  part 20. on Problematic Posteriors,  on  and some threads ( https://groups.google.com/forum/#!topic/stan-users/rKxrs2gsPr4.) on Forum. But I still can’t understand if Stan is suitable for this inference task.

The problem can be easily solved analytically:

pdfOf(a | a+b== sumValue)By(x) = pdfA(x).*pdfB(2-sumValue-x) / k, where k = Integral(pdfA(x).*pdfB(sumValue-x)).

See matlab code for this:

x = [-2:0.01:8];

pdfA = @(x)normpdf(x, 1, 1);

pdfB = @(x)normpdf(x, 2, 1.5);

sumValue = 4;

 

pdfA_condNoAdj = @(x)pdfA(x).*pdfB(sumValue-x);

pdfA_condTotal =  integral(pdfA_condNoAdj ,-10,10);

pdfA_cond = @(x)pdfA_condNoAdj(x)/pdfA_condTotal;

 

figure(1);clf

hold on;

plot(x, pdfA(x));

plot(x, pdfB(x));

plot(x, pdfA_cond(x));

hold off;

grid on;

legend('pdfA','pdfB','pdfA cond');

 

So we see that posterior distribution for A is easily calculated.

Now let’s try to that in infer.net

        static void Main(string[] args)

        {

            Variable<double> a = Variable.GaussianFromMeanAndVariance(1, 1).Named("a");

            Variable<double> b = Variable.GaussianFromMeanAndVariance(2, 1.5).Named("b");

            Variable<double> c;

 

            c = a + b;

            c.ObservedValue = 4;

 

            InferenceEngine engine = new InferenceEngine();

            engine.Algorithm = new ExpectationPropagation();

            Console.WriteLine("Dist over x=" + engine.Infer(a));

            Console.ReadKey();

        }

The engine makes inference that is close to Matlab’s version, but not identical. Posterior PdfA  is normal with mu=1.4, std = 0.6.

How can I get posterior  distribution for a in Stan?

I try to do something like following:

data {   real cInput; // s.e. of effect estimates   }

parameters {   real a;    real b; }

model {

    a ~ normal(0,1);

    b ~ normal(0,1);

    cInput<-a + b;

}

And I would NOT like use model where C = a+b + noise, because it would change the problem.

Can Stan solve the problem above the way Matlab and Infer.Net did?

How can I infer posterior probability density function? 

Is there some way to assign observed value to linear combination (or any function) of parameters?

Regards, 

Nikolay

[Bob Carpenter]

Your model's only identified by the prior, but that's no big deal
in this case. There's also only a single degree of freedom, because
given a, b = 4 - a. From your note, you're assuming the density

p(a) propto normal(a | 1, 1) * normal(1 - a | 2, 1.5)

In Stan, you can do that directly and define b along the way or not:

parameters {
real a;
}
transformed parameters {
real b;
b <- 4 - a;
}
model {
a ~ normal(1, 1);
b ~ normal(2, 1.5);
}

You could've also skipped the transformed parameters block and
written the second line of the model as:

(4 - a) ~ normal(2, 1.5);

Either way, given your intent, you can ignore the Jacobian warning.

This behaves well and produces a posterior sample with mean
and quantiles:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a 1.30 0.02 0.83 -0.30 0.73 1.27 1.84 2.93 1470 1
b 2.70 0.02 0.83 1.07 2.16 2.73 3.27 4.30 1470 1
lp__ -0.65 0.02 0.71 -2.59 -0.82 -0.38 -0.20 -0.15 1980 1

That doesn't match your report from Infer.NET, but I have no idea
what Infer.NET does to fit.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Nikolay Evdokimov]

Thanks a lot. That worked and helped. I tried the same for recovering random walk with noisy observations and that worked.

data {

  vector[10] yObserved ; // s.e. of effect estimates 

}

parameters {

  vector[10] dx; 

}

transformed parameters

{

   vector[10] yTrue; 

   vector[10] yObservationErrors;

   yTrue <-  cumulative_sum(dx); 

   yObservationErrors <- yObserved  -yTrue;

}

суббота, 2 января 2016 г., 19:38:24 UTC+3 пользователь Nikolay Evdokimov написал:

[Bob Carpenter]

I don't know what you mean by random walk here.

The Stan program needs a model block to constrain
dx in some way or the posterior will be improper.

- Bob

[Nikolay Evdokimov]

Ah, sure. I didn't copy paste that part from model:

model {

    yObservationErrors ~ normal(0,1);

    dx ~ normal(0,1);

}

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/Qhw0A2gr67g/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

> On Jan 2, 2016, at 5:21 PM, Nikolay Evdokimov <kolyaev...@gmail.com> wrote:
>
> Ah, sure. I didn't copy paste that part from model:
>
> model {
> yObservationErrors ~ normal(0,1);
> dx ~ normal(0,1);
> }

That's an unorthodox way of writing the model, but it should
be OK --- Stan only cares about the density that's defined.

- Bob