RPS_120118_g1015

2498 days ago by trybik

%html <a href="http://www.mimuw.edu.pl/~trybik/edu/0809/rps/r-skrypt.pdf">Skrypt do R</a> 
       
%html <p>Zadanie 1.6.5. Technikalia:</p> <ul> <li>stat. opisowa: mean, median, var (tab. 1.1)</li> <li>próbkowanie z rozkładów: rnorm, rt, rexp, runif (rodz. 1.4)</li> <li>iterowanie (sapply, replicate) i własne funkcje (rodz. 1.5)</li> </ul> 
       

Zadanie 1.6.5. Technikalia:

  • stat. opisowa: mean, median, var (tab. 1.1)
  • próbkowanie z rozkładów: rnorm, rt, rexp, runif (rodz. 1.4)
  • iterowanie (sapply, replicate) i własne funkcje (rodz. 1.5)
mean(1:10) 
       
[1] 5.5
v = rnorm(10) v f = function(x) 2*x sapply(v,f) 2*v 
       
 [1]  0.3304776 -0.3065795 -1.5273660 -1.8512027  1.1813446 
0.0996762 -0.4043618 -2.0172786
 [9] -0.3786503  0.6958622

 [1]  0.6609551 -0.6131590 -3.0547319 -3.7024053  2.3626893 
0.1993524 -0.8087236 -4.0345571
 [9] -0.7573007  1.3917243
 [1]  0.6609551 -0.6131590 -3.0547319 -3.7024053  2.3626893 
0.1993524 -0.8087236 -4.0345571
 [9] -0.7573007  1.3917243
sapply(rep(1,10), f) replicate(10, f(1)) 
       
 [1] 2 2 2 2 2 2 2 2 2 2
 [1] 2 2 2 2 2 2 2 2 2 2
zad5 = function(n=1000, m=100, mean=10){ rfuns = c( function(m) rnorm(m, mean=mean), function(m) mean+rt(m, 2), function(m) rexp(m, rate=1/mean) ) ret = sapply(rfuns, function(rfun) var(replicate(n,mean(rfun(m))))/var(replicate(n,median(rfun(m)))) ) names(ret) = c( paste("N(",mean,",1)", sep=""), paste("t(2)+",mean, sep=""), paste("Exp(1/",mean,")", sep="") ) ret } replicate(3,zad5()) 
       
               [,1]      [,2]      [,3]
N(10,1)   0.7211740 0.6325867 0.5541744
t(2)+10   6.2535260 5.7617489 4.9647831
Exp(1/10) 0.9425525 0.9557639 0.8921562
%html <p>Zadanie 1.6.10. Technikalia:</p> <ul> <li>f-cje std: max, seq (operator :)</li> <li>próbkowanie z d-punktowego jednostajnego: sample(d, n)</li> <li>ew. for i print</li> </ul> 
       

Zadanie 1.6.10. Technikalia:

  • f-cje std: max, seq (operator :)
  • próbkowanie z d-punktowego jednostajnego: sample(d, n)
  • ew. for i print
zad10 = function(n=100, d=n^2, m=10000) { .est1 = function(x) 2*mean(x)-1 .est2 = function(x) { l = length(x); max(x)*(l+1)/l } test.est = function(est.fun) { est.d.vec = replicate(m, est.fun(sample(d,n,replace=TRUE))) print(paste("Mean:", mean(est.d.vec), "vs.", d, "; Var:", var(est.d.vec))) } test.est(.est1) test.est(.est2) } zad10() zad10(n=10,d=10000) 
       
[1] "Mean: 9999.61877 vs. 10000 ; Var: 336584.306557183"
[1] "Mean: 10001.114738 vs. 10000 ; Var: 9537.98295108646"
[1] "Mean: 9988.97524 vs. 10000 ; Var: 3290215.7965986"
[1] "Mean: 9996.42545 vs. 10000 ; Var: 854651.900545352"
%html <p>Zadanie 1.6.7 + wizualizacja zbieżności do teoretycznej gęstości.</p> 
       

Zadanie 1.6.7 + wizualizacja zbieżności do teoretycznej gęstości.

zad7 = function(n, rate=1) { inv.cdf = function(x) -log(1-x)/rate hist.test = function(n) { s = inv.cdf(runif(n)) h = hist(s, freq=FALSE, breaks=10,main=paste(n,"samples"),ylim=c(0,rate),xlim=c(0,(6/rate))) rug(s, col="red") lines(density(s), col="red") curve(dexp(x, rate=rate), col="blue", add=TRUE) h } old.par = par(no.readonly=TRUE) par(cex=0.75) #par(mfrow=c(2,2)) #layout.show() #sapply(n.vec, hist.test) h = hist.test(n) par(old.par) h } hs = sapply(2^(7:10),function(n) zad7(n, rate=1/5)) png_output=dev.off() 
       
Warning messages:
1: In rug(s, col = "red") : some values will be clipped
2: In rug(s, col = "red") : some values will be clipped
3: In rug(s, col = "red") : some values will be clipped