Predictive distributions: let y be the number of 6’s in 1000 independent rolls of a par- ticular real die, which may be unfair. Let θ be the probability that the die lands on ‘6.’ Suppose your prior distribution for θ is as follows:
Pr(θ = 1/12) = 0.25, Pr(θ = 1/6) = 0.5, Pr(θ = 1/4)= 0.25.
(a) Using the normal approximation for the conditional distributions, p(y|θ), sketch your approximate prior predictive distribution for y
Since y is summation of a large number (1000) of independent variables, according to CLT it can be apporoximated by a normal distribution:
\[ p(y|\theta) \sim N(n\theta, n\theta(1-\theta))\] Thus we can formulate the approximate prior predictive distribution for y, as follows:
\[ p(y)= \sum_\theta p(y,\theta)=\sum_\theta p(y|\theta)p(\theta) = 0.25N(y|1000\times\frac{1}{4}, 1000\times\frac{1}{4}\times\frac{3}{4}) + 0.5N(y|1000\times\frac{1}{6}, 1000\times\frac{1}{6}\times\frac{5}{6}) + 0.25N(y|1000\times\frac{1}{12}, 1000\times\frac{1}{12}\times\frac{11}{12})\]
Using the following R code I will sketch the distribution:
theta <- c(1/12,1/6,1/4)
prob_theta <- c(0.25,0.5,0.25)
n=1000
means <- theta*n
sds <- sqrt(n*theta*(1-theta))
x <- seq(-4,4,length.out = 10000)*sum(sds) + mean(means)
y <- prob_theta[1]*dnorm(x,means[1],sds[1]) + prob_theta[2]*dnorm(x,means[2],sds[2]) + prob_theta[3]*dnorm(x,means[3],sds[3])
plot(x,y,type = "l",lty=1)
(b) Give approximate 5%, 25%, 50%, 75%, and 95% points for the distribution of y. (Be careful here: y does not have a normal distribution, but you can still use the normal distribution as part of your analysis.)
As we can see in the plot, 25%, 50% and the other 25% of the density of this distribution are in almost exclusive pieces of x axis and each piece have a normal distribution with certain parameters. So we can approximate the quantiles this way:
\[ q_{5\%}= q_{20\%-N1}= 75.98 \] \[ q_{25\%}= q_{99\%-N1}= 103.66\] \[ q_{50\%}= q_{50\%-N2}= 166.667\] \[ q_{75\%}= q_{99\%-N2}= 194.83\]
\[ q_{95\%}= q_{80\%-N3}= 261.52\]
The answers are calculated using these lines of code:
q20=qnorm(0.2,means[1],sds[1])
print(q20)
## [1] 75.9775
q25=qnorm(0.99,means[1],sds[1])
print(q25)
## [1] 103.6658
q50=qnorm(0.5,means[2],sds[2])
print(q50)
## [1] 166.6667
q75=qnorm(0.99,means[2],sds[2])
print(q75)
## [1] 194.0829
q95=qnorm(0.8,means[3],sds[3])
print(q95)
## [1] 261.5244
Discrete sample spaces: suppose there are N cable cars in San Francisco, numbered sequentially from 1 to N. You see a cable car at random; it is numbered 203. You wish to estimate N. (See Goodman, 1952, for a discussion and references to several versions of this problem, and Jeffreys, 1961, Lee, 1989, and Jaynes, 2003, for Bayesian treatments.)
(a) Assume your prior distribution on N is geometric with mean 100; that is,
p(N) = (1/100)(99/100)N−1, for N = 1,2,….
What is your posterior distribution for N?
\[ p(N|y)=\frac{p(y|N)p(N)}{p(y)}\] \[p(y|N)=\begin{cases} 0 & y>N \\ \frac{1}{N} & y\leq N \end{cases}\] \[ p(y)=\sum_N p(y,N)=\sum_{N=y}^\infty p(y|N)p(N)=\sum_{N=203}^\infty (\frac{1}{N}\times\frac{1}{100}\times(\frac{99}{100})^{N-1})=0.00047\] Thus: \[ p(N|y=203)=\frac{1}{0.047}(\frac{1}{N}\times(\frac{99}{100})^{N-1})\] (b) What are the posterior mean and standard deviation of N? (Sum the infinite series analytically or approximate them on the computer.) \[E(N|y)=\sum_N Np(N|y)=\frac{1}{0.047}\sum_{N=203}^\infty(0.99)^{N-1}=279.39\] \[ Var(N|y)=\sum_N (N-279.39)^2p(N|y)=\frac{1}{0.047}\sum_{N=203}^\infty \frac{(N-279.39)^2}{N}(0.99)^{N-1}=6401.340\\ \sigma_{N|y}=80.01\] (c) Choose a reasonable ‘noninformative’ prior distribution for N and give the resulting posterior distribution, mean, and standard deviation for N.
At first I tried uniform distribution but it does not work. Since in this case the denominator of posterior distribution does not converge. suppose we choose poisson distribution with parameter 100 as a noninformative prior:
\[p(N)=\frac{100^Ne^{-100}}{N!}\] \[p(N|y)=\frac{\frac{1}{N}\times \frac{100^Ne^{-100}}{N!}}{p(y)}\] \[p(y)=\sum_{N=203}^\infty \frac{1}{N}.\frac{100^Ne^{-100}}{N!}=k\] This summation converges and we assue the answer is k \[P(N|y)=\frac{1}{k} \times \frac{1}{N}\times \frac{100^Ne^{-100}}{N!}\] \[E(N|y)=\sum_N Np(N|y)=\frac{1}{k}\sum_{N=203}^\infty \frac{100^Ne^{-100}}{N!}=\frac{1}{k}\times5.453751e-20\] \[ Var(N|y)=\sum_N (N-\frac{1}{k}\times5.453751e-20)^2p(N|y)=\frac{1}{k}\sum_{N=203}^\infty \frac{(N-\frac{1}{k}\times5.453751e-20)^2}{N}\times\frac{100^Ne^{-100}}{N!}\] Since I do not know k, I cannot calculate amount of the standard deviation
You wish to use a binomial model to make inferences about the odds of a ‘success’,ω,and you decide to use a uniform prior on the odds. Is this a proper or an improper prior? Why? If you were to estimate the probability of a success, θ, what prior would you have to use on θ so your inferences would be equivalent to those based on the odds?
\[ \omega = \frac{\theta}{1-\theta}\] This is an improper prior because range of the omega is zero to positive infinity. So every amount you choose for C (as the density of uniform distribution), the integration will be greater than 1 and the function will not meet the conditions of a probabiltiy density function.
\[p_{\omega}(\omega)=C\] \[p_{\theta}(\theta)=|J|\times p_{\omega}(\frac{\theta}{1-\theta})\] |J| simply equals to: \[\frac{\partial \omega}{\partial \theta}=\frac{1}{(1-\theta)^2}\] So: \[ p_{\theta}(\theta)=\frac{1}{(1-\theta)^2}C\] If we use the prior above, both inferences will be equivalent.