Simulasyon

Simülasyon neden yapılır?

Gerçek veri toplamak zaman alır, pahalıdır veya bazen veriyi toplamak mümkün değildir. Simülasyonla istediğimiz dağılımdan, istediğimiz büyüklükte veri üretebiliriz ve sonuçların gerçekte nasıl görünmesi gerektiğini bildiğimiz için yorumlamak da kolaylaşır.

R’da temel simülasyon fonksiyonları:

rnorm(n, mean, sd)- Normal runif(n, min, max) - Düzgün (Uniform) rpois(n, lambda) - Poisson rbinom(n, size, prob) - Binom

Normal dağılım

örnek olarak ortalaması 170, standart sapması 10 olan 500 tane boy uzunluğu verisi üretelim.

set.seed(1)

boy <- rnorm(n = 500, mean = 170, sd = 10)

cat("Ortalama:", round(mean(boy), 2), "\n")

## Ortalama: 170.23

cat("Standart Sapma:", round(sd(boy), 2), "\n")

## Standart Sapma: 10.12

cat("Min:", round(min(boy), 2), "| Max:", round(max(boy), 2), "\n")

## Min: 139.92 | Max: 208.1

hist(boy,
     breaks = 30,
     probability = TRUE,
     col = "#AED6F1",
     border = "white",
     main = "Normal Dağılımdan Üretilen Boy Verisi (n=500)",
     xlab = "Boy (cm)",
     ylab = "Yoğunluk")

curve(dnorm(x, mean = mean(boy), sd = sd(boy)),
      col = "#1A5276",
      lwd = 2,
      add = TRUE)

Üretilen 500 gözlemlik boy verisi, ortalama 169.74 cm ve standart sapma 9.72 cm ile normal dağılıma yakın bir yapı sergilemektedir. Histogramın simetrik görünümü ve teorik normal dağılım eğrisine olan yakınlığı, rnorm() fonksiyonunun beklenen dağılımı başarıyla ürettiğini göstermektedir. Gözlenen hafif sapmalar, örneklem kaynaklı rastlantısallığın doğal bir sonucudur. Örneklem büyüklüğü artırıldığında teorik eğriye olan uyumun iyileşmesi beklenmektedir.

Uniform Dağılım

0 ile 100 arasında eşit olasılıkla 500 not değeri üretelim.

set.seed(1)

not <- runif(n = 500, min = 0, max = 100)

cat("Ortalama:", round(mean(not), 2), "\n")

## Ortalama: 49.57

cat("Standart Sapma:", round(sd(not), 2), "\n")

## Standart Sapma: 28.33

Grafiğini çizelim

hist(not,
     breaks = 20,
     probability = TRUE,
     col = "#A9DFBF",
     border = "white",
     main = "uniform Dağılımdan Üretilen Not Verisi (n=500)",
     xlab = "Not",
     ylab = "Yoğunluk")

abline(h = 1/100, col = "#1D8348", lwd = 2, lty = 2)

legend("topright", legend = "Teorik Yoğunluk (1/100)",
       col = "#1D8348", lwd = 2, lty = 2)

500 gözlemlik not verisi, 0 ile 100 arasında düzgün (uniform) dağılımdan üretilmiştir. Ortalama 49.57 ve standart sapma 28.33 olarak hesaplanmıştır. Histogramda çubukların teorik yoğunluk değeri olan 0.01 civarında rastlantısal dalgalanmalar gösterdiği görülmektedir. Normal dağılımdan farklı olarak belirli bir değer aralığında yoğunlaşma gözlemlenmemekte, tüm değerler eşit olasılıkla üretilmektedir.

Uniform (düzgün) dağılım, diğer dağılımları simüle etmek için bir araç olarak kullanılır. Örneğin, elimizde sadece runif() olsa bile, ondan ürettiğimiz sayıları dönüştürerek normal, üstel, hatta Poisson dağılımı elde edebiliriz. Bu yüzden simülasyon dünyasında “her şeyin temeli” olarak bilinmektedir. Kısacası düzgün dağılım, “hiçbir değerin diğerinden ayrıcalıklı olmadığı” durumları modeller. (Örnek olarak, Piyango veya tombalada her sayının çekilme olasılığı eşittir)

Poisson dağılımı

Günde ortalama 5 müşteri gelen bir mağaza için 500 günlük veri üretelim.

set.seed(1)

musteri <- rpois(n = 500, lambda = 5)

cat("Ortalama:", round(mean(musteri), 2), "\n")

## Ortalama: 5.01

cat("Varyans:", round(var(musteri), 2), "\n")

## Varyans: 4.72

ortalama ve varyansın birbirine çok yakın olması poisson dağılımının en önemli özelliğidir. örneklem büyüklüğü arttıkça eşit çıkması beklenir.

barplot(table(musteri),
        col = "#F9E79F",
        border = "white",
        main = "Poisson Dagilimindan uretilen Musteri Sayısı (λ=5, n=500)",
        xlab = "Gunluk Musteri Sayısı",
        ylab = "Frekans")

500 günlük müşteri sayısı verisi, ortalama günlük müşteri sayısı λ=5 olan Poisson dağılımından üretilmiştir. Elde edilen ortalama 5.01 ve varyans 4.72 olup bu değerler Poisson dağılımının temel özelliği olan ortalama-varyans eşitliğini desteklemektedir. Grafikte en yüksek frekansın 4 ve 5 değerlerinde gözlemlendiği, sağa doğru hafif bir kuyruk yapısının bulunduğu görülmektedir.

çok değişkenli normal dağılım

library(MASS)

Ortalama Vektörü Tek değişkenlide tek bir ortalama vardı. Ve ortalama = 170 olarak ayarlanmıştı. Çok değişkenlide her değişkenin kendi ortalaması vardır. Bunları bir arada yazıyoruz:

mu <- c(170, 70)  # boy ortalaması 170, ağırlık ortalaması 70

mvrnorm ile veri üretimi

library(MASS)

set.seed(1)

# Ortalama vektörü
mu <- c(170, 70)  # boy: 170 cm, ağırlık: 70 kg

# Kovaryans matrisi
sigma <- matrix(c(100, 30,
                   30, 64), nrow = 2)

# Çok değişkenli normal dağılımdan veri üret
veri <- mvrnorm(n = 500, mu = mu, Sigma = sigma)

# Sütun isimlerini tanımla
colnames(veri) <- c("boy", "agirlik")

# İlk 6 satıra bak
head(veri)

##           boy  agirlik
## [1,] 176.1574 72.87715
## [2,] 167.2686 70.79268
## [3,] 173.8677 81.51314
## [4,] 155.0232 61.43133
## [5,] 170.2486 62.32746
## [6,] 183.1072 64.86158

# Her değişken için ortalama
cat("Boy ortalaması:", round(mean(veri[, "boy"]), 2), "\n")

## Boy ortalaması: 169.63

cat("Ağırlık ortalaması:", round(mean(veri[, "agirlik"]), 2), "\n")

## Ağırlık ortalaması: 70.15

# Kovaryans matrisi
cat("\nGözlenen kovaryans matrisi:\n")

## 
## Gözlenen kovaryans matrisi:

round(cov(veri), 2)

##            boy agirlik
## boy     106.28   27.14
## agirlik  27.14   66.07

# Korelasyon matrisi
cat("\nKorelasyon matrisi:\n")

## 
## Korelasyon matrisi:

round(cor(veri), 2)

##          boy agirlik
## boy     1.00    0.32
## agirlik 0.32    1.00

plot(veri[, "boy"], veri[, "agirlik"],
     main = "Boy ve Agirlik Iliskisi",
     xlab = "Boy (cm)",
     ylab = "Agirlik (kg)",
     col = "#2E86C1",
     pch = 16,
     cex = 0.7)

abline(lm(veri[, "agirlik"] ~ veri[, "boy"]),
       col = "#E74C3C",
       lwd = 2)

Çok değişkenli normal dağılımdan üretilen 500 gözlemlik boy ve ağırlık verisi incelendiğinde, gözlenen ortalamaların (boy: 169.63 cm, ağırlık: 70.15 kg) hedef değerlere yakın olduğu görülmektedir. Gözlenen kovaryans matrisi, tanımlanan teorik yapıyla büyük ölçüde örtüşmektedir. Saçılım grafiği, iki değişken arasında pozitif yönlü orta düzeyde bir ilişki (r = 0.32) olduğunu ortaya koymaktadır; boy arttıkça ağırlığın da artma eğiliminde olduğu gözlemlenmektedir.

Korelasyon yapısını tanımlayarak veri üretimi

Boy=ortalama 170, standart sapma 10; kilo= ortalama 70, standart sapma 8; yaş= ortalama 35 standart sapma 10 ve korelasyonlar boy-ağırlık:0,60 boy-yaş: 0,10 ağırlık-yaş: 0,30

değerlerine sahip olan veri üretelim.

library(MASS)

set.seed(42)

# Korelasyon matrisi
R <- matrix(c(1.00, 0.60, 0.10,
              0.60, 1.00, 0.30,
              0.10, 0.30, 1.00), nrow = 3)

# Standart sapmalar
sd_vektor <- c(10, 8, 10)

# Korelasyondan kovaryans matrisine dönüştür
sigma <- diag(sd_vektor) %*% R %*% diag(sd_vektor)

round(sigma, 2)

##      [,1] [,2] [,3]
## [1,]  100   48   10
## [2,]   48   64   24
## [3,]   10   24  100

veriyi üretelim

# Ortalama vektörü
mu <- c(170, 70, 35)

# Veri üret
veri3 <- mvrnorm(n = 500, mu = mu, Sigma = sigma)
colnames(veri3) <- c("boy", "agirlik", "yas")

# Gözlenen korelasyon matrisine bak
cat("Gözlenen korelasyon matrisi:\n")

## Gözlenen korelasyon matrisi:

round(cor(veri3), 2)

##          boy agirlik  yas
## boy     1.00    0.60 0.04
## agirlik 0.60    1.00 0.28
## yas     0.04    0.28 1.00

pairs(veri3,
      main = "Degiskenler Arasi Iliskiler",
      col = "#2E86C1",
      pch = 16,
      cex = 0.5)

Üç değişkenli çok değişkenli normal dağılımdan üretilen 500 gözlemlik veri incelendiğinde, tanımlanan korelasyon yapısının büyük ölçüde korunduğu görülmektedir. Boy ile ağırlık arasındaki korelasyon 0.60 ile en güçlü ilişkiyi sergilerken, ağırlık ile yaş arasında 0.28 düzeyinde orta, boy ile yaş arasında ise 0.04 gibi oldukça zayıf bir ilişki gözlemlenmiştir. Saçılım grafikleri bu örüntüyü görsel olarak da doğrulamaktadır.

Madde Tepki Kuramı (MTK) Simülasyonu

MTK, bir testteki maddelerin özelliklerini ve bireylerin yetenek düzeylerini aynı anda modellemeye çalışır. Klasik test teorisinden farkları: bir maddenin güçlük düzeyi, onu çözen gruba bağlı değildir yetenek ve madde özellikleri birbirinden bağımsız tahmin edilir

1PL Modeli

# 1PL olasılık fonksiyonu
p_1pl <- function(theta, b) {
  1 / (1 + exp(-(theta - b)))
}

# Birkaç örnek hesaplama
cat("Yetenek=0, Gucluk=0 :", round(p_1pl(0, 0), 3), "\n")   # Eşit düzey

## Yetenek=0, Gucluk=0 : 0.5

cat("Yetenek=2, Gucluk=0 :", round(p_1pl(2, 0), 3), "\n")   # Yetenekli, kolay soru

## Yetenek=2, Gucluk=0 : 0.881

cat("Yetenek=0, Gucluk=2 :", round(p_1pl(0, 2), 3), "\n")   # Zayıf, zor soru

## Yetenek=0, Gucluk=2 : 0.119

cat("Yetenek=-1, Gucluk=1 :", round(p_1pl(-1, 1), 3), "\n") # Zayıf, zor soru

## Yetenek=-1, Gucluk=1 : 0.119

✓Yetenek=0, Güçlük=0 → 0.50 → öğrencinin yeteneği ile maddenin zorluğu eşit, yani %50 olasılıkla maddeyi doğru yanıtlar. ✓Yetenek=2, Güçlük=0 → 0.881 → yetenek düzeyi yüksek öğrenci, kolay madde → %88 olasılıkla maddeyi doğru yanıtlar. ✓Yetenek=0, Güçlük=2 → 0.119 → yetenek düzeyi düşük öğrenci, zor madde → %12 olasılıkla maddeyi doğru yanıtlar ✓Yetenek=-1, Güçlük=1 → 0.119 → hem öğrenci zayıf hem madde zor → %12 olasılıkla maddeyi doğru yanıtlar

set.seed(1)

# 1000 öğrencinin yetenek düzeyi
# MTK'da yetenek genellikle ortalaması 0, sd=1 olan normal dağılımdan üretilir
theta <- rnorm(n = 1000, mean = 0, sd = 1)

cat("Ortalama:", round(mean(theta), 2), "\n")

## Ortalama: -0.01

cat("Standart Sapma:", round(sd(theta), 2), "\n")

## Standart Sapma: 1.03

Madde Karakteristik Eğrisi (MKE)

# Yetenek ekseninde -4 ile +4 arası değerler
theta_seq <- seq(-4, 4, by = 0.1)

# Üç farklı güçlükte madde
p_kolay  <- p_1pl(theta_seq, b = -1)
p_orta   <- p_1pl(theta_seq, b =  0)
p_zor    <- p_1pl(theta_seq, b =  1)

# Grafik
plot(theta_seq, p_orta,
     type = "l", lwd = 2, col = "#2E86C1",
     ylim = c(0, 1),
     main = "1PL Madde Karakteristik Egrileri",
     xlab = "Yetenek (theta)",
     ylab = "Dogru Yanitlama Olasiligi")

lines(theta_seq, p_kolay, lwd = 2, col = "#27AE60")
lines(theta_seq, p_zor,   lwd = 2, col = "#E74C3C")

abline(h = 0.5, lty = 2, col = "gray")

legend("topleft",
       legend = c("Kolay (b=-1)", "Orta (b=0)", "Zor (b=1)"),
       col = c("#27AE60", "#2E86C1", "#E74C3C"),
       lwd = 2)

1PL modeline göre çizilen madde karakteristik eğrileri incelendiğinde, üç maddenin aynı ayırt edicilik eğimine sahip olduğu, yalnızca güçlük parametresi bakımından farklılaştığı görülmektedir. Kolay madde (b=-1) olasılık eğrisi sola, zor madde (b=1) ise sağa kaymıştır. Her eğri, θ = b noktasında 0.50 olasılık değerini almaktadır.

2PL Modeli

# 2PL olasılık fonksiyonu
p_2pl <- function(theta, a, b) {
  1 / (1 + exp(-a * (theta - b)))
}

cat("a=1, b=0, theta=0 :", round(p_2pl(0, a=1, b=0), 3), "\n")

## a=1, b=0, theta=0 : 0.5

cat("a=2, b=0, theta=0 :", round(p_2pl(0, a=2, b=0), 3), "\n")

## a=2, b=0, theta=0 : 0.5

cat("a=0.5, b=0, theta=1 :", round(p_2pl(1, a=0.5, b=0), 3), "\n")

## a=0.5, b=0, theta=1 : 0.622

cat("a=2, b=0, theta=1 :", round(p_2pl(1, a=2, b=0), 3), "\n")

## a=2, b=0, theta=1 : 0.881

✓a=1, b=0, theta=0 → 0.50 ve a=2, b=0, theta=0 → 0.50 → ikisi de 0.50 çıktı çünkü θ = b olduğunda a ne olursa olsun olasılık her zaman 0.50. a parametresi eğrinin şeklini değiştiriyor, konumunu değil.

✓a=0.5, b=0, theta=1 → 0.622 → ayırt edicilik düşük, θ b’den 1 birim yüksek olmasına rağmen olasılık sadece %62

✓a=2, b=0, theta=1 → 0.881 → ayırt edicilik yüksek, aynı 1 birim fark olasılığı %88’e çıkardı

a büyüdükçe yetenek farkına daha duyarlı hale geliyor.

MKE Grafiği

theta_seq <- seq(-4, 4, by = 0.1)

# Aynı güçlük (b=0), farklı ayırt edicilik
p_dusuk  <- p_2pl(theta_seq, a = 0.5, b = 0)
p_orta   <- p_2pl(theta_seq, a = 1.0, b = 0)
p_yuksek <- p_2pl(theta_seq, a = 2.0, b = 0)

plot(theta_seq, p_orta,
     type = "l", lwd = 2, col = "#2E86C1",
     ylim = c(0, 1),
     main = "2PL Madde Karakteristik Egrileri",
     xlab = "Yetenek (theta)",
     ylab = "Dogru Yanitlama Olasiligi")

lines(theta_seq, p_dusuk,  lwd = 2, col = "#27AE60")
lines(theta_seq, p_yuksek, lwd = 2, col = "#E74C3C")

abline(h = 0.5, lty = 2, col = "gray")
abline(v = 0,   lty = 2, col = "gray")

legend("topleft",
       legend = c("Dusuk (a=0.5)", "Orta (a=1.0)", "Yuksek (a=2.0)"),
       col = c("#27AE60", "#2E86C1", "#E74C3C"),
       lwd = 2)

2PL modelinde üretilen madde karakteristik eğrileri incelendiğinde, güçlük parametresi sabit tutulurken ayırt edicilik parametresinin (a) eğrinin dikliğini doğrudan etkilediği görülmektedir. Yüksek ayırt ediciliğe sahip madde (a=2.0) yetenek farklılıklarına çok daha duyarlı bir yapı sergilemekte; düşük ayırt ediciliğe sahip madde (a=0.5) ise yetenek düzeyinden bağımsız biçimde düz bir seyir izlemektedir. Tüm eğriler b=0 noktasında 0.50 olasılık değerinde kesişmektedir.

3PL Modeli

# 3PL olasılık fonksiyonu
p_3pl <- function(theta, a, b, c) {
  c + (1 - c) * (1 / (1 + exp(-a * (theta - b))))
}

# Test edelim
cat("Cok dusuk yetenek (theta=-4):\n")

## Cok dusuk yetenek (theta=-4):

cat("1PL:", round(p_1pl(-4, b=0), 3), "\n")

## 1PL: 0.018

cat("2PL:", round(p_2pl(-4, a=1, b=0), 3), "\n")

## 2PL: 0.018

cat("3PL:", round(p_3pl(-4, a=1, b=0, c=0.25), 3), "\n")

## 3PL: 0.263

1PL ve 2PL modelde p değeri 0.018 bulunmuştur. çok düşük yetenekli öğrenci için maddeyi doğru yanıtlama olasılıkğı neredeyse sıfır bulunmuştur. 3PL modelde ise p değeri 0.263 olarak elde edilmiştir. Aynı öğrenci şans parametresi modele dahil edildiğinde (c=0.25) ile %26 olasılığa sahiptir. yani şans eseri doğru yanıtlayabilir.