## com25b of september 2021,
#  stk 4021, BDA exercise 2.10, street cars in San Francisco 

#  first: a single observation 
#  afterwards i extend this to say n = 3 observations 
#  then i ask students to do with with n = 3 and with 203, 157, 222 
#  later we also extend the exercise to allow My Own Prior 

y1 <- 203  
ymore <- c(203, 157, 222) 
ymax <- max(ymore) 

# I edit out and in and note the results

Nmax <- 10^4 # more than enough 
Nval <- 1:Nmax
p0 <- 0.01
prior <- p0*(1-p0)^(Nval-1)

# have a look at the prior: 
plot(Nval[1:500],prior[1:500])

uposterior1 <- prior*(1/Nval)*(Nval >= y1)
const1 <- sum(uposterior1)
posterior1 <- uposterior1/const1

uposterior3 <- prior*(1/Nval)^3*(Nval >= ymax)
const3 <- sum(uposterior3)
posterior3 <- uposterior3/const3

par(cex.lab=1.4142)
par(bty="l")

matplot(Nval[1:600], 
   cbind(prior[1:600], posterior1[1:600], posterior3[1:600]),
   xlab="N", ylab="prior, post1, post3",
   type="l",lwd=2.0,col=c(1,2,4))

N0 <- sum(Nval*prior)
sigma0 <- sqrt( sum((Nval - N0)^2*prior) )
print(c(N0, sigma0))

Npost1 <- sum(Nval*posterior1)
sigmapost1 <- sqrt( sum((Nval - Npost1)^2*posterior1) )

Npost3 <- sum(Nval*posterior3)
sigmapost3 <- sqrt( sum((Nval - Npost3)^2*posterior3) )

print(c(N0, sigma0))
print(c(Npost1, sigmapost1))
print(c(Npost3, sigmapost3))

# looking at cumulative prior, cumulative posterior,
# also to find quantiles: 

cumprior <- cumsum(prior)
cumposterior1 <- cumsum(posterior1) 
cumposterior3 <- cumsum(posterior3) 

low0 <- min(Nval[cumprior >= 0.05])
mid0 <- min(Nval[cumprior >= 0.50]) 
up0 <- min(Nval[cumprior >= 0.95])

low1 <- min(Nval[cumposterior1 >= 0.05])
mid1 <- min(Nval[cumposterior1 >= 0.50]) 
up1 <- min(Nval[cumposterior1 >= 0.95])

low3 <- min(Nval[cumposterior3 >= 0.05])
mid3 <- min(Nval[cumposterior3 >= 0.50]) 
up3 <- min(Nval[cumposterior3 >= 0.95])

check0 <- c(low0,mid0,up0,N0,sigma0)
check1 <- c(low1,mid1,up1,Npost1,sigmapost1)
check3 <- c(low3,mid3,up3,Npost3,sigmapost3)

print( round(check0,3) )
print( round(check1,3) )
print( round(check3,3) )

#  low med  up   mean     stdev 
#   6   69  299  100.000  99.499 before data 
# 206  254  438  279.089  79.965 after 1 data point
# 224  255  381  272.578  55.090 after 3 data points