Bu haftaki konumuz veri üretimi.
set.seed(41)
x <- rnorm(100,100,20)
mean(x)
## [1] 103.8063
sd(x)
## [1] 20.56458
set.seed(410)
x<- rnorm(100,100,20)
mean(x)
## [1] 99.62622
sd(x)
## [1] 19.06762
n_k <-10
b_k <-50
tekrar <- 100
kucuk_orn <- numeric(tekrar)
buyuk_orn <- numeric(tekrar)
for(i in 1:tekrar){
kucuk_orn[i] <- mean(rnorm(n=n_k, mean=35, sd=15))
buyuk_orn[i] <- mean(rnorm(n=b_k, mean=35, sd=15))
}
hist(kucuk_orn)
hist(buyuk_orn)
dir.create("Simulasyon_k")
## Warning in dir.create("Simulasyon_k"): 'Simulasyon_k' zaten var
dir.create("Simulasyon_b")
## Warning in dir.create("Simulasyon_b"): 'Simulasyon_b' zaten var
for(i in 1:5) {
write.table(rnorm(n=n_k, mean=35, sd=15),
file=paste0("Simulasyon_k/simulasyon_", i, ".txt", sep=""))
write.table(rnorm(n=b_k, mean=35, sd=15),
file=paste0("Simulasyon_b/simulasyon_", i, ".txt", sep=""))
}
kucuk_orn_ort <- list()
buyuk_orn_ort <- list()
#for(i in 1:50){
#buyuk_orn_ort[[i]] <-read.table(file=paste0("Simulasyon_b/simulasyon_", i, ".txt",sep=""))[,1]
#kucuk_orn_ort[[i]] <-read.table(file=paste0("Simulasyon_k/simulasyon_", i, ".txt",sep=""))[,1]
#}
for döngüsünde anlayamadığım bir şekilde hata aldığım için diyezledim.
bunları da kullanabilirsin:
lapply(buyuk_orn_ort, mean)
sapply(buyuk_orn_ort, mean)
sapply(kucuk_orn_ort, mean)
library(psych)
library(irtoys)
## Warning: package 'irtoys' was built under R version 4.5.3
## Zorunlu paket yükleniyor: sm
## Warning: package 'sm' was built under R version 4.5.3
## Package 'sm', version 2.2-6.0: type help(sm) for summary information
## Zorunlu paket yükleniyor: ltm
## Warning: package 'ltm' was built under R version 4.5.3
## Zorunlu paket yükleniyor: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:sm':
##
## muscle
## Zorunlu paket yükleniyor: msm
## Warning: package 'msm' was built under R version 4.5.3
## Zorunlu paket yükleniyor: polycor
##
## Attaching package: 'polycor'
## The following object is masked from 'package:psych':
##
## polyserial
##
## Attaching package: 'ltm'
## The following object is masked from 'package:psych':
##
## factor.scores
##
## Attaching package: 'irtoys'
## The following object is masked from 'package:psych':
##
## sim
library(mirt)
## Zorunlu paket yükleniyor: stats4
## Zorunlu paket yükleniyor: lattice
##
## Attaching package: 'mirt'
## The following object is masked from 'package:ltm':
##
## Science
hangi model, kaç madde, madde ve yetenek par
3pl 8 maddelik
set.seed(41)
madde <- 8
maddepar <- cbind(
rnorm(madde, mean=0, sd=0.75),
rnorm(madde, mean=0.30, sd=0.50),
rnorm(madde, mean=0.16, sd=0.05)
)
set.seed(41)
birey <- 1000
yetenek <- rnorm(birey, 0, 1)
yetenek[1:10]
## [1] -0.7943683 0.1972575 1.0017043 1.2888254 0.9057534 0.4936675
## [7] 0.5992858 -1.5796070 1.0006207 2.1880077
veri <- irtoys::sim(ip=maddepar, x=yetenek)
head(veri)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 1 1 0 0 0 1 1
## [2,] 1 0 1 1 1 1 0 0
## [3,] 1 1 1 1 0 1 0 1
## [4,] 1 1 0 1 0 1 1 0
## [5,] 0 1 0 1 1 1 1 1
## [6,] 1 0 1 0 0 1 0 0
colnames(veri) <- paste0("madde", 1:madde)
head(veri)
## madde1 madde2 madde3 madde4 madde5 madde6 madde7 madde8
## [1,] 0 1 1 0 0 0 1 1
## [2,] 1 0 1 1 1 1 0 0
## [3,] 1 1 1 1 0 1 0 1
## [4,] 1 1 0 1 0 1 1 0
## [5,] 0 1 0 1 1 1 1 1
## [6,] 1 0 1 0 0 1 0 0
ip=maddepar
x=yetenek
sim
## function (ip, x = NULL)
## {
## if (is.list(ip))
## ip = ip$est
## i = irf(ip = ip, x = x)
## d = dim(i$f)
## u = runif(d[1] * d[2])
## dim(u) = d
## return(ifelse(i$f > u, 1, 0))
## }
## <bytecode: 0x000001e816680690>
## <environment: namespace:irtoys>
#if (is.list(ip))
#ip=ip$est
i=irf(ip=ip, x=x)
head(i$x)
## [1] -0.7943683 0.1972575 1.0017043 1.2888254 0.9057534 0.4936675
head(i$f)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.7749947 0.5380812 0.5095941 0.5029471 0.3356915 0.5706020 0.4451587
## [2,] 0.6682795 0.5668132 0.6623887 0.6702723 0.4641372 0.6372562 0.5432340
## [3,] 0.5724039 0.5904107 0.7737769 0.7992191 0.5835847 0.6913506 0.6237813
## [4,] 0.5382811 0.5988614 0.8074882 0.8369551 0.6262072 0.7101708 0.6516678
## [5,] 0.5839099 0.5875885 0.7617225 0.7854700 0.5691916 0.6849854 0.6143199
## [6,] 0.6333260 0.5754880 0.7059016 0.7208425 0.5074758 0.6573343 0.5731340
## [,8]
## [1,] 0.8453039
## [2,] 0.6450350
## [3,] 0.4492294
## [4,] 0.3881618
## [5,] 0.4712631
## [6,] 0.5716393
d=dim(i$f)
u=runif(d[1]*d[2])
dim(u)=d
head(ifelse(i$f >u, 1, 0))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 0 0 1 1 1
## [2,] 1 1 0 0 0 1 1 1
## [3,] 1 1 1 1 1 1 0 0
## [4,] 1 0 1 0 1 1 1 1
## [5,] 1 1 1 1 1 1 1 1
## [6,] 1 0 0 1 1 1 1 0
madde parametresi üretsin, yetenek parametresi üretsin, madde cevapları üretsin madde par, yetenek par, madde cevap
veri_uretimi <- function (maddesay, bireysay, seed){
# seet fonksiyonu
set.seed(seed)
# madde par üretimi
print("madde parametreleri uretildi")
maddepar <- cbind(
rnorm(madde, mean=1.13, sd=0.25)*1.702,
rnorm(madde, mean=0.30, sd=0.50)*1.702,
rnorm(madde, mean=0.16, sd=0.05)*1.702
)
# yetenek par
yetenek <- rnorm(bireysay,0,1)
print("yetenek parametreleri uretildi")
}
#3pl modele göre veri üretimi
cevaplar <-irtoys::sim(maddepar, yetenek)
colnames(cevaplar) <- paste0("madde",1: madde)
print("bilesenler birlestirildi")
## [1] "bilesenler birlestirildi"
veri <- list(maddepar=maddepar,
yetenek=yetenek,
cevaplar=cevaplar)
#return(veri)
veri_1 <-veri_uretimi(maddesay=8, bireysay= 1000, seed=66)
## [1] "madde parametreleri uretildi"
## [1] "yetenek parametreleri uretildi"
kestirilen par isimli fonksiyon; veri_üretimi fonksiyonunun cevaplar bileşeni girdi olarak
kestirilen_par <- function(veri, par=3){
if(par==3){
model <- mirt::mirt(veri,
1,
itemtype="3PL",
verbose= FALSE)
} else if(par==2){
model <- mirt::mirt(veri,
1,
itemtype="2PL",
verbose= FALSE)
}else{
model <- mirt::mirt(veri,
1,
itemtype="1PL",
verbose= FALSE)
}
kestirim <- mirt::coef(model, simplify=TRUE, IRTpars=TRUE)$item[,1:3]
kestirim
}
### kestirilen_par(veri_1$cevaplar)
### veri_1$maddepar
### kestirilen madde parametreleri ile gerçek parametreler arasındaki farka göre
### hata RMSEA^2, BIAS^2 (sistematik hatayı temsil eder), SE^2 (standart hata-kestirilen parametre değerinin standart sapmasıdır)
### RMSEA^2 = BIAS^2 + SE^2
hata <- function(kestirilen, gercek){
result <- data.frame(parametreler =c("a","b","c"),
rmse=sapply(1:3, function(i)sqrt(mean((kestirilen[,i]- gercek[,i])^2))),
bias=sapply(1:3, function(i)mean(kestirilen[,i]-gercek[,i])),
se= sapply(1:3, function(i) sd(kestirilen[,i])))
result
}
### hata(kestirilen_par(veri_1$cevaplar), veri_1$maddepar)
Yeterli sonuçlara, minimum hatalara “1 tekrar” bizi götürmüyor genellikle
veri_uretimi <- function(maddesay, bireysay, seed=NULL) {
if(!is.null(seed)){
set.seed(seed)
}else{
seed <- sample.int(10000,1)
set.seed(seed)
}
#seet fonksiyonu
cat("atanan seed", seed, " n")
#madde par üretimi
print("madde parametreleri uretildi")
maddepar <- cbind(
rlnorm(maddesay, meanlog=0, sdlog=0.5)*1.702,
rlnorm(maddesay, meanlog=0.30, sdlog=0.50)*1.702,
rlnorm(maddesay, meanlog=0.16, sdlog=0.50)
)
#yetenek par
yetenek <- rnorm(bireysay,0,1)
print("yetenek parametreleri uretildi")
#3 pl modele göre veri seti
cevaplar <- irtoys::sim(maddepar, yetenek)
colnames(cevaplar) <- paste0("madde", 1:maddesay)
print("bilesenler birlestiriliyor")
veri <- list(maddepar=maddepar,
yetenek=yetenek,
seed=seed,
cevaplar=cevaplar)
return(veri)
}
#cat("atanan seed degeri")
#veri2 <- veri_uretimi (10,1000) #(madde, tekrar, seed değeri)
#veriler <- replicate(5L, veri_uretimi (10,1000))
#lapply(1:5, function(i) veriler[,i][1])
#sapply(1:5, function(i) veriler[,i][1])
Kaç çekirdekle çalıştığı, paralel analiz yapacağız.
#library(doParallel)
#detectCores()
#cl <- makeCluster(3) #en fazla n-2 olsun
#registerDoParallel(cl)
#foreach(i=1:4)%dopar% sqrt(i)
#foreach(i=1:4, j=1:4) %dopar% { sqrt(i * j) }
#stopCluster(cl)
hızlı işlem yaptığı için tercih ediyoruz.
for(i in 1:5){
sum(rnorm(1e8))
}
#cl <- makeCluster(3)
#registerDoParallel(cl)
#foreach(i=1:5) %dopar% {sum(rnorm(1e8))}
#stopCluster(cl)
#cl <- makeCluster(8)
#registerDoParallel(cl)
tekrar=5L
maddesay=10
bireysay=1000
#simulasyon <-foreach(i=1:5,
# .packages = c("mirt", "irtoys", "doParallel"),
# .combine = rbind) %dopar% {
# adim1 <- veri_uretimi(maddesay = maddesay,
# bireysay = bireysay)
# adim2 <- kestirilen_par(adim1$cevaplar)
# hata(adim2, adim1$maddepar)
# }
##task 1 failed - "The following items have only one response category and cannot be estimated: madde1 madde3 madde4 madde5 madde7 madde9 madde10 "
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.5.3
## Warning: package 'dplyr' was built under R version 4.5.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::%+%() masks psych::%+%()
## ✖ ggplot2::alpha() masks psych::alpha()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#simulasyon_v1 <- simulasyon %>%
# group_by(parametreler) %>%
# summarise(rmse=round(mean(rmse),2),
# bias= round(mean(bias),2),
# se= round(mean(se),2)) %>%
# mutate(maddesay=maddesay,
# bireysay=bireysay) %>%
# as.data.frame()
#simulasyon_v1
#stopCluster(cl)
dlnorm’dan da üretebiliriz, rl.norm.
Bir sürü hata aldım, aslında derste anlamıştım. Bazı hatalarla baş edemeyince diyezleyerek sisteme yüklüyorum.