0

I am reading a statistics textbook Introduction to Statistics for Engineers by Sheldon Ross, p.275 and trying to re-do its examples on paper and in Octave. I am not able to replicate many Bayes calculations in Octave when it comes to the integration part. Please advise how to go about replicating below calculation in Octave? Below is a simple Bayes estimator example which naturally becomes a symbolic integration problem, where I often encounter difficulty doing in Octave.

[Clarification: This calculation is from a textbook and I understand it by hand. What I don't understand is how should one approach such statistical computing exercises in practice. This question relates to statistical/scientific computing, not coding or statistics per se.]

Suppose enter image description here are independent Bernoulli variables, having pdf enter image description here

p is the unknown variable enter image description here . Compute the Bayes estimator for p.

We know that enter image description here

The conditional pdf of p given X is then

enter image description here

It can be shown that, enter image description here ---(1)

Using (1) and letting enter image description here, the conditional pdf becomes

enter image description here

Recall Bayes estimator is enter image description here.

Therefore, Bayes estimator for p is:

enter image description here

Now, I try to replicate these steps using Octave as below and failed (integration took 40mins on my $2500 Dell desktop). Can you show my confused soul how do you do the above steps in Octave or Matlab or R to arrive at the same Bayes estimator?

#Use Octave to derive the above Bayes estimator
pkg load symbolic;
syms p n x;
f = (p^x) * (1-p)^(n-x);
F = int(f, p, [0, 1]); #integrate f, which gives the conditional pdf denominator
f_conditional = f/F; #the conditional pdf
integrand = p * f_conditional; # the integrand to derive Bayes estimator
estimator = int(integrand, p, [0, 1]);
#this integration takes forever, how else should I replicate the above in Octave?
limestreetlab
  • 173
  • 1
  • 11
  • This doesn't sound like a coding problem. You may want to have a look at [How to ask?](https://stackoverflow.com/help/how-to-ask) and provide a [minimal reproducible example](https://stackoverflow.com/help/minimal-reproducible-example) to help people with limited statistical knowledge or who are tired of formulas^^ – max Jul 12 '20 at 10:19
  • 1
    The question and presented calculations are confusing. What is the actual question? It sounds like you want to find the posterior probability for `p`, where the likelihood is the binomial distribution; this posterior is known to be beta distributed, for a suitable beta (possibly uninformative) prior. Then you can plug this into a Bayesian Estimator; in your case it seems you are assuming a quadratic loss; therefore the Bayes Estimate is simply the mean of the beta distribution, which has the simple analytical form `a/(a+b)`. So what is the question, and why isn't the answer simply a/(a+b)? – Tasos Papastylianou Jul 12 '20 at 12:46
  • Also, I'm usually the first to complain about pedantic "this isn't matlab"-style comments, but this _really_ isn't Matlab related at all; matlab's symbolic toolbox gives entirely different results to octave/sympy here. Furthermore, as max pointed out, this doesn't seem to be a coding problem either; if you're more interested in the theoretical part of the question (or the best practical way to implement such theory), then you might want to ask this over at [Cross-Validated](https://stats.stackexchange.com/) instead. As it stands there's no bug to help with here. – Tasos Papastylianou Jul 12 '20 at 12:49
  • @Tasos, the question isn't statistics per se. Those equations are straight from a textbook example, which I understand. But I want to do examples by hand for understanding and then in Octave for practice. I don't understand the quadratic loss etc (this a 2nd-year undergrad textbook). The question is 1. calculating a Bayes estimator involves conditional pdf, so we encounter naturally symbolic integration, right?, 2. given symbolic integration seems slow in Matlab/Octave, how should I approach doing pdf integration then? and 3. how would you replicate this example in Octave? – limestreetlab Jul 12 '20 at 13:23

0 Answers0