# ex2.R
# Demonstrate simple statistics in R
# to be used for copy-paste into browser at https://rdrr.io/snippets/

# 2019-02-03 15:56 (Gunnar Taraldsen)
# 2019-02-03 21:15 (GT)

######## Bernoulli ++ ###########################

require(distr)
iMaxI = 4
nProb = 0.5
aI = seq(2,iMaxI,1)

Binom = Binom(prob=nProb, size=1)
USum = Binom(prob=nProb, size=1)/2

for (i in aI) {

  B = Binom(prob=nProb, size=1)

  Binom = Binom + B
  USum = USum + B/2^i

}

USum

# frame()
u = seq(0, 1, .01)
# lines(u, p.r(USum)(u),lty=2, col="red")
plot(u, p.r(USum)(u),lty=2, type="l")

plot(Binom)
plot(USum)



######## Arbitrary density ###########################

library(distr)

gtpdf = function(x) (1/3) * (1 + 4*x)  # probability density function
gtdist = AbscontDistribution(d=gtpdf, low1=0.0, up1=1.0)  # signature for a dist with pdf ~ gtpdf

# gtpdf = function(x) (1/3) * (1 + 4*x) * (x > 0) * (x < 1)  # probability density function
# gtdist = AbscontDistribution(d=gtpdf)  # signature for a dist with pdf ~ gtpdf

plot(gtdist)

rgtdist = r(gtdist)                 # function to create random variates from p

set.seed(1)                      # for reproduceable example
X = rgtdist(1000)                 # sample from X ~ p
x = seq(0, 1, .01)
hist(X, freq=F, breaks=50, xlim=c(0,1))
lines(x,gtpdf(x),lty=2, col="red")


######## Joint from marginal and conditional  ###########################

library(distr)

# Marginal distribution
xpdf = function(x) (1/3) * (1 + 4*x)  # probability density function
xdist = AbscontDistribution(d=xpdf, low1=0.0, up1=1.0)  # signature for a dist with pdf ~ xpdf
plot(xdist)
rxdist = r(xdist)                 # function to create random variates from xdist

# set.seed(1)                      # for reproduceable example
aX = rxdist(1000)                 # sample from X ~ p
aY = aX; # To have something
i=0

for (x in aX) {
  i=i+1;
  ypdfCond = function(y) (2*y + 4*x)/(1 + 4*x)  # conditional probability density function
  ydistCond = AbscontDistribution(d=ypdfCond, low1=0.0, up1=1.0)  # signature for a dist with pdf ~ ypdf
  rydistCond = r(ydistCond)                 # function to create random variates from ydist
  aY[i] = rydistCond(1) 
}

ypdf = function(y) (2/3)*(y + 1)  # conditional probability density function
y = seq(0, 1, .01)
hist(aY, freq=F, breaks=50, xlim=c(0,1))
lines(y,ypdf(y),lty=2, col="red")

# EOF ex2.R