Statistical Models

Amy Jiang

Janurary 4, 2015

---

A learning script for plotting each statistical models’ pdf and cdf, together with some basic statistical measures like mean, variance, std Variance in comparison with the known statistical functions.

---

BETA Distribution

Beta distribution PDF Chart

set.seed(5)
x <- seq(0, 1, length=200)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

data1 <- dbeta(x, 0.5, 0.5)
data2 <- dbeta(x, 5, 1)
data3 <- dbeta(x, 1, 3)
data4 <- dbeta(x, 2, 2)
data5 <- dbeta(x, 2, 5)

plot(x, data1, ylim=c(0, 2.5), ylab="PDF", type="l", lwd=2, 
     col="red", yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, data2, col="blue", lwd="2")
lines(x, data3, col="green", lwd="2")
lines(x, data4, col="brown", lwd="2")
lines(x, data5, col="purple", lwd="2")

legend(0.4, 2.5, c(expression(paste(alpha, "=", beta, "=0.5")),
                   expression(paste(alpha, "=5, ", beta, "=1")),
                   expression(paste(alpha, "=1, ", beta, "=3")),
                   expression(paste(alpha, "=2, ", beta, "=2")),
                   expression(paste(alpha, "=2, ", beta, "=5"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=.75)

Beta distribution CDF Chart

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

pdata1 <- pbeta(x, 0.5, 0.5)
pdata2 <- pbeta(x, 5, 1)
pdata3 <- pbeta(x, 1, 3)
pdata4 <- pbeta(x, 2, 2)
pdata5 <- pbeta(x, 2, 5)

plot(x, pdata1, ylab="CDF", type="l", lwd=2, col="red", yaxs="i", 
     xaxs="i", frame.plot=FALSE)
lines(x, pdata2, col="blue", lwd="2")
lines(x, pdata3, col="green", lwd="2")
lines(x, pdata4, col="brown", lwd="2")
lines(x, pdata5, col="purple", lwd="2")

legend(0, 1, c(expression(paste(alpha, "=", beta, "=0.5")),
                   expression(paste(alpha, "=5, ", beta, "=1")),
                   expression(paste(alpha, "=1, ", beta, "=3")),
                   expression(paste(alpha, "=2, ", beta, "=2")),
                   expression(paste(alpha, "=2, ", beta, "=5"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=0.75)

Beta distribution Mean, Variance, Standard Deviation

betaMean  <- function(a, b) {
    a/(a+b)
}
betaVar  <- function(a, b){
    a*b/((a+b)^2*(a+b+1))
}
betaSd <- function(a, b){
    sqrt(a*b/((a+b)^2*(a+b+1)))
}

a <- 0.5; b <- 0.5
modelMean <- betaMean(a, b) 
realMean <- mean(rbeta(1000, a, b)) 
modelVar <- betaVar(a, b)  
realVar <- var(rbeta(1000, a, b)) 
modelSd <- betaSd(a, b) # 
realSd <- sd(rbeta(1000, a, b)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ", realSd)

## Modeled mean:  0.5 ; Real mean:  0.4948387 
## Modeled variance:  0.125 ; Real variance:  0.1224283 
## Modeled standard variation:  0.3535534 ; Real standard variation:  0.3637433

---

Binomial Distribution

Binomial distribution PDF Chart

myBinom  <- function(n, p) {
    x <- seq(0, n, length=50)    
    y  <- (factorial(n)/(factorial(n-x)*factorial(x)))*p^x*(1-p)^(n-x) 
    df <- data.frame(x, y)
    colnames(df)  <- c("x", "y")
    return(df)
}

data1 <- myBinom(20, 0.5)
data2 <- myBinom(20, 0.7)
data3 <- myBinom(40, 0.5)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 
plot(data1, ylim=c(0, 0.25), xlim=c(0,40), ylab="PDF", type="n", yaxs="i", 
     xaxs="i", bty="n", frame.plot=FALSE)
points(data1$x, data1$y, col="red", pch=21, bg="red")
points(data2$x, data2$y, col="green", pch=21, bg="green")
points(data3$x, data3$y, col="blue", pch=21, bg="blue")

legend(25, 0.25, c("p=0.5, n=20","p=0.7, n=20", "p=0.5, n=40"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green'), cex=.75)

Beta distribution CDF Chart

pdata1 <- cumsum(data1$y); pdata1 <- pdata1/max(pdata1)
pdata2 <- cumsum(data2$y); pdata2 <- pdata2/max(pdata2)
pdata3 <- cumsum(data3$y); pdata3 <- pdata3/max(pdata3)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(seq(0, 20, length=50), pdata1, ylab="CDF", xlim=c(0, 40), 
     ylim=c(0.0, 1.1), lwd=2, col="red", pch=21, bg="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
points(seq(0,20, length=50), pdata2, col="green", pch=21, bg="green")
points(seq(0,40, length=50), pdata3, col="blue", pch=21, bg="blue")

legend(30, 0.4, c("p=0.5, n=20","p=0.7, n=20", "p=0.5, n=40"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green'), cex=.75)

Binomial distribution Mean, Variance, Standard Deviation

binomMean  <- function(n, p) {
    n*p
}
binomVar  <- function(n, p){
    n*p*(1-p)
}
binomSd <- function(a, b){
    sqrt(n*p*(1-p))
}

n <- 20; p <- 0.5
modelMean <- binomMean(n, p) 
realMean <- mean(rbinom(1000, n, p)) 
modelVar <- binomVar(n, p)  
realVar <- var(rbinom(1000, n, p)) 
modelSd <- binomSd(n, p) 
realSd <- sd(rbinom(1000, n, p)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ", realSd)

## Modeled mean:  10 ; Real mean:  10.023 
## Modeled variance:  5 ; Real variance:  5.258775 
## Modeled standard variation:  2.236068 ; Real standard variation:  2.218007

---

Chi-squared Distribution

Chi-squared distribution PDF Chart

x <- seq(0, 8, length=200)

data1 <- dchisq(x, 1)
data2 <- dchisq(x, 2)
data3 <- dchisq(x, 4)
data4 <- dchisq(x, 6)
data5 <- dchisq(x, 9)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(x, data1, ylim=c(0, 0.5), ylab="PDF", type="l", lwd=2, col="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, data2, col="blue", lwd="2")
lines(x, data3, col="green", lwd="2")
lines(x, data4, col="brown", lwd="2")
lines(x, data5, col="purple", lwd="2")

legend(6, 0.5, c("k=1","k=2","k=4","k=6","k=9"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=.75)

Chi-squared distribution CDF Chart

pdata1 <- pchisq(x, 1)
pdata2 <- pchisq(x, 2)
pdata3 <- pchisq(x, 4)
pdata4 <- pchisq(x, 6)
pdata5 <- pchisq(x, 9)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(x, pdata1, ylab="CDF", ylim=c(0, 1), type="l", lwd=2, col="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, pdata2, col="blue", lwd="2")
lines(x, pdata3, col="green", lwd="2")
lines(x, pdata4, col="brown", lwd="2")
lines(x, pdata5, col="purple", lwd="2")

legend(6.5, 0.4, c("k=1","k=2","k=4","k=6","k=9"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=0.75)

Chi-squared distribution Mean, Variance, Standard Deviation

chisqMean  <- function(k) {
    k
}
chisqVar  <- function(k){
    2*k
}
chisqSd <- function(k){
    sqrt(2*k)
}

k <- 1
modelMean <- chisqMean(k) 
realMean <- mean(rchisq(1000, k)) 
modelVar <- chisqVar(k)  
realVar <- var(rchisq(1000, k)) 
modelSd <- chisqSd(k) # 
realSd <- sd(rchisq(1000, k)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ", realSd)

## Modeled mean:  1 ; Real mean:  0.9498187 
## Modeled variance:  2 ; Real variance:  2.109575 
## Modeled standard variation:  1.414214 ; Real standard variation:  1.413127

---

Exponential Distribution

Exponential distribution PDF Chart

x <- seq(0, 5, length=200)

data1 <- dexp(x, 0.5)
data2 <- dexp(x, 1)
data3 <- dexp(x, 1.5)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(x, data1, ylim=c(0, 1.6), ylab="PDF", type="l", lwd=2, col="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, data2, col="blue", lwd="2")
lines(x, data3, col="green", lwd="2")

legend(3.5, 1.5, c(expression(paste(lambda, "=0.5")),
                 expression(paste(lambda, "=1.0")),
                 expression(paste(lambda, "=1.5"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=.75)

Chi-squared distribution CDF Chart

pdata1 <- pexp(x, 0.5)
pdata2 <- pexp(x, 1)
pdata3 <- pexp(x, 1.5)


par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(x, pdata1, ylab="CDF", ylim=c(0, 1), type="l", lwd=2, col="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, pdata2, col="blue", lwd="2")
lines(x, pdata3, col="green", lwd="2")

legend(3.5, 0.4, c(expression(paste(lambda, "=0.5")),
                 expression(paste(lambda, "=1.0")),
                 expression(paste(lambda, "=1.5"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=0.75)

Chi-squared distribution Mean, Variance, Standard Deviation

expMean  <- function(lambda) {
    1/lambda
}
expVar  <- function(lambda){
    1/lambda^2
}
expSd <- function(lambda){
    1/lambda
}

lambda <- 1
modelMean <- expMean(lambda) 
realMean <- mean(rexp(1000, lambda)) 
modelVar <- expVar(lambda)  
realVar <- var(rexp(1000, lambda)) 
modelSd <- expSd(lambda) # 
realSd <- sd(rexp(1000, lambda)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ", realSd)

## Modeled mean:  1 ; Real mean:  1.038569 
## Modeled variance:  1 ; Real variance:  0.9252989 
## Modeled standard variation:  1 ; Real standard variation:  1.024291

---

Gamma Distribution

Gamma distribution PDF Chart

x <- seq(0, 20, length=200)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

data1 <- dgamma(x, 1.0, 0.5)
data2 <- dgamma(x, 2.0, 0.5)
data3 <- dgamma(x, 3.0, 0.5)
data4 <- dgamma(x, 5.0, 1.0)
data5 <- dgamma(x, 9.0, 2)

plot(x, data1, ylim=c(0, 0.5), ylab="PDF", type="l", lwd=2, 
     col="red", yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, data2, col="blue", lwd="2")
lines(x, data3, col="green", lwd="2")
lines(x, data4, col="brown", lwd="2")
lines(x, data5, col="purple", lwd="2")

legend(15, 0.5, c(expression(paste(alpha, "=1, ", beta, "=0.5")),
                   expression(paste(alpha, "=2, ", beta, "=0.5")),
                   expression(paste(alpha, "=3, ", beta, "=0.5")),
                   expression(paste(alpha, "=5, ", beta, "=1")),
                   expression(paste(alpha, "=9, ", beta, "=0.2"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=.75)

Gammaa distribution CDF Chart

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

pdata1 <- pgamma(x, 1, 0.5)
pdata2 <- pgamma(x, 2, 0.5)
pdata3 <- pgamma(x, 3, 0.5)
pdata4 <- pgamma(x, 5, 1)
pdata5 <- pgamma(x, 9, 2)

plot(x, pdata1, ylab="CDF", ylim=c(0,1), type="l", lwd=2, col="red", yaxs="i", 
     xaxs="i", frame.plot=FALSE)
lines(x, pdata2, col="blue", lwd="2")
lines(x, pdata3, col="green", lwd="2")
lines(x, pdata4, col="brown", lwd="2")
lines(x, pdata5, col="purple", lwd="2")

legend(15, 0.5, c(expression(paste(alpha, "=1, ", beta, "=0.5")),
                   expression(paste(alpha, "=2, ", beta, "=0.5")),
                   expression(paste(alpha, "=3, ", beta, "=0.5")),
                   expression(paste(alpha, "=5, ", beta, "=1")),
                   expression(paste(alpha, "=9, ", beta, "=0.2"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"),
       cex=0.75)

Gamma distribution Mean, Variance, Standard Deviation

gammaMean  <- function(a, b) {
    a/b
}
gammaVar  <- function(a, b){
    a/b^2
}
gammaSd <- function(a, b){
    sqrt(a/b^2)
}

a <- 1; b <- 0.5
modelMean <- gammaMean(a, b) 
realMean <- mean(rgamma(1000, a, b)) 
modelVar <- gammaVar(a, b)  
realVar <- var(rgamma(1000, a, b)) 
modelSd <- gammaSd(a, b) # 
realSd <- sd(rgamma(1000, a, b)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ", realSd)

## Modeled mean:  2 ; Real mean:  1.896239 
## Modeled variance:  4 ; Real variance:  4.570306 
## Modeled standard variation:  2 ; Real standard variation:  1.847329

---

Geometric Distribution

Geometric distribution PDF Chart

myGeom  <- function(x, p) {
    x <- seq(0, x, length=25)    
    y  <- (1-p)^x*p
    df <- data.frame(x, y)
    colnames(df)  <- c("x", "y")
    return(df)
}

data1 <- myGeom(10, 0.2)
data2 <- myGeom(10, 0.5)
data3 <- myGeom(10, 0.8)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 
plot(data1, ylim=c(0, 1.0), xlim=c(0,10), ylab="PDF", type="n", yaxs="i",
     xaxs="i", bty="n", frame.plot=FALSE)
points(data1$x, data1$y, col="red", pch=21, bg="red")
lines(data1$x, data1$y, col="red")
points(data2$x, data2$y, col="green", pch=21, bg="green")
lines(data2$x, data2$y, col="green")
points(data3$x, data3$y, col="blue", pch=21, bg="blue")
lines(data3$x, data3$y, col="blue")

legend(6, 1.0, c("k=10, p=0.2","k=10, p=0.5", "k=10, p=0.8"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green'), cex=.75)

Geometric distribution CDF Chart

pdata1 <- cumsum(data1$y); pdata1 <- pdata1/max(pdata1)
pdata2 <- cumsum(data2$y); pdata2 <- pdata2/max(pdata2)
pdata3 <- cumsum(data3$y); pdata3 <- pdata3/max(pdata3)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(seq(0, 10, length=25), pdata1, ylab="CDF", xlim=c(0, 10), 
     ylim=c(0.0, 1.1), lwd=2, col="red", pch=21, bg="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
lines(seq(0,10, length=25), pdata1, col="red")
points(seq(0,10, length=25), pdata2, col="green", pch=21, bg="green")
lines(seq(0,10, length=25), pdata2, col="green")
points(seq(0,10, length=25), pdata3, col="blue", pch=21, bg="blue")
lines(seq(0,10, length=25), pdata3, col="blue")

legend(6, 0.4, c("k=10, p=0.2","k=10, p=0.5", "k=10, p=0.8"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green'), cex=.75)

Geometric distribution Mean, Variance, Standard Deviation

geomMean  <- function(p) {
    (1-p)/p
}
geomVar  <- function(p){
    (1-p)/p^2
}
geomSd <- function(p){
    sqrt((1-p)/p^2)
}

p <- 0.2
modelMean <- geomMean(p) 
realMean <- mean(rgeom(1000, p)) 
modelVar <- geomVar(p)  
realVar <- var(rgeom(1000, p)) 
modelSd <- geomSd(p) 
realSd <- sd(rgeom(1000, p)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ",
    realSd)

## Modeled mean:  4 ; Real mean:  3.848 
## Modeled variance:  20 ; Real variance:  21.52527 
## Modeled standard variation:  4.472136 ; Real standard variation:  4.281752

---

Normal Distribution

Normal distribution PDF Chart

x <- seq(-5, 5, length=200)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

data1 <- dnorm(x, 0, sqrt(0.2))
data2 <- dnorm(x, 0, sqrt(1))
data3 <- dnorm(x, 0, sqrt(5))
data4 <- dnorm(x, -2, sqrt(0.5))

plot(x, data1, ylim=c(0, 1), xlim=c(-5,5), ylab="PDF", type="l", lwd=2, 
     col="red", yaxs="i", xaxs="i", frame.plot=FALSE)
lines(x, data2, col="blue", lwd="2")
lines(x, data3, col="green", lwd="2")
lines(x, data4, col="brown", lwd="2")

legend(2, 1, c(expression(paste(mu, "=0, ", sigma^2, "=0.2")),
                   expression(paste(mu, "=0, ", sigma^2, "=1")),
                   expression(paste(mu, "=0, ", sigma^2, "=5")),
                   expression(paste(mu, "=-2, ", sigma^2, "=0.5"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"), cex=.75)

Normal distribution CDF Chart

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

pdata1 <- pnorm(x, 0, sqrt(0.2))
pdata2 <- pnorm(x, 0, sqrt(1))
pdata3 <- pnorm(x, 0, sqrt(5))
pdata4 <- pnorm(x, -2, sqrt(0.5))

plot(x, pdata1, ylab="CDF", ylim=c(0,1), xlim=c(-5,5), type="l", lwd=2, col="red", yaxs="i", 
     xaxs="i", frame.plot=FALSE)
lines(x, pdata2, col="blue", lwd="2")
lines(x, pdata3, col="green", lwd="2")
lines(x, pdata4, col="brown", lwd="2")

legend(2, 0.4, c(expression(paste(mu, "=0, ", sigma^2, "=0.2")),
                   expression(paste(mu, "=0, ", sigma^2, "=1")),
                   expression(paste(mu, "=0, ", sigma^2, "=5")),
                   expression(paste(mu, "=-2, ", sigma^2, "=0.5"))),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green',' brown', "purple"),
       cex=0.75)

Normal distribution Mean, Variance, Standard Deviation

normalMean  <- function(mu) {
    mu
}
normalVar  <- function(sigma){
    sigma ^ 2
}
normalSd <- function(sigma){
    sigma
}

mu <- 10; sigma <- 1
modelMean <- normalMean(mu) 
realMean <- mean(rnorm(1000, mu, sigma)) 
modelVar <- normalVar(sigma)  
realVar <- var(rnorm(1000, mu, sigma)) 
modelSd <- normalSd(sigma) 
realSd <- sd(rnorm(1000, mu, sigma)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ", realSd)

## Modeled mean:  10 ; Real mean:  9.990457 
## Modeled variance:  1 ; Real variance:  0.9818876 
## Modeled standard variation:  1 ; Real standard variation:  1.031042

---

Poisson Distribution

Piosson distribution PDF Chart

myPoisson  <- function(k, lambda) {
    x <- seq(0, k, length=k+1)    
    y  <- lambda^x*exp(-lambda)/factorial(x)
    df <- data.frame(x, y)
    colnames(df)  <- c("x", "y")
    return(df)
}

data1 <- myPoisson(20, 1)
data2 <- myPoisson(20, 4)
data3 <- myPoisson(20, 10)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 
plot(data1, ylim=c(0, 0.4), xlim=c(0,20), ylab="PDF", type="n", yaxs="i",
     xaxs="i", bty="n", frame.plot=FALSE)
points(data1$x, data1$y, col="red", pch=21, bg="red")
lines(data1$x, data1$y, col="red")
points(data2$x, data2$y, col="green", pch=21, bg="green")
lines(data2$x, data2$y, col="green")
points(data3$x, data3$y, col="blue", pch=21, bg="blue")
lines(data3$x, data3$y, col="blue")

legend(12, 0.4, c("k=20, lambda=1","k=20, lambda=4", "k=20, lambda=10"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green'), cex=.75)

Poisson distribution CDF Chart

pdata1 <- cumsum(data1$y); pdata1 <- pdata1/max(pdata1)
pdata2 <- cumsum(data2$y); pdata2 <- pdata2/max(pdata2)
pdata3 <- cumsum(data3$y); pdata3 <- pdata3/max(pdata3)

par(oma=c(1, 1, 1, 1))
par(mar=c(2, 4, 2, 8)+0.1) 

plot(seq(0, 20, length=21), pdata1, ylab="CDF", xlim=c(0, 20), 
     ylim=c(0.0, 1.1), lwd=2, col="red", pch=21, bg="red", 
     yaxs="i", xaxs="i", frame.plot=FALSE)
lines(seq(0,20, length=21), pdata1, col="red")
points(seq(0,20, length=21), pdata2, col="green", pch=21, bg="green")
lines(seq(0,20, length=21), pdata2, col="green")
points(seq(0,20, length=21), pdata3, col="blue", pch=21, bg="blue")
lines(seq(0,20, length=21), pdata3, col="blue")

legend(12, 0.4, c("k=20, lambda=1","k=20, lambda=4", "k=20, lambda=10"),
       x.intersp = 0.5, y.intersp = 1, bty="n",
       lty=c(1, 1, 1, 1), col=c('red', 'blue', 'green'), cex=.75)

Poisson distribution Mean, Variance, Standard Deviation

poissonMean  <- function(lambda) {
    lambda
}
poissonVar  <- function(lambda){
    lambda
}
poissonSd <- function(lambda){
    sqrt(lambda)
}

lambda <- 1
modelMean <- poissonMean(lambda) 
realMean <- mean(rpois(1000, lambda)) 
modelVar <- poissonVar(lambda)  
realVar <- var(rpois(1000, lambda)) 
modelSd <- poissonSd(lambda) 
realSd <- sd(rpois(1000, lambda)) 

cat("Modeled mean: ", modelMean, "; Real mean: ", realMean, "\n")
cat("Modeled variance: ", modelVar, "; Real variance: ", realVar, "\n")
cat("Modeled standard variation: ", modelSd, "; Real standard variation: ",
    realSd)

## Modeled mean:  1 ; Real mean:  0.939 
## Modeled variance:  1 ; Real variance:  1.161417 
## Modeled standard variation:  1 ; Real standard variation:  1.072775