[過去ログ] 高校数学の質問スレ(医者・東大卒専用) Part438 (1002レス)
上下前次1-新
抽出解除 必死チェッカー(本家) (べ) 自ID レス栞 あぼーん
このスレッドは過去ログ倉庫に格納されています。
次スレ検索 歴削→次スレ 栞削→次スレ 過去ログメニュー
リロード規制です。10分ほどで解除するので、他のブラウザへ避難してください。
861: 06/02(月)11:14 ID:GMuHFUYr(1/2) AAS
 x = c(-0.86, -0.3, -0.05, 0.73) 
 n = c(5, 5, 5, 5) 
 y = c(0, 1, 3, 5) 
 (data = cbind(x, n, y)) 
 (response = cbind(y, n - y) ) 
 results = glm(response ~ x, family = binomial) 
 #summary(results) 
  
 -results$coef[1]/results$coef[2] 
 library(MASS)  # mvrnorm を使うため 
  
 # 推定された係数と共分散行列 
 beta_hat = coef(results) 
 (vcov_matrix = vcov(results)) 
  
 # 多変量正規乱数を生成(β0, β1) 
 set.seed(123)  # 再現性のため 
 samples = mvrnorm(n = 10000, mu = beta_hat, Sigma = vcov_matrix) 
  
 # 各サンプルから LD50 を計算 
 LD50_samples = -samples[,1] / samples[,2] 
  
 # 信頼区間(95%) 
 CI = quantile(LD50_samples, probs = c(0.025, 0.975)) 
  
 # 結果表示 
 cat("シミュレーションによるLD50の95%信頼区間:\n") 
 print(CI) 
862: 06/02(月)12:04 ID:GMuHFUYr(2/2) AAS
 # データ 
 x <- c(-0.86, -0.3, -0.05, 0.73) 
 n <- c(5, 5, 5, 5) 
 y <- c(0, 1, 3, 5) 
  
 # JAGSモデル 
 model_string <- " 
 model { 
   for (i in 1:N) { 
     y[i] ~ dbin(p[i], n[i]) 
     logit(p[i]) <- beta0 + beta1 * x[i] 
   } 
   # 事前分布(非情報的) 
   beta0 ~ dnorm(0.0, 0.001) 
   beta1 ~ dnorm(0.0, 0.001) 
  
   # LD50の定義 
   LD50 <- -beta0 / beta1 
 } 
 " 
  
 # JAGSに渡すデータ 
 data_jags <- list( 
   x = x, 
   n = n, 
   y = y, 
   N = length(y) 
 ) 
  
 # 初期値 
 inits <- function() { 
   list(beta0 = rnorm(1, 0, 1), beta1 = rnorm(1, 0, 1)) 
 } 
  
 # モデル作成と実行 
 model <- jags.model(textConnection(model_string), data = data_jags, inits = inits, n.chains = 3) 
 update(model, 1000)  # バーンイン 
  
 # サンプリング 
 samples <- coda.samples(model, variable.names = c("beta0", "beta1", "LD50"), n.iter = 10000) 
  
 # 結果表示(LD50の95%信用区間) 
 summary(samples) 
 LD50_samples <- as.matrix(samples)[, "LD50"] 
 quantile(LD50_samples, probs = c(0.025, 0.975)) 
 plot(samples) 
上下前次1-新書関写板覧索設栞歴
スレ情報 赤レス抽出 画像レス抽出 歴の未読スレ
ぬこの手 ぬこTOP 0.029s