Hi:
A univariate t-distribution with 1 df is equivalent to the standard Cauchy
distribution, for
which it is well established that the population mean does not exist. You're
basically
simulating a vector version of the Cauchy distribution (an IID 8-dimensional
version),
so the same problem is likely to arise. This is exhibited in the tail
behavior of the
(simulated) distribution of means.
It's a little more convenient to do this in terms of arrays rather than
lists:
library(mvtnorm)
# Generate 10 multivariate t samples in 8D and replicate the process 1000
times
# result is a 3D array of size 10 x 8 x 1000.
rs <- array(replicate(1000, rmvt(10, sigma = diag(8), df = 1)), c(10, 8,
1000))
# Average over the 10 samples in each replicate, yielding an 8 x 1000 matrix
-
# i.e., average over the first dimension of the array
rsm <- apply(rs, c(2, 3), mean)
# Transpose the matrix, convert it to a data frame, and melt it using the
# reshape package; this generates an 8000 x 2 data frame, where the first
# variable is an indicator for each component of the vector (V1 - V8) and
# the value represents the sample values.
library(reshape)
rr <- melt(as.data.frame(t(rsm)))
# Use the histogram function in lattice to generate the histograms for each
# component of the sample vector:
library(lattice)
histogram(~ value | variable, data = rr, xlab = 'x', type = 'density')
Since your original post averaged over the entire matrix in each replicate,
let's do that - just because...
rmns <- apply(rs, 3, mean)
hist(rmns, breaks = 20, probability = TRUE)
#or in lattice,
histogram(~ rmns, breaks = 20, density = TRUE)
[Hint: If you use the median rather than the mean as an estimator, its
sampling distribution will behave somewhat closer to what you
may expect...]
HTH,
Dennis
On Sun, May 9, 2010 at 10:37 PM, Shant Ch <sha1...@yahoo.com> wrote:
> Hello,
>
> I want to draw a histogram of the mean of sample observations drawn from
> multivariate t distribution. I am getting the following error corresponding
> to the code I used. Though I am getting the graph, but I am curious to know
> the warning message.
>
> Warning messages:
> 1: In if (freq) x$counts else { :
> the condition has length > 1 and only the first element will be used
> 2: In if (!freq) "Density" else "Frequency" :
> the condition has length > 1 and only the first element will be used
>
> -----------------------------------------------------------------------------------
> library(mvtnorm);
> s12<-c(1:1000);
> var21<-lapply(s12,function(x){
> rs<-rmvt(10, sigma=diag(8), df=1);
> my<-mean(rs);
> sy<-sqrt(var(rs))
> return(cbind(my,sy))
> });
> data1<-do.call(rbind,var21);
> dataMat<-data.frame(data1);
> W<-dataMat$my;
> hist(W,breaks=20,probability=T)
>
> -----------------------------------------------------------------------------------
>
> Thanks.
>
> Shant
>
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
