MSA 2013

Literatura

Rao, C. R.: Lineární metody statstické indukce a jejich aplikace

Haerdle, W., Simar, L.: Applied Multivariate Statistical Analysis

Mardia, K. V., et al.: Multivariate Analysis

Christensen, R. Linear Models for Multivariate, Time Series, and Spatial Data

Haerdle, W., Hlávka, Z.: Multivariate Statistics: Exercises and Solutions

Mnohorozměrná rozdělení

Clear[μμ,ΣΣ,xx,μ,Σ,σ,x]
μμ[p_]:=Array[(Subscript[μ, ##]) &,p]
ΣΣ[p_]:=Array[(Subscript[σ, ##]) &,{p,p}]*Boole[Array[LessEqual,{p,p}]]+
Transpose[Array[(Subscript[σ, ##]) &,{p,p}]*Boole[Array[Less,{p,p}]]]
xx[p_]:=Array[(Subscript[x, ##]) &,p]
xmatrix[p_]:=Array[(Subscript[x, ##]) &,{p,p}]*Boole[Array[LessEqual,{p,p}]]+
Transpose[Array[(Subscript[x, ##]) &,{p,p}]*Boole[Array[Less,{p,p}]]]

Boole[Array[Greater,{5,5}]]//MatrixForm

μμ[3]/.μμ[3]->{1,2,3}
(μμ[3]*μμ[3])/.{Thread[μμ[3]->{1,2,3}]}
(μμ[3]. μμ[3])/.{Thread[μμ[3]->{1,2,3}]}

μμ[3]

ΣΣ[4]//MatrixForm

?*Multi*

?*Multivar*

Normální rozdělení

μμ[3]

MultinormalDistribution

(xx[p]/.p:>  4)//Evaluate

PDF[MultinormalDistribution[μμ[3],ΣΣ[3]],xx[3]]//Simplify//TraditionalForm

PDF[MultinormalDistribution[{(Subscript[μ, 1]),(Subscript[μ, 2])},{{((Subscript[σ, 1])^2),ρ (Subscript[σ, 1]) (Subscript[σ, 2])},{ρ (Subscript[σ, 1]) (Subscript[σ, 2]),((Subscript[σ, 2])^2)}}],{(Subscript[x, 1]),(Subscript[x, 2])}]//TraditionalForm

Manipulate[Plot3D[PDF[MultinormalDistribution[{-1,1},{{1,ρ 2},{ρ 2,4}}],{x,y}],{x,-3,3},{y,-6,6},PlotRange->All,PlotPoints->50],{{ρ,0},-1,1},SaveDefinitions->True]

Wishartovo rozdělení

Needs["MultivariateStatistics`"]

PDF[WishartDistribution[ΣΣ[3],10],xmatrix[3]]

Mean[WishartDistribution[ΣΣ[3],10]]//MatrixForm

Variance[WishartDistribution[ΣΣ[3],10]]//MatrixForm

Hotellingovo (T^2)

?*Distr*

Mean[FRatioDistribution[m,n]]//TraditionalForm

Mean[HotellingTSquareDistribution[m,n]]//TraditionalForm

Table[Mean[HotellingTSquareDistribution[10,m]],{m,{15,20,30}}]

Plot[Evaluate@Table[PDF[HotellingTSquareDistribution[10,m],x],{m,{15,20,30}}],{x,0,40},Filling->Axis]

Manipulate[Plot[PDF[HotellingTSquareDistribution[m,n],x]//Evaluate,{x,0,40},PlotRange->All],{m,3,10},{n,15,30},SaveDefinitions->True]

Mnohorozměrné t-rozdělení

Normální rozdělení 3

Kopuly

Testy dobré shody

Hyperlink[Button["TimeSeriesProcesses"],"paclet:guide/TimeSeriesProcesses"]

Row[{Hyperlink[Button["DistributionFitTest"],"paclet:ref/DistributionFitTest"]
Hyperlink[Button["KuiperTest"],"paclet:ref/KuiperTest"]}]

?*Fit*

Analýza hlavních komponent (Principal Components)

I když to není v návodu (v systému Mathematica) přímo napsáno, ze zkušenosti se zdá, že Eigenvalues (Subscript[λ, 1])≥(Subscript[λ, 2])≥...≥(Subscript[λ, n]) symetrické positivně semidefinitní matice A (n×n) jsou setříděné sestupně podle velikosti a jim odpovídající vlastní vektory (Subscript[x, 1]),...,(Subscript[x, n]) jsou ortonormální. Pak je odpovídající spektrální rozklad

Spektrální rozklad varianční matice Σ (p×p):

Σ=(∑_(i=1)^p)(Subscript[λ, i])(Subscript[x, i])(Subscript[x, i]^T).

Data are in the form 𝒳=((((Subscript[X, 1]),...,(Subscript[X, n])))^T), in our case,  n=10, p=4. The built-in function   PrincipalComponents[data]  returns the transform of the original data to the principal components,

((Subscript[X, 1]),...,(Subscript[X, n])) ↦ ((Subscript[Y, 1]),...,(Subscript[Y, n])),

where

(Subscript[Y, i])=((((Subscript[x, 1]^T)(Subscript[X, i]),...,(Subscript[x, p]^T)(Subscript[X, i])))^T),..., i=1,...,n,

so that the principal components are calculated for every observation.

(data = {{13.2, 200, 58, 21.2}, {10, 263, 48, 44.5}, {8.1, 294, 80, 31}, {8.8, 190, 50, 19.5}, {9, 276, 91, 40.6}, {7.9, 204, 78, 38.7}, {3.3, 110, 77, 11.1}, {5.9, 238, 72, 15.8}, {15.4, 335, 80, 31.9}, {17.4, 211, 60, 25.8}})//TableForm
Dimensions[data]

Here we give the basic characteristics:

pc=PrincipalComponents[data];
cov=Covariance[data];
{Variance[pc],Total[Variance[pc]],Tr[cov]}
Eigenvalues[cov]
Eigenvalues[Correlation[data]]

Závisí na jednotkách:

data2=data;
data2[[All,1]]=1000 data2[[All,1]]

Variance[PrincipalComponents[data2]]

Eigenvectors[cov]//MatrixForm
Eigenvectors[Covariance[data2]]//MatrixForm

Eigenvectors[Correlation[data]]//MatrixForm
Eigenvectors[Correlation[data2]]//MatrixForm

Tr[cov]

{Dimensions[pc],Dimensions[cov]}

Correlation[pc]//Chop//MatrixForm

Comparison of PC’s calculating from the covariance matrix and correlation matrix. How the PC’s change?

dataStandStandardize=Standardize[data];
Mean[dataStandStandardize]
Variance[dataStandStandardize]

Diskriminační analýza

Lososi

SetDirectory[NotebookDirectory[]]

$ImportFormats

lososCelySoubor=Import["salmon.txt","Table"];
losos=Most[Rest[lososCelySoubor]];
lososAlaska=Select[losos,#[[2]]=="Alaska"&]
lososCanada=Select[losos,#[[2]]=="Canada"&]

lososAlaska[[1]]

{(covarianceAlaska=Covariance[Map[Take[#,{3,4}]&,lososAlaska]])//MatrixForm,
(covarianceCanada=Covariance[Map[Take[#,{3,4}]&,lososCanada]])//MatrixForm}//N

Možné též stáhnout přímo:

(i2=Import["https://www.scss.tcd.ie/~houldinb/Index/MVA_files/salmon.txt","Table"])//FullForm

Test hypotézy o rovnosti dvou variančních matic

Test poměrem věrohodností (Likelihood Ratio Test, LRT)

Covariance[{{1,2,3},{10,12,14}}]

Clear[EqualCovarinacesLRT];
EqualCovarinacesLRT[data1_,data2_]:=Module[{n1,n2,n,s1,s2,spooled,p,statistic,degreesoffreedom},
n1=Length[data1];n2=Length[data2];n=n1+n2;
p=Length[data1//Transpose];
s1=Covariance[data1//N];s2=Covariance[data2//N];
spooled=(1./n)(n1 s1+n2 s2);
statistic=n Log[Det[spooled]]-(n1 Log[Det[s1]]+n2 Log[Det[s2]]);
degreesoffreedom=(1/2)p(p+1);
Text@Grid[{{"Covariance matrices:"},{spooled//MatrixForm,s1//MatrixForm,s2//MatrixForm},{"\n"},{statistic,p},{"P-value ",SurvivalFunction[ChiSquareDistribution[degreesoffreedom],statistic]}}]
  	]

EqualCovarinacesLRT[Map[Take[#,{3,4}]&,lososAlaska],Map[Take[#,{3,4}]&,lososCanada]]

Faktorová analýza

Shluková analýza

Hierarchické metody

Help ⟶ Documentation Center ⟶ HierarchicalClustering/tutorial/HierarchicalClustering

Hyperlink[Button["HierarchicalClustering"],"paclet:HierarchicalClustering/tutorial/HierarchicalClustering"]

Needs["HierarchicalClustering`"]

data={{0.7253055971477931`,-1.1302536982240747`},{0.5810146986977534`,1.4463710601068265`},{0.8460694069050189`,2.0329569918374406`},{0.2169129449840166`,-0.33368074105304846`},{-0.02675560596423776`,1.1506670168024957`},{-0.38070042430357054`,-0.24750818092033833`},{2.4423867332670177`,1.9380008134953024`},{-0.3152439519516439`,-0.2760546508444408`},{1.1616304061543707`,0.794607673348773`},{-0.25171211458875775`,2.0210234868823256`}};

(dc=DistanceMatrix[data])//MatrixForm

DendrogramPlot[data,LeafLabels->Automatic,Orientation->Left]

GraphicsRow[{DendrogramPlot[data,LeafLabels->Automatic,Orientation->Left,Linkage->"Single"],DendrogramPlot[data,LeafLabels->Automatic,Orientation->Left,Linkage->"Complete",PlotStyle->{Orange}]}]

Nehierarchické metody

Stojí na superfunkci FindClusters

Hyperlink[Button["Nehierarchické shlukování"],"paclet:tutorial/PartitioningDataIntoClusters"]

Demonstrace

Help ⟶ Documentation Center ⟶ tutorial/PartitioningDataIntoClusters

http://demonstrations.wolfram.com/ClusterAnalysis/

FindClusters

FindClusters

Generovani 4 vyberu o rozsahu 100 z dvourozmerneho normalniho rozdeleni:

GaussianRandomData[n_Integer,p_,sigma_]:=Table[p+{Re[#],Im[#]}&[RandomReal[NormalDistribution[0,sigma]](E^(I RandomReal[{0,2π}]))],{n}];
datapairs=BlockRandom[
SeedRandom[1234];
Join[GaussianRandomData[100,{2,1},.3],GaussianRandomData[100,{1,1.5},.2],
GaussianRandomData[100,{1,1.1},.4],GaussianRandomData[100,{1.75,1.75},0.1]]];

c1=FindClusters[datapairs];
p1=ListPlot[c1,PlotStyle->{PointSize[0.02]}]

c1rot=FindClusters[RotateLeft[datapairs,50]];
p1rot=ListPlot[c1rot,PlotStyle->{PointSize[0.012]}]

Show[p1,p1rot]

c2=FindClusters[datapairs,Method->"Agglomerate"];
ListPlot[c2,PlotStyle->{PointSize[0.02]}]

c3=FindClusters[datapairs,DistanceFunction->CanberraDistance];
ListPlot[c3,PlotStyle->{PointSize[0.02]}]

c4=FindClusters[datapairs,5,Method->"Agglomerate"];
ListPlot[c4,PlotStyle->{PointSize[0.02]}]

Zadavani dat pomoci transformaci

Data {1,2,10,...}

FindClusters[{1->a,2->b,10->c,12->d,3->e,1->f,13->g,25->h}]

Pozorovani : {Jmeno, Prijmeni, Vaha, Mira}

datarecords={{"Joe", "Smith",158,64.4},{"Mary", "Davis",137,64.4},{"Bob", "Lewis",141,62.8},{"John", "Thompson",235,71.1},{"Lewis", "Black",225,71.4},{"Sally", "Jones",168,62.},{"Tom", "Smith",243,70.9},{"Jane", "Doe",225,71.4}};

Shlukovani pouze na zaklade vahy a vysky

FindClusters[Drop[datarecords,None,{1,2}]->datarecords]//MatrixForm

Mnohorozměrné škálování (Multidimensional scaling)

http://demonstrations.wolfram.com/MultidimensionalScaling/

JH: Discriminant Analysis and Multidimensional Scaling:
http://library.wolfram.com/infocenter/Conferences/5777/

ANOVA a porovnávání průměrů

<<ANOVA`

onewaydata={{1,7.0},{1,5.3},{1,5.9},{1,6.6},{1,4.9},{2,4.4},{2,6.8},{2,7.7},{2,8.3},{2,6.6},{3,8.1},{3,10.4},{3,8.0},{3,6.8},{3,9.2},{4,5.7},{4,3.9},{4,6.2},{4,5.5},{4,6.2}};

Map[Function[j,{j,Select[onewaydata,#[[1]]==j&][[All,2]]//Mean}],Union[onewaydata[[All,1]]]]
Sort[%,#1[[2]]<#2[[2]]&]

ANOVA[onewaydata,PostTests->Tukey]//TableForm

Skupiny prumeru, ktere se vyznamne nelisi, podtrhneme:

_______________________________
                                              ____________________

ANOVA[onewaydata,PostTests->{Duncan,Tukey,Bonferroni,StudentNewmanKeuls},SignificanceLevel->.006,CellMeans->False]//TableForm

Maximální věrohodnost

Likelihood
LogLikelihood

Testy dobré shody (tests of goodness of fit)

DistributionFitTest

({{      , "DistanceToBoundary", uniformity, based on distance to uniform boundaries}, {      , "MardiaCombined", normality, combined Mardia skewness and kurtosis}, {      , "MardiaKurtosis", normality, based on multivariate kurtosis}, {      , "MardiaSkewness", normality, based on multivariate skewness}, {      , "SzekelyEnergy", data, based on Newton's potential energy}})

SeedRandom[13];
data=RandomVariate[NormalDistribution[],2500];
ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]

SeedRandom[13];
data=RandomVariate[CauchyDistribution[0,1],2500];
ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]
{Mean[data],Median[data],StandardDeviation[data]}

{Min[data],Max[data]}

Histogram[data,{-10,10,1},"PDF",Epilog->First@Plot[PDF[NormalDistribution[0,1],x],{x,-4,4},PlotStyle->Red]]

Histogram

SeedRandom[13];data=RandomVariate[StudentTDistribution[100],2500];
ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]

SeedRandom[13];data=RandomVariate[StudentTDistribution[50],2500];
ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]