R A MATLAB programcsomag alkalmazása valószínűségszámítási és statisztikai feladatokhoz
Tóth László
R MATLAB bevezető 1.
R MATLAB bevezető 2.
Mátrix létrehozása és invertálása: >> A=[1 2; 3 4] A = 1 2 3 4 >> inv(A) ans = -2.0000 1.5000
1.0000 -0.5000
Tömb feltöltése: >> s1=[1:10] s1 = 1 2
3
4
5
6
7
8
9
1
Binomiális eloszlás számítása 2. Jelölje X egy p valószínűségű A esemény bekövetkezéseinek számát. Ekkor X eloszlása binomiális: n k pk = P(X = k) = p (1 − p)n−k , k = 0, 1, . . . , n. k pk -t értékét adja meg binopdf(k,n,p) a függvény. BINOPDF - Binomial probability density function Az eloszlásfüggvény X F (x) = pk ,
−∞ < x < ∞
0≤k≤x
értékét binocdf(x,n,p) adja meg (Binomial cumulative distribution function)
Binomiális eloszlás számítása 1. Az eloszlásfüggvény F −1 (s) 0 ≤ s ≤ 1 inverze bioinv(s,n,p). A várható értéket és a szórásnégyzetet [kozep,variacia]=binostat(n,p) adja meg. Hasonló függvények léteznek a fontos diszkrét és folytonos eloszlásokra. geopdf(k,p), hygepdf(k,N,M,n), nbinpdf(k,r,p), poissonpdf(k,lampda), unidpdf(k,n) betapdf(x,a,b), unifpdf(x,a,b), evpdf(x,mu,sigma), lognpdf(x,mu,sigma), raylpdf(x,c), wblpdf(x,a,b)
1. Példa Egy gyárban egy minőségellenőr napi 400 terméket vizsgál meg. Tudjuk, hogy a termékek a (gyártástechnológiákól adódóan 1%-ban hibásak. Mi a valószínűsége, hogy az ellenőr egy adott nap egy hibásat sem talál? Binomiális eloszlás, n = 400, p = 0.01, k = 0 paraméterekkel! >> p=binopdf(0,400,0.01) p = 0.0180
2. Példa Hány hibás termék találása a legvalószínűbb? Megvizsgáljuk k összes szóbajöhető értékét, és ezek közül kiválasztjuk a legnagyobb valószínűségűt, és kiolvassuk a hozzájuk tartozó k-t. >> >> >> >>
hibasak=0:400; y = binopdf(hibasak,400,.01); [x,i]=max(y); db=hibasak(i)
db = 4
3. Példa Ha egy focicsapat egy idényben 90 mérkőzést játszik, és bármelyiken 50% a nyerési esély, mi a valószínűsége, hogy 55-nál több mérkőzést nyer? >> p2=1-binocdf(55,90,0.5) p2 = 0.0132 Az esetek túlnyomó többségében (98%) milyen intervallumban mozog a nyert mérkőzések száma? >> i=binoinv([0.01 0.99],90,0.5) i = 34
56
Részlet a dokumentációból 1.
Hipergeometriai eloszlás - 4. Példa Hipergeometriai eloszlás sűrűségfüggvénye: K M−K p(X = x) =
x
N−x M N
Kiszámítása: p=hygecdf(X,M,K,N) Van 80 külsőre egyforma merevlemezünk, melyek közül 25 hibás. Ha kiválasztunk 15 darabot, mi a valószínűsége, hogy 0..5 hibás lesz köztük? >> p=hygepdf(0:5,80,25,15) p = 0.0018 0.0164 0.0656 0.2342
0.1521
0.2281
Eloszlások ábrázolása - 5. példa Generáljunk 500 elemű mintát n = 20, p = 0.3 paraméterű binomiális eloszlásból (úgy, hogy annak a véletlen kísérlettel való definiálását használjuk)! Ábrázoljuk együtt az empirikus és az elméleti eloszlást! >> >> >> >>
X=rand(20,500); y=sum(X<0.3);[n,z]=hist(y,0:20); figure,bar(0:20,binopdf(0:20,20,0.3)), hold on; plot(z,n./500,’r*’);
>> >> >> >>
X=rand(20,500); y=sum(X<0.3);[n,z]=hist(y,0:20); figure,bar(0:20,binopdf(0:20,20,0.3)), hold on; plot(z,n./500,’r*’);
Normális eloszlás - 6. Példa A µ várható értékű, σ szórású normális eloszlású X sűrűségfüggvénye: 1 (x − µ)2 √ f (x) = exp − , −∞ < x < ∞ 2σ 2 σ 2π Ezen f (x) előállítása: normpdf(x,mu,sigma). Azonos várható értékű, különböző szórású normális eloszások előállítása: figure; % azonos varhato erteku normalis plot(-4:0.05:4,normpdf(-4:0.05:4,0,0.7),’r-’);hold on; plot(-4:0.05:4,normpdf(-4:0.05:4,0,1),’g--’);hold on; plot(-4:0.05:4,normpdf(-4:0.05:4,0,1.5),’k-’);
Részlet a dokumentációból 2.
figure; % azonos varhato erteku normalis plot(-4:0.05:4,normpdf(-4:0.05:4,0,0.7),’r-’);hold on; plot(-4:0.05:4,normpdf(-4:0.05:4,0,1),’g--’);hold on; plot(-4:0.05:4,normpdf(-4:0.05:4,0,1.5),’k-’);
Empirikus jellemzők számítása mintákból „Szokásos” jellemzők: mean,median,var,std, átlag, medián, várható érték, szórásnégyzet range,skewness,kurtosis,cov terjedelem, ferdeség, lapultság, kovariancia További jellemzők: iqr : empirikus kvartilisek távolsága (25% és 75%-os kvartilisek távolsága) prctile(X,p) : az X empirikus p-kvantilise, azaz azon mintaelem, amelynél a mintaelemek éppen p%-a kisebb-egyenlő. trimmean(X,p) : olyan számtani közép, melyben a minta alsó és felső (p/2)%-át nem vesszük figyelembe. corrcoef : empirikus korrelációs mátrix
7. Példa Generáljunk 299 √ elemű mintát a [−1, 1] intervallumon egyenletes, valamint a (0, 1/ 3) paraméterű normális eloszlásból! Ezek elméleti várható értékei és szórásai megegyeznek. Az empirikus jellemzők egy része csak kicsit fog eltérni, viszont a két minta függetlensége miatt kicsi korrelációt kell kapnunk x=unifrnd(-1,1,299,1); y=normrnd(0,1/(sqrt(3)),299,1); z=[x y]; kozep=mean(z), csonkakozep=trimmean(z,10), medi=median(z), mertanikozep=geomean(abs(z)), variancia=var(z), szoras=std(z), interkvar=iqr(z), terjedelem=range(z), m4=moment(z,4), lapultsag=kurtosis(z), ferdeseg=skewness(z), alkvantilis=prctile(z,5), felkvantilis=prctile(z,95), kovariancia=cov(z), korrelacio=corrcoef(z)
Számított empirikus jellemzők 1. kozep =
0.0830
-0.0067
medi =
0.1261
0.0402
csonkakozep =
0.0898
-0.0032
mertanikozep =
0.3604
0.3358
variancia =
0.3214
0.3349
szoras =
0.5670
0.5787
terjedelem =
1.9875
3.4518
interkvar =
0.9529
0.8080
Számított empirikus jellemzők 2. m4 =
0.1910
0.3092
lapultsag =
1.8609
2.7753
ferdeseg =
-0.1151
-0.1202
alkvantilis =
-0.8624
-0.9565
felkvantilis =
0.9372
0.8992
kovariancia =
0.3214 -0.0194
-0.0194 0.3349
1.0000 -0.0590
-0.0590 1.0000
korrelacio =
8. Példa - Statisztikai próbák 1. Legyen x 100 elemű minta standard normális eloszlásból, és y 100 elemű minta egyenletes eloszlásból úgy, hogy az elméleti várható értékek, illetve szórások megegyezzenek. Vizsgálódjunk 95%-os szinten! Ellenőrizzük Wilcoxon-féle rangösszeg próbával, hogy x és y azonos eloszlásból származik-e! >> x=normrnd(0,1,100,1); m=mean(x), s=std(x); >> y=unifrnd(-sqrt(3),sqrt(3),100,1); >> [pranksum,hranksum]=ranksum(x,y,0.05) pranksum = 0.9134 hranksum = 0 Láthatjuk, hogy a próba nem tudta megkülönböztetni a két eloszlást (a próba főleg eltolásra érzékeny).
9. Példa - Statisztikai próbák 2. Ellenőrizzük u próbával, hogy a H0 : E (x) = 0 hipotézis igaz-e? [hztest,sigztest,ciztest]=ztest(x,0,1,0.05) hztest = 0 sigztest = 0.2218 ciztest = -0.3182 0.0738 Az u próba a hipotézist elfogadta.
Előző példa folytatása Ábrázoljuk az u statisztika, valamint az alsó és felső kritikus értékek elhelyezkedését. u1=norminv(0.975,0,1); a=u1:0.02:4; figure; plot(-4:0.02:4, normpdf(-4:0.02:4,0,1),’k-’); hold on; fill([-u1 -a],[0 normpdf(-a)],’c’); fill([u1 a], [0 normpdf(a)],’c’); u=m*sqrt(50); plot(u,0:0.025:0.5,’k*’);
A besatírozott részek mutatják, hogy hová kéne esnie ahhoz az u statisztikának, hogy 0.05 szinten elvessük H0 -at.
10. Példa - Statisztikai próbák 3. Generáljunk egy új z mintát normális eloszlásból, és ellenőrizük kétmintás t-próbával, hogy a H0 : E (x) = E (z) igaz-e? [httest2,sigttest2,cittest2]=ttest2(x,z,0.05) httest2 = sigttest2 = cittest2 = -0.3371 0.2590
0 0.7963
Szórásanalízis - egyszeres ANOVA Három különböző takarmány hatását mérték 9,6, illetve 8 kísérleti állaton. Egyetlen tényezőként a takarmányt vizsgáljuk, aminek 3 szintje van. x=... [ 9.40 9.48 7.56 11.52 11.56 12.12 11.36 4.60 14.48 22.84 15.32 11.04 17.92 19.68 26.20 ... 17.35 16.36 15.88 14.28 18.60 19.32 14.20 17.52]; szam=... [1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3]; Az x vektorban tároljuk a tömegnövekedéseket A szam vektorban adjuk meg. hogy melyik mintalelem hányas számú takarmányhoz tartozik.
szor=[1 1 1]; for i=1:3 szor(i)=var(x(szam==i)); figure, hist(x(szam==i),3); figure, normplot(x(szam==i)); end; p=anova1(x, szam);
Eredményül a fenti szorzófelbontó táblát kapjuk A {Prob>F}=p=.0002 az jelenti, hogy elvetjük a takarmányok egyforma hatását.
Boxplot is keletkezik, a szintek eltérő hatása innen is látszik.
Főkomponens analízis Generáljunk 1000 elemű mintát kétdimenziós normális eloszlásból, és ábrázoljuk a két főkomponenst! m=1000; D2=[4 1; 1 2]; X2=mvnrnd(zeros(2,1),D2,m); [Fi2,Fk2,Saj2,Hott2]=princomp(X2); v=zeros(2,301); figure; axis([-8 8 -8 8]); hold on; plot(X2(:,1),X2(:,2),’c.’);hold on; for i=1:2 for j=1:2 v(j,:)=(0:Fi2(j,i)/100:3*Fi2(j,i))*sqrt(Saj2(i)); end plot(v(1,:),v(2,:),’b-’); hold on; end;
Főkomponens analízis folytatás Generáljuk adott szórásmátrixú 3-dimenziós normális eloszlásból 150 elemű mintát! Becsüljük a főkomponenseket a mintából! Az [Fi,Fk,Saj,Hott]=princomp(X) az X minta mátrix főirányait (Fi), főkomponenseit (Fk), sajátértékeit (Saj), és a Hotelling-féle T 2 statisztikát adja meg. A [pc,var,exlp]=pcacov(D) a megadott D elméleti szórásmátrix esetén számítja ki a főirányokat, a szórás mátrix sajátértékeit, és azt, hogy az egyes főkomponensek a szórás hány százalékát magyarázzák. A k főkomponens levonása utáni maradékot a maradek=pcares(X,k) adja meg.
>> >> >> >> >> >> >>
n=150; a1=[1 -1 1]’; a1=a1./sqrt(3); a2=[1 0 -1]’; a2=a2./sqrt(2); a3=[1 2 1]’; a3=a3./sqrt(6); A=[a1 a2 a3]; D=A*diag([4 1 0.1])*A’; X=mvnrnd(zeros(3,1),D,n); [Fi,Fk,Saj,Hott]=princomp(X); maradek=pcares(X,2); [pc,var,expl]=pcacov(D);
Az eredményül kapott főirányok és sajátértékek: >> Fi = -0.5863 0.5800 -0.5656 >> Saj = 4.0298 0.9814 0.1035
0.6867 -0.0145 -0.7268
0.4297 0.8145 0.3898
pc = -0.5774 0.5774 -0.5774 var = 4.0000 1.0000 0.1000 expl = 78.4314 19.6078 1.9608
-0.7071 -0.0000 0.7071
0.4082 0.8165 0.4082
Irodalom MATLAB - Frissített kiadás Numerikus módszerek, grafika, statisztika, eszköztárak Stoyan Gisbert (szerk.) TypoTeX 2005, Alkalmazott matematika sorozat Statistical Analysis Techniques in Particle Physics: Fits, Density Estimation, and Supervised Learning Narsky / Porter Wiley-VCH Verlag GmbH, 2013 Probability, Statistics, and Random Processes for Engineers, 4e Stark / Woods Pearson Education Inc, 2012 Nonparametric Statistics with Applications to Science and Engineering Kvam / Vidakovic John Wiley & Sons, Inc. 2007
Köszönöm a figyelmet!
