[R] montecarlo simulation

[Vito Ricci]

Ciao Francesca,

here are some examples or R code using montecarlo
simulation concernig parameters estimates of sever
random variables:


Estimate of parameter theta of a uniform variable with
ML method
Stima del parametro theta di una v.a. unif(0,theta)
con metodo
dei momenti e confronto con la stima di massima
verosimiglianza (M.V.)


theta<-5
n<-50
x<-runif(n,0,theta)
est.theta<-2*mean(x) ## stima con metodo dei momenti
var.theta<-(est.theta^3)/(3*n) ## stima della varianza
di theta
mlest.theta<-max(x) ## stima di theta con il metodo
della M.V.
acc1<-(theta-est.theta)/theta
acc2<-(theta-mlest.theta)/theta
est.theta-mlest.theta

Simulazione di una v.a. uniforme (0,t), stima di theta
con metodo dei
momenti e della M.V. e stampa dei relativi grafici

t<-2.5 ## valore parametro
n<-20  ## dimensione campionaria
k<-200 ## numero di campioni da simulare
X<-matrix(runif(n*k,0,t),nrow=n)
est.tmom<-2*apply(X,2,mean)
est.tml<-apply(X,2,max)
hist(est.tmom,nclass=10)
hist(est.tml, nclass=10,add=T,col="red")
var.x<-apply(X,2,var) ## varianza della v.a.
var.mom<-var(est.tmom)## varianza della stima com
metodo dei momenti
var.ml<-var(est.tml)  ## varianza della stima com
metodo della M.V.

medie.stime<-c(mean(est.tml),mean(est.tmom))
var.stime<-c(var.ml,var.mom)
tabella<-data.frame(medie.stime,var.stime,metodo<-c("Stima
ML","Stima Momenti"))
names(tabella)[3]<-"Metodo"
tabella

Weibull simulation

n<-50
k<-500
X<-matrix(rweibull(n*k,shape=10,scale=5),nrow=n)
mediana<-apply(X,2,median)
avg.med<-mean(mediana)
var.med<-var(mediana)
hist(mediana,freq=T)
media<-apply(X,2,mean)
avg.media<-mean(media)
var.media<-var(media)
hist(media,freq=T)
tabella<-data.frame(c(avg.med,avg.media),c(var.med,var.media))
names(tabella)<-c("statistica","varianza")
qqnorm(mediana)
qqnorm(media)
shapiro.test(mediana)
shapiro.test(media)
qqplot(mediana,media)


Beta pdf simulation
Simulazione di una variabile beta

N<-1000 ##numero osservazioni totali
n<-50 ##numero di elementi per campione
k<-N/n ##numero dei campioni
for (i in 1:n)
{
X<-matrix(rbeta(N,10,20),nrow=n)
}

n<-20
k<-200
X<-matrix(rbeta(n*k,10,20),nrow=n)

Estimate of parameters of a beta pdf with ML method
using nlm()
Stima dei parametri di una v.a. Beta con metodo della
massima verosimiglianza 
usando la funzione nlm() per un campione di ampiezza
n=100 e dei relativi errori standard

 x<-rbeta(100,10,20)
 f<-function(p)
100*log(beta(p[1],p[2]))+(1-p[1])*sum(log(x))+(1-p[2])*sum(log(1-x))
 out<-nlm(f,p=c(20,20),hessian=TRUE)


Estimate of parameters of a Gamma pdf with ML method
using nlm().
Stima di massima verosimiglianza dei parametri di una
v.a. Gamma da cui è stato estratto un campione di
ampiezza n=200 con l'uso della funzione nlm()

x<-rgamma(200,10,5)
f<-function(p)
(-200*p[1]*log(p[2])+200*log(gamma(p[1]))+(1-p[1])*sum(log(x))+p[2]*sum(x))
out<-nlm(f,p=c(15,10),hessian=TRUE)
out

Best
Vito



Hy!
I would like to know how run a montecarlo simulation
with R.
Thank you!!!!
Francesca Matalucci

=====
Diventare costruttori di soluzioni

Visitate il portale http://www.modugno.it/
e in particolare la sezione su Palese http://www.modugno.it/archivio/cat_palese.shtml


		
___________________________________

http://it.seriea.fantasysports.yahoo.com/