- こちらで多次元データを扱った。
- そこでは、k次元立方格子状のデータ格納をarray()を使って行うとともに
- それをハンドリングするためのn進数ハンドリングを考えた。
- 多項分布では項数-1の次元で格納することができるはず・・・
- お試し中ソース
k<-3 # カテゴリ数
# 自由変数の数はk-1
n<-10 # 総数
library(MCMCpack)
prob<-rdirichlet(1,rep(1,k)) # カテゴリの生起確率
aaa<-array(0,rep((n+1),(k-1))) # 計算結果格納用アレイ
# (n+1)^自由変数の数 の値がとりうる値
keta<-0:(k-2) # n進数計算のための道具
keta<-(n+1)^keta # n進数のアレイ内番地計算用の道具
counter<-rep(0,k-1) # 同上
counter[1]<--1
loop<-TRUE
while(loop){
counter[1]<-counter[1]+1 #順々にn進数をその一桁目から増やす
# くらい上がりの処理
for(i in 1:(k-2)){
if(counter[i]==n+2){
counter[i]<-0
counter[i+1]<-counter[i+1]+1
}
}
# 1列のベクトルで言うと、これは何番目?
banchi<-sum(keta*counter)+1
vvv<-counter+1
# vvvはそれぞれの自由変数の値
# 不自由な変数の値は
lastval<-n-sum(vvv)
vvv2<-sum(vvv)-k+1
tmpval<-c(vvv,vvv2)
tmpprob<-prob^tmpval
aaa[banchi]<-prod(tmpprob)*exp(lgamma(n+1)-sum(lgamma(tmpval+1)))
if(sum(counter)==(n+1)*(k-1)){
loop<-FALSE
}
}