تمرین هفتم درس R

تمرین درس R سارا اکبر ی در R  توابع بسیاری برای محاسبه چندک ها مقادیر , مقادیر احتمال توزیع های مختلف و تولید اعداد تصادفی وجود دارد . به طوری که کافیست در توابع rfunx(n,param), qfunx(p,param),pfunx(x,param) , dfunx(x,param) کد توزیع و پارامترهای مربوطه را به جای عبارات fun و param بنویسیم . در زیر یک سری از توابع احتمل همراه با مثال را نوشته ایم : تابع هندسی – دستور زیر توزیع تابع هندسی رابررسی می کند:

Geom

به مثال زیر برای تابع هندسی است توجه کنید :

Examples

qgeom((1:9)/10, prob = .2)

[1]  0  0  1  2  3  4  5  7 10

Ni <- rgeom(20, prob = 1/4); table(factor(Ni, 0:max(Ni)))

0 1 2 3 4 5 6 7 8 9

8 2 2 2 2 1 1 1 0 1

تابع فوق هندسی - دستور زیر توزیع تابع فوق هندسی رابررسی می کند:

Hyper

به مثال زیر برای تابع فوق هندسی است توجه کنید :

Examples

m <- 10; n <- 7; k <- 8

x <- 0:(k+1)

rbind(phyper(x, m, n, k), dhyper(x, m, n, k))

all(phyper(x, m, n, k) == cumsum(dhyper(x, m, n, k)))  # FALSE

## but error is very small:

signif(phyper(x, m, n, k) - cumsum(dhyper(x, m, n, k)), digits = 3)

تابع لجستیک - دستور زیر توزیع تابع لجستیک رابررسی می کند:

Logis

به مثال زیر برای تابع لجستیک است توجه کنید

Examples

var(rlogis(4000, 0, scale = 5))  # approximately (+/- 3)

pi^2/3 * 5^2

تابع نرمالlog- دستور زیر توزیع تابع نرمال log  رابررسی می کند:

Lnorm

به مثال زیر برای تابع نرمالlog است توجه کنید

Examples

(x1 <- cbind(1, 1:10))

norm(x1)

norm(x1, “I”)

norm(x1, “M”)

stopifnot(all.equal(norm(x1, “F”),

                    sqrt(sum(x1^2))))

hilbert <- function(n) { i <- 1:n; 1 / outer(i – 1, i, “+”) }

h9 <- hilbert(9)

## all 5 types of norm:

(nTyp <- eval(formals(base::norm)$type))

sapply(nTyp, norm, x = h9)

تابع دو جمله ای منغی - دستور زیر توزیع تابع دو جمله ای منفی رابررسی می کند:

Nbinom

به مثال زیر برای تابع دوجمله ای منفی است توجه کنید

Examples

require(graphics)

# Compute P(45 < X < 55) for X Binomial(100,0.5)

sum(dbinom(46:54, 100, 0.5))

## Using "log = TRUE" for an extended range :

n <- 2000

k <- seq(0, n, by = 20)

plot (k, dbinom(k, n, pi/10, log = TRUE), type = "l", ylab = "log density",

      main = "dbinom(*, log=TRUE) is better than  log(dbinom(*))")

lines(k, log(dbinom(k, n, pi/10)), col = "red", lwd = 2)

## extreme points are omitted since dbinom gives 0.

mtext("dbinom(k, log=TRUE)", adj = 0)

mtext("extended range", adj = 0, line = -1, font = 4)

mtext("log(dbinom(k))", col = "red", adj = 1)

 محاسبه احتمال مقادیر در توزیع فراوانی ها توزیع فراوانی هر متغیرتصادفی نشان می دهد که پراکندگی مقادیر آن متغیر تصادفی چگونه است . برای محاسبه احتمال وقوع یک متغیر تصادفی دارای توزیع خاص ، توابع گوناگونی در نرم افزار R وجود دارد . مقادیر این احتمالات معمولا در تجزیه های آماری مورد نیاز هستند و در پایان کتب آماری به صورت جداول که مورد استفاده فراوان قرار می گیرند ، جداول Z،T، X^2 و F می باشند . توزیع نرمال فرض کنید فراوانی نمرات پایان ترم در آزمون یکی از دروس دارای توزیع نرمال باشد پنانپه میانگین دانشجویان برابر با 72 و انحراف معیار کلاس 2/15 باشد و بخواهیم درصد دانشجویانی که نمره بیش از 84 کسب نموده اند را محاسبه کنیم باید از تابع pnorm() توزیع Z به صورت زیر استفاده نماییم .

> pnorm(84,mean=72,sd=15.2,lower.tail=FALSE)

[1] 0.2149176

همانگونه که مشاهده می شود حدود 5/21 درصد دانشجویان کلاس نمراتی بیش از 84 کسب نموده اند . برهان lower.tail=FALSE  به برنامه می گوید تا احتمال مقادیر بیش از 84 را محاسبه نمایید . اگر این برهان برابر TRUE شود احتمال مقادیر کوچکتر و مساوی 84 را محاسبه خواهد نمود . از توابع دیگر مجموعه توابع نرمال می توان به dnorm() (محاسبه چندک برای احتمالات خاص ) Qnorm() (ایجاد تابع چندکی ) ، rnorm() (ایجاد داده های نرمال تصادفی ) اشاره نمود . به عنوان مثال چنانچه بخواهیم 30 عدد تصادفی که از یک توزیع نرمال با میانگین 17 و انحراف معیار 3 استخراج شده باشند به دست آورید از تابع rnorm() به صورت زیر استفاده می نمایید :

> rnorm(30,mean=17,sd=3)

 [1] 17.460853 15.083290 19.179011 15.886316 25.164220 19.188159 14.817927

 [8] 17.193516 16.317907 22.447368 19.077821 14.528341 17.389222 14.649513

[15] 21.529489 19.409211 20.729602 15.585371 19.445909 17.016549  9.041588

[22] 14.774410 14.370598 16.132691 19.803648 16.077984 22.278385 12.827155

[29] 17.640852 16.770788

تابع نرمال - دستور زیر توزیع تابع نرمال  رابررسی می کند:

Norm   

به مثال زیر برای تابع نرمال است توجه کنید

Examples

(x1 <- cbind(1, 1:10))

norm(x1)

norm(x1, "I")

norm(x1, "M")

stopifnot(all.equal(norm(x1, "F"),

                    sqrt(sum(x1^2))))

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }

h9 <- hilbert(9)

## all 5 types of norm:

(nTyp <- eval(formals(base::norm)$type))

sapply(nTyp, norm, x = h9)

تابع پواسون- دستور زیر توزیع تابع پواسون رابررسی می کند:

Pois

به مثال زیر برای تابع پواسون است توجه کنید

Examples

require(graphics)

-log(dpois(0:7, lambda = 1) * gamma(1+ 0:7)) # == 1

Ni <- rpois(50, lambda = 4); table(factor(Ni, 0:max(Ni)))

1 - ppois(10*(15:25), lambda = 100)  # becomes 0 (cancellation)

    ppois(10*(15:25), lambda = 100, lower.tail = FALSE)  # no cancellation

par(mfrow = c(2, 1))

x <- seq(-0.01, 5, 0.01)

plot(x, ppois(x, 1), type = "s", ylab = "F(x)", main = "Poisson(1) CDF")

plot(x, pbinom(x, 100, 0.01), type = "s", ylab = "F(x)",

     main = "Binomial(100, 0.01) CDF")

آزمون t برای انجام مقایسه آماری بین میانگین دو گروه از آزمون t استفاده می شود . بسته به نوع مطالعه و استنباط های لازم آزمون t به صورت یک نمونه ای و یا دو نمونه ای  صورت می گیرد . در انواع یک نمونه ای  آزمون t نشان می دهد که آیا بین میانگین نمونه و یک مقدار فرضی تفاوت معنی دار وجود دارد یا خیر ؟ برای انجام انواع آزمون t در نرم افزار R از تابع t.test() وبرهان هایش استفاده می شود . برهان های این تابع به کاربر امکان می دهد که آزمون های t با فرضیات و شرایط گوناگون را اجرا برساند . تابعt-student- دستور زیر توزیع تابع t-student رابررسی می کند: t به مثال زیر برای تابع t است توجه کنید

Examples

a <- matrix(1:30, 5, 6)ta <- t(a) ##-- i.e.,  a[i, j] == ta[j, i] for all i,j :for(j in seq(ncol(a)))  if(! all(a[, j] == ta[j, ])) stop("wrong transpose")

weibull

Examples

x <- c(0, rlnorm(50))

all.equal(dweibull(x, shape = 1), dexp(x))

all.equal(pweibull(x, shape = 1, scale = pi), pexp(x, rate = 1/pi))

## Cumulative hazard H():

all.equal(pweibull(x, 2.5, pi, lower.tail = FALSE, log.p = TRUE),

          -(x/pi)^2.5, tolerance = 1e-15)

all.equal(qweibull(x/11, shape = 1, scale = pi), qexp(x/11, rate = 1/pi))

تابع یکنواخت- دستور زیر توزیع تابع یکنواخت رابررسی می کند: Unif به مثال زیر برای تابع یکنواخت است توجه کنید

Examples

u <- runif(20)

## The following relations always hold :

punif(u) == u

dunif(u) == 1

var(runif(10000))  #- ~ = 1/12 = .08333

تابع ویلکاکسون - دستور زیر توزیع تابع ویلکاکسون رابررسی می کند:

Wilcoxon

به مثال زیر برای تابع ویلکاکسون  است توجه کنید

Examples

require(graphics)

x <- -1L4*6 + 1)

fx <- dwilcox(x, 4, 6)

Fx <- pwilcox(x, 4, 6)

layout(rbind(1,2), widths = 1, heights = c(3,2))

plot(x, fx, type = “h”, col = “violet”,

     main =  “Probabilities (density) of Wilcoxon-Statist.(n=6, m=4)”)

plot(x, Fx, type = “s”, col = “blue”,

     main =  “Distribution of Wilcoxon-Statist.(n=6, m=4)”)

abline(h = 0:1, col = “gray20”, lty = 2)

layout(1) # set back

N <- 200

hist(U <- rwilcox(N, m = 4,n = 6), breaks = 0:25 – 1/2,

     border = “red”, col = “pink”, sub = paste(“N =”,N))

mtext(“N * f(x),  f() = true \”density\””, side = 3, col = “blue”)

 lines(x, N*fx, type = “h”, col = “blue”, lwd = 2)

points(x, N*fx, cex = 2)

## Better is a Quantile-Quantile Plot

qqplot(U, qw <- qwilcox((1:N – 1/2)/N, m = 4, n = 6),

       main = paste(“Q-Q-Plot of empirical and theoretical quantiles”,

                     “Wilcoxon Statistic,  (m=4, n=6)”, sep = “\n”))

n <- as.numeric(names(print(tU <- table(U))))

text(n+.2, n+.5, labels = tU, col = “red”)

با استفاده از تابع density() می توان تابع چگالی یک مجموعه داده را به دست آورد . برای رسم نمودار تابع چگالی احتمال نیز از تابع plot() برای خروجی density() استفاده می کنیم :

density(VADeaths)

Call:

        density.default(x = VADeaths)

Data: VADeaths (20 obs.);       Bandwidth 'bw' = 9.64

       x                 y           

 Min.   :-20.520   Min.   :2.748e-05 

 1st Qu.:  9.615   1st Qu.:1.690e-03 

 Median : 39.750   Median :7.787e-03 

 Mean   : 39.750   Mean   :8.286e-03 

 3rd Qu.: 69.885   3rd Qu.:1.333e-02 

 Max.   :100.020   Max.   :2.063e-02

به مثال زیر برای Density است توجه کنید

Examples

require(graphics)

plot(density(c(-20, rep(0,98), 20)), xlim = c(-4, 4))  # IQR = 0

# The Old Faithful geyser data

d <- density(faithful$eruptions, bw = "sj")

d

plot(d)

plot(d, type = "n")

polygon(d, col = "wheat")

## Missing values:

x <- xx <- faithful$eruptions

x[i.out <- sample(length(x), 10)] <- NA

doR <- density(x, bw = 0.15, na.rm = TRUE)

lines(doR, col = "blue")

points(xx[i.out], rep(0.01, 10))

## Weighted observations:

fe <- sort(faithful$eruptions) # has quite a few non-unique values

## use 'counts / n' as weights:

dw <- density(unique(fe), weights = table(fe)/length(fe), bw = d$bw)

utils::str(dw) ## smaller n: only 126, but identical estimate:

stopifnot(all.equal(d[1:3], dw[1:3]))

## simulation from a density() fit:

# a kernel density fit is an equally-weighted mixture.

fit <- density(xx)

N <- 1e6

x.new <- rnorm(N, sample(xx, size = N, replace = TRUE), fit$bw)

plot(fit)

lines(density(x.new), col = "blue")

(kernels <- eval(formals(density.default)$kernel))

## show the kernels in the R parametrization

plot (density(0, bw = 1), xlab = "",

      main = "R's density() kernels with bw = 1")

for(i in 2:length(kernels))

   lines(density(0, bw = 1, kernel =  kernels[i]), col = i)

legend(1.5,.4, legend = kernels, col = seq(kernels),

       lty = 1, cex = .8, y.intersp = 1)

## show the kernels in the S parametrization

plot(density(0, from = -1.2, to = 1.2, width = 2, kernel = "gaussian"),

     type = "l", ylim = c(0, 1), xlab = "",

     main = "R's density() kernels with width = 1")

for(i in 2:length(kernels))

   lines(density(0, width = 2, kernel =  kernels[i]), col = i)

legend(0.6, 1.0, legend = kernels, col = seq(kernels), lty = 1)

##-------- Semi-advanced theoretic from here on -------------

(RKs <- cbind(sapply(kernels,

                     function(k) density(kernel = k, give.Rkern = TRUE))))

100*round(RKs["epanechnikov",]/RKs, 4) ## Efficiencies

bw <- bw.SJ(precip) ## sensible automatic choice

plot(density(precip, bw = bw),

     main = "same sd bandwidths, 7 different kernels")

for(i in 2:length(kernels))

   lines(density(precip, bw = bw, kernel = kernels[i]), col = i)

## Bandwidth Adjustment for "Exactly Equivalent Kernels"

h.f <- sapply(kernels, function(k)density(kernel = k, give.Rkern = TRUE))

(h.f <- (h.f["gaussian"] / h.f)^ .2)

## -> 1, 1.01, .995, 1.007,... close to 1 => adjustment barely visible..

plot(density(precip, bw = bw),

     main = "equivalent bandwidths, 7 different kernels")

for(i in 2:length(kernels))

   lines(density(precip, bw = bw, adjust = h.f[i], kernel = kernels[i]),

         col = i)

legend(55, 0.035, legend = kernels, col = seq(kernels), lty = 1)

برای رسم نمودار تابع چگالی احتمال نیز از تابع plot() برای خروجی density() استفاده می کنیم : plot(density(VADeaths))

Examples

require(stats) # for lowess, rpois, rnormplot(cars)lines(lowess(cars)) plot(sin, -pi, 2*pi) # see ?plot.function ## Discrete Distribution Plot:plot(table(rpois(100, 5)), type = "h", col = "red", lwd = 10,     main = "rpois(100, lambda = 5)") ## Simple quantiles/ECDF, see ecdf() {library(stats)} for a better one:plot(x <- sort(rnorm(47)), type = "s", main = "plot(x, type = \s\")"")points(x, cex = .5, col = ""dark red"")